From: bangerth Date: Tue, 12 Feb 2013 15:20:15 +0000 (+0000) Subject: Follow the documented advice of not looping over all elements of a matrix but do... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=502b1c06459077779a38da8e9cfaf54240e10ef6;p=dealii-svn.git Follow the documented advice of not looping over all elements of a matrix but do so row by row. git-svn-id: https://svn.dealii.org/trunk@28325 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/source/lac/sparse_direct.cc b/deal.II/source/lac/sparse_direct.cc index 4da14da2ad..eb42a4d19e 100644 --- a/deal.II/source/lac/sparse_direct.cc +++ b/deal.II/source/lac/sparse_direct.cc @@ -1726,16 +1726,21 @@ factorize (const Matrix &matrix) // not yet written to std::vector row_pointers = Ap; - for (typename Matrix::const_iterator p=matrix.begin(); - p!=matrix.end(); ++p) + // loop over the elements of the matrix row by row, as suggested + // in the documentation of the sparse matrix iterator class + for (unsigned int row = 0; row < matrix.n(); ++row) { - // write entry into the first - // free one for this row - Ai[row_pointers[p->row()]] = p->column(); - Ax[row_pointers[p->row()]] = p->value(); + for (typename Matrix::const_iterator p=matrix.begin(row); + p!=matrix.end(row); ++p) + { + // write entry into the first + // free one for this row + Ai[row_pointers[row]] = p->column(); + Ax[row_pointers[row]] = p->value(); - // then move pointer ahead - ++row_pointers[p->row()]; + // then move pointer ahead + ++row_pointers[row]; + } } // at the end, we should have @@ -1955,15 +1960,20 @@ void SparseDirectMUMPS::initialize_matrix (const Matrix &matrix) unsigned int index = 0; - for (typename Matrix::const_iterator ptr = matrix.begin (); - ptr != matrix.end (); ++ptr) - if (std::abs (ptr->value ()) > 0.0) - { - a[index] = ptr->value (); - irn[index] = ptr->row () + 1; - jcn[index] = ptr->column () + 1; - ++index; - } + // loop over the elements of the matrix row by row, as suggested + // in the documentation of the sparse matrix iterator class + for (unsigned int row = 0; row < matrix.n(); ++row) + { + for (typename Matrix::const_iterator ptr = matrix.begin (row); + ptr != matrix.end (row); ++ptr) + if (std::abs (ptr->value ()) > 0.0) + { + a[index] = ptr->value (); + irn[index] = row + 1; + jcn[index] = ptr->column () + 1; + ++index; + } + } id.n = n; id.nz = nz;