From: kronbichler Date: Tue, 21 Apr 2009 08:02:28 +0000 (+0000) Subject: Provide a faster implementation of precondition_SSOR. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c880b9bcb97685dbe45b09a7a92d9ae8401c2110;p=dealii-svn.git Provide a faster implementation of precondition_SSOR. git-svn-id: https://svn.dealii.org/trunk@18670 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/lac/include/lac/precondition.h b/deal.II/lac/include/lac/precondition.h index 57227b9a59..d0194995c4 100644 --- a/deal.II/lac/include/lac/precondition.h +++ b/deal.II/lac/include/lac/precondition.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006 by the deal.II authors +// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2009 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -459,6 +459,26 @@ template > class PreconditionSSOR : public PreconditionRelaxation { public: + + /** + * A typedef to the base class. + */ + typedef PreconditionRelaxation BaseClass; + + + /** + * Initialize matrix and + * relaxation parameter. The + * matrix is just stored in the + * preconditioner object. The + * relaxation parameter should be + * larger than zero and smaller + * than 2 for numerical + * reasons. It defaults to 1. + */ + void initialize (const MATRIX &A, + typename BaseClass::AdditionalData parameters = BaseClass::AdditionalData()); + /** * Apply preconditioner. */ @@ -474,6 +494,14 @@ class PreconditionSSOR : public PreconditionRelaxation */ template void Tvmult (VECTOR&, const VECTOR&) const; + + private: + /** + * An array that stores for each matrix + * row where the first position after + * the diagonal is located. + */ + std::vector pos_right_of_diagonal; }; @@ -893,13 +921,46 @@ PreconditionSOR::Tvmult (VECTOR &dst, const VECTOR &src) const //--------------------------------------------------------------------------- +template +inline void +PreconditionSSOR::initialize (const MATRIX &rA, + typename BaseClass::AdditionalData parameters) +{ + this->PreconditionRelaxation::initialize (rA, parameters); + + // calculate the positions first after + // the diagonal. + const std::size_t * rowstart_ptr = + this->A->get_sparsity_pattern().get_rowstart_indices(); + const unsigned int * const colnums = + this->A->get_sparsity_pattern().get_column_numbers(); + const unsigned int n = this->A->n(); + pos_right_of_diagonal.resize(n); + for (unsigned int row=0; row template inline void PreconditionSSOR::vmult (VECTOR &dst, const VECTOR &src) const { Assert (this->A!=0, ExcNotInitialized()); - this->A->precondition_SSOR (dst, src, this->relaxation); + this->A->precondition_SSOR (dst, src, this->relaxation, pos_right_of_diagonal); } diff --git a/deal.II/lac/include/lac/sparse_decomposition.templates.h b/deal.II/lac/include/lac/sparse_decomposition.templates.h index 42b0b536fa..3d58ed0ae9 100644 --- a/deal.II/lac/include/lac/sparse_decomposition.templates.h +++ b/deal.II/lac/include/lac/sparse_decomposition.templates.h @@ -1,5 +1,5 @@ //--------------------------------------------------------------------------- -// Copyright (C) 2002, 2003, 2004, 2005, 2006 by the deal.II authors +// Copyright (C) 2002, 2003, 2004, 2005, 2006, 2009 by the deal.II authors // by the deal.II authors and Stephen "Cheffo" Kolaroff // // This file is subject to QPL and may not be distributed diff --git a/deal.II/lac/include/lac/sparse_ilu.h b/deal.II/lac/include/lac/sparse_ilu.h index b5938cbcec..8c4f2ceaf3 100644 --- a/deal.II/lac/include/lac/sparse_ilu.h +++ b/deal.II/lac/include/lac/sparse_ilu.h @@ -1,5 +1,5 @@ //--------------------------------------------------------------------------- -// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2008 by the deal.II authors +// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2008, 2009 by the deal.II authors // by the deal.II authors and Stephen "Cheffo" Kolaroff // // This file is subject to QPL and may not be distributed @@ -46,7 +46,7 @@ DEAL_II_NAMESPACE_OPEN * @; others can be generated in application programs (see the * section on @ref Instantiations in the manual). * - * @author Wolfgang Bangerth, 2008; unified interface: Ralf Hartmann + * @author Wolfgang Bangerth, 2008, 2009; unified interface: Ralf Hartmann */ template class SparseILU : public SparseLUDecomposition diff --git a/deal.II/lac/include/lac/sparse_ilu.templates.h b/deal.II/lac/include/lac/sparse_ilu.templates.h index c19807c119..01b4a7164d 100644 --- a/deal.II/lac/include/lac/sparse_ilu.templates.h +++ b/deal.II/lac/include/lac/sparse_ilu.templates.h @@ -1,5 +1,5 @@ //--------------------------------------------------------------------------- -// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2008 by the deal.II authors +// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2008, 2009 by the deal.II authors // by the deal.II authors and Stephen "Cheffo" Kolaroff // // This file is subject to QPL and may not be distributed diff --git a/deal.II/lac/include/lac/sparse_matrix.h b/deal.II/lac/include/lac/sparse_matrix.h index b6d3480613..f3948e2ee0 100644 --- a/deal.II/lac/include/lac/sparse_matrix.h +++ b/deal.II/lac/include/lac/sparse_matrix.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name: $ // -// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors +// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2009 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -714,17 +714,20 @@ class SparseMatrix : public virtual Subscriptor /** * Return the number of actually - * nonzero elements of this - * matrix. + * nonzero elements of this matrix. It + * is possible to specify the parameter + * threshold in order to count + * only the elements that have absolute + * value greater than the threshold. * - * Note, that this function does - * (in contrary to - * n_nonzero_elements()) not - * count all entries of the - * sparsity pattern but only the - * ones that are nonzero. - */ - unsigned int n_actually_nonzero_elements () const; + * Note, that this function does (in + * contrary to n_nonzero_elements()) + * not count all entries of the + * sparsity pattern but only the ones + * that are nonzero (or whose absolute + * value is greater than threshold). + */ + unsigned int n_actually_nonzero_elements (const double threshold = 0.) const; /** * Return a (constant) reference @@ -1472,9 +1475,10 @@ class SparseMatrix : public virtual Subscriptor * src. */ template - void precondition_SSOR (Vector &dst, - const Vector &src, - const number om = 1.) const; + void precondition_SSOR (Vector &dst, + const Vector &src, + const number om = 1., + const std::vector&pos_right_of_diagonal=std::vector()) const; /** * Apply SOR preconditioning diff --git a/deal.II/lac/include/lac/sparse_matrix.templates.h b/deal.II/lac/include/lac/sparse_matrix.templates.h index d81775ef71..a9e7711739 100644 --- a/deal.II/lac/include/lac/sparse_matrix.templates.h +++ b/deal.II/lac/include/lac/sparse_matrix.templates.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors +// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2009 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -1155,9 +1155,10 @@ SparseMatrix::precondition_Jacobi (Vector &dst, template template void -SparseMatrix::precondition_SSOR (Vector &dst, - const Vector &src, - const number om) const +SparseMatrix::precondition_SSOR (Vector &dst, + const Vector &src, + const number om, + const std::vector &pos_right_of_diagonal) const { // to understand how this function works // you may want to take a look at the CVS @@ -1175,6 +1176,55 @@ SparseMatrix::precondition_SSOR (Vector &dst, const std::size_t *rowstart_ptr = &cols->rowstart[0]; somenumber *dst_ptr = &dst(0); + // case when we have stored the position + // just right of the diagonal (then we + // don't have to search for it). + if (pos_right_of_diagonal.size() != 0) + { + Assert (pos_right_of_diagonal.size() == dst.size(), + ExcDimensionMismatch (pos_right_of_diagonal.size(), dst.size())); + + // forward sweep + for (unsigned int row=0; rowcolnums[j]); + + // divide by diagonal element + *dst_ptr = s * om / val[*rowstart_ptr]; + }; + + 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) + { + const unsigned int end_row = *(rowstart_ptr+1); + const unsigned int first_right_of_diagonal_index + = pos_right_of_diagonal[row]; + double s = *dst_ptr; + for (unsigned int j=first_right_of_diagonal_index; jcolnums[j]); + + *dst_ptr = s * om / val[*rowstart_ptr]; + }; + return; + } + + // case when we need to get the position + // of the first element right of the + // diagonal manually for each sweep. // forward sweep for (unsigned int row=0; row:: precondition_SSOR (Vector &, const Vector &, - const S1) const; + const S1, + const std::vector&) const; template void SparseMatrix:: precondition_SOR (Vector &,