From: guido Date: Wed, 23 Apr 2003 08:13:15 +0000 (+0000) Subject: new permuted SOR, must be tested X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=26bf56085845f8564a0c44f0bb748290115eae02;p=dealii-svn.git new permuted SOR, must be tested git-svn-id: https://svn.dealii.org/trunk@7446 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/lac/include/lac/precondition.h b/deal.II/lac/include/lac/precondition.h index 62e6652285..2e2d88b862 100644 --- a/deal.II/lac/include/lac/precondition.h +++ b/deal.II/lac/include/lac/precondition.h @@ -355,6 +355,87 @@ class PreconditionSSOR : public PreconditionRelaxation }; +/** + * Permuted SOR preconditioner using matrix built-in function. The MATRIX + * class used is required to have functions + * @p{PSOR(VECTOR&, const VECTOR&, double)} and + * @p{TPSOR(VECTOR&, const VECTOR&, double)}. + * + * + * @sect2{Usage example} + * @begin{itemize} + * // Declare related objects + * + * SparseMatrix A; + * Vector x; + * Vector b; + * SolverCG<> solver(...); + * + * //...initialize and build A + * + * std::vector permutation(x.size()); + * + * //...fill permutation with reasonable values + * + * // Define and initialize preconditioner + * + * PreconditionPSOR > precondition; + * precondition.initialize (A, permutation, .6); + * + * solver.solve (A, x, b, precondition); + * @end{itemize} + * + * @author Guido Kanschat, 2003 + */ +template > +class PreconditionPSOR : public PreconditionRelaxation +{ + public: + /** + * Initialize matrix and + * relaxation parameter. The + * matrix is just stored in the + * preconditioner object. + * + * The permutation vector is + * stored as a + * pointer. Therefore, it has to + * be assured that the lifetime + * of the vector exceeds the + * lifetime of the + * preconditioner. + * + * The relaxation parameter + * should be larger than zero and + * smaller than 2 for numerical + * reasons. It defaults to 1. + */ + void initialize (const MATRIX& A, + const std::vector& permutation, + typename PreconditionRelaxation::AdditionalData + parameters = AdditionalData()); + + /** + * Apply preconditioner. + */ + template + void vmult (VECTOR&, const VECTOR&) const; + + /** + * Apply transpose + * preconditioner. + */ + template + void Tvmult (VECTOR&, const VECTOR&) const; + private: + /** + * Storage for the permutation vector. + */ + const std::vector* permutation; +}; + + + /** * Preconditioner using an iterative solver. This preconditioner uses * a fully initialized LAC iterative solver for the approximate @@ -624,6 +705,43 @@ PreconditionSSOR::Tvmult (VECTOR& dst, const VECTOR& src) const } +//----------------------------------------------------------------------// + +template +inline void +PreconditionPSOR::initialize ( + const MATRIX &rA, + const std::vector& p, + typename PreconditionRelaxation::AdditionalData parameters) +{ + permutation = &p; + PreconditionRelaxation::initalize(rA, parameters); +} + + +template +template +inline void +PreconditionPSOR::vmult (VECTOR& dst, const VECTOR& src) const +{ + Assert (this->A!=0, ExcNotInitialized()); + dst = src; + this->A->PSOR (dst, *permutation, this->relaxation); +} + + + +template +template +inline void +PreconditionPSOR::Tvmult (VECTOR& dst, const VECTOR& src) const +{ + Assert (this->A!=0, ExcNotInitialized()); + dst = src; + this->A->TPSOR (dst, *permutation, this->relaxation); +} + + //----------------------------------------------------------------------// diff --git a/deal.II/lac/include/lac/sparse_matrix.h b/deal.II/lac/include/lac/sparse_matrix.h index d5d6680460..e69fc38545 100644 --- a/deal.II/lac/include/lac/sparse_matrix.h +++ b/deal.II/lac/include/lac/sparse_matrix.h @@ -38,7 +38,7 @@ template class FullMatrix; * corresponding ``.templates.h'' file and instantiate the respective * class yourself. * - * @author Original version by Roland Becker, Guido Kanschat, Franz-Theo Suttmeier; lots of enhancements, reorganisation and documentation by Wolfgang Bangerth 1998 + * @author several, 1994-2003 */ template class SparseMatrix : public Subscriptor @@ -802,6 +802,46 @@ class SparseMatrix : public Subscriptor */ template void TSOR (Vector &v, + const number om = 1.) const; + + /** + * Perform a permuted SOR + * preconditioning in-place. + * + * The standard SOR method is + * applied in the order + * prescribed by @p{permutation}, + * that is, first the row + * @p{permutation[0]}, then + * @p{permutation[1]} and so + * on. + * + * @p{omega} is the relaxation + * parameter. + */ + template + void PSOR (Vector &v, + const std::vector permutation, + const number om = 1.) const; + + /** + * Perform a transposed permuted SOR + * preconditioning in-place. + * + * The transposed SOR method is + * applied in the order + * prescribed by @p{permutation}, + * that is, first the row + * @p{permutation[m()-1]}, then + * @p{permutation[m()-2]} and so + * on. + * + * @p{omega} is the relaxation + * parameter. + */ + template + void TPSOR (Vector &v, + const std::vector permutation, const number om = 1.) const; /** diff --git a/deal.II/lac/include/lac/sparse_matrix.templates.h b/deal.II/lac/include/lac/sparse_matrix.templates.h index 3464c7da3e..f5cf298da6 100644 --- a/deal.II/lac/include/lac/sparse_matrix.templates.h +++ b/deal.II/lac/include/lac/sparse_matrix.templates.h @@ -993,6 +993,69 @@ SparseMatrix::TSOR (Vector& dst, const number om) const } +template +template +void +SparseMatrix::PSOR ( + Vector& dst, + const std::vector permutation, + const number om) const +{ + Assert (cols != 0, ExcMatrixNotInitialized()); + Assert (val != 0, ExcMatrixNotInitialized()); + Assert (m() == n(), ExcMatrixNotSquare()); + Assert (m() == dst.size(), ExcDimensionMismatch(m(),dst.size())); + Assert (m() == permutation.size(), + ExcDimensionMismatch(m(), permutation.size())); + + for (unsigned int urow=0; urowrowstart[row]; jrowstart[row+1]; ++j) + { + const unsigned int col = permutation[cols->colnums[j]]; + if (col < row) + s -= val[j] * dst(col); + } + + dst(row) = s * om / val[cols->rowstart[row]]; + } +} + + +template +template +void +SparseMatrix::TPSOR ( + Vector& dst, + const std::vector permutation, + const number om) const +{ + Assert (cols != 0, ExcMatrixNotInitialized()); + Assert (val != 0, ExcMatrixNotInitialized()); + Assert (m() == n(), ExcMatrixNotSquare()); + Assert (m() == dst.size(), ExcDimensionMismatch(m(),dst.size())); + Assert (m() == permutation.size(), + ExcDimensionMismatch(m(), permutation.size())); + + for (unsigned int urow=m(); urow != 0;) + { + --urow; + const unsigned int row = permutation[urow]; + somenumber s = dst(row); + for (unsigned int j=cols->rowstart[row]; jrowstart[row+1]; ++j) + { + const unsigned int col = permutation[cols->colnums[j]]; + if (col > row) + s -= val[j] * dst(col); + } + + dst(row) = s * om / val[cols->rowstart[row]]; + } +} + + template template void