From 87277ad833c18356d100f22a08729c5bc47b43db Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Mon, 20 Apr 2009 18:08:55 +0000 Subject: [PATCH] The SSOR kernel was implemented inefficiently. Improved it a bit by avoiding 'if' statements at the innermost loop and by not wasting floating point multiplications. Could still do better if we saved the position of the diagonal in an extra array at creation time of a preconditioner. git-svn-id: https://svn.dealii.org/trunk@18658 0785d39b-7218-0410-832d-ea1e28bc413d --- .../lac/include/lac/sparse_matrix.templates.h | 30 +++++++++++-------- 1 file changed, 18 insertions(+), 12 deletions(-) diff --git a/deal.II/lac/include/lac/sparse_matrix.templates.h b/deal.II/lac/include/lac/sparse_matrix.templates.h index 70e32b125e..d81775ef71 100644 --- a/deal.II/lac/include/lac/sparse_matrix.templates.h +++ b/deal.II/lac/include/lac/sparse_matrix.templates.h @@ -224,11 +224,15 @@ SparseMatrix::n_nonzero_elements () const template unsigned int -SparseMatrix::n_actually_nonzero_elements () const +SparseMatrix::n_actually_nonzero_elements (const double threshold) const { Assert (cols != 0, ExcNotInitialized()); - return std::count_if(&val[0], &val[n_nonzero_elements ()], - std::bind2nd(std::not_equal_to(), 0)); + Assert (threshold >= 0, ExcMessage ("Negative threshold!")); + unsigned int nnz = 0; + for (unsigned int i=0; i threshold) + ++nnz; + return nnz; } @@ -1188,12 +1192,13 @@ SparseMatrix::precondition_SSOR (Vector &dst, row) - &cols->colnums[0]); - + + double s = *dst_ptr; for (unsigned int j=(*rowstart_ptr)+1; jcolnums[j]); + s -= val[j] * dst(cols->colnums[j]); // divide by diagonal element - *dst_ptr /= val[*rowstart_ptr]; + *dst_ptr = s * om / val[*rowstart_ptr]; }; rowstart_ptr = &cols->rowstart[0]; @@ -1206,16 +1211,17 @@ SparseMatrix::precondition_SSOR (Vector &dst, dst_ptr = &dst(n-1); for (int row=n-1; row>=0; --row, --rowstart_ptr, --dst_ptr) { + const unsigned int end_row = *(rowstart_ptr+1); const unsigned int first_right_of_diagonal_index = (internals::SparsityPatternTools::optimized_lower_bound (&cols->colnums[*rowstart_ptr+1], - &cols->colnums[*(rowstart_ptr+1)], + &cols->colnums[end_row], static_cast(row)) - &cols->colnums[0]); - for (unsigned int j=first_right_of_diagonal_index; j<*(rowstart_ptr+1); ++j) - if (cols->colnums[j] > static_cast(row)) - *dst_ptr -= om* val[j] * dst(cols->colnums[j]); + double s = *dst_ptr; + for (unsigned int j=first_right_of_diagonal_index; jcolnums[j]); - *dst_ptr /= val[*rowstart_ptr]; + *dst_ptr = s * om / val[*rowstart_ptr]; }; } @@ -1276,7 +1282,7 @@ SparseMatrix::SOR (Vector& dst, { const unsigned int col = cols->colnums[j]; if (col < row) - s -= val[j] * dst(col); + s -= val[j] * dst(col); } dst(row) = s * om / val[cols->rowstart[row]]; -- 2.39.5