From 7aca368b123b3de3eacae66f9b5b954ba31761d5 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Tue, 7 Jul 1998 14:13:59 +0000 Subject: [PATCH] Implement a version of SSOR which is 40 per cent faster than the old version. git-svn-id: https://svn.dealii.org/trunk@438 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/lac/source/dsmatrix.cc | 60 +++++++++++++++++++++++++--------- 1 file changed, 44 insertions(+), 16 deletions(-) diff --git a/deal.II/lac/source/dsmatrix.cc b/deal.II/lac/source/dsmatrix.cc index 1b9df3fadc..100f9f3a9e 100644 --- a/deal.II/lac/source/dsmatrix.cc +++ b/deal.II/lac/source/dsmatrix.cc @@ -586,31 +586,59 @@ void dSMatrix::precondition_SSOR (dVector& dst, const dVector& src, const double om) const { + // to understand how this function works + // you may want to take a look at the CVS + // archives to see the original version + // which is much clearer... Assert (cols != 0, ExcMatrixNotInitialized()); Assert (val != 0, ExcMatrixNotInitialized()); Assert (m() == n(), ExcMatrixNotSquare()); - const unsigned int n = src.size(); - unsigned int j; + const unsigned int n = src.size(); + const unsigned int *rowstart_ptr = &cols->rowstart[0]; + double *dst_ptr = &dst(0); - for (unsigned int row=0; rowrowstart[row]; jrowstart[row+1] ;j++) - if ((unsigned int)cols->colnums[j] < row) - dst(row) -= om* val[j] * dst(cols->colnums[j]); - dst(row) /= val[cols->rowstart[row]]; - } - for (unsigned int row=0; rowrowstart[row]]; + *dst_ptr = src(row); + // find the first element in this line + // which is on the right of the diagonal. + // we need to precondition with the + // elements on the left only. + // note: the first entry in each + // line denotes the diagonal element, + // which we need not check. + const unsigned int first_right_of_diagonal_index + = (lower_bound (&cols->colnums[*rowstart_ptr+1], + &cols->colnums[*(rowstart_ptr+1)], + static_cast(row)) - + &cols->colnums[0]); + + for (unsigned int j=(*rowstart_ptr)+1; jcolnums[j]); + *dst_ptr /= val[*rowstart_ptr]; + }; - for (int row=n-1; row>=0; --row) + rowstart_ptr = &cols->rowstart[0]; + dst_ptr = &dst(0); + for (unsigned int row=0; rowrowstart[n-1]; + dst_ptr = &dst(n-1); + for (int row=n-1; row>=0; --row, --rowstart_ptr, --dst_ptr) { - for (j=cols->rowstart[row];jrowstart[row+1];j++) + const unsigned int first_right_of_diagonal_index + = (lower_bound (&cols->colnums[*rowstart_ptr+1], + &cols->colnums[*(rowstart_ptr+1)], + static_cast(row)) - + &cols->colnums[0]); + for (unsigned int j=first_right_of_diagonal_index; j<*(rowstart_ptr+1); ++j) if (cols->colnums[j] > row) - dst(row) -= om* val[j] * dst(cols->colnums[j]); - dst(row) /= val[cols->rowstart[row]]; - } + *dst_ptr -= om* val[j] * dst(cols->colnums[j]); + + *dst_ptr /= val[*rowstart_ptr]; + }; } void -- 2.39.5