From: Guido Kanschat Date: Fri, 16 Jun 2000 20:03:23 +0000 (+0000) Subject: new preconditioner interface X-Git-Tag: v8.0.0~20383 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=e2053739578d156adc9c8757f7c0103141f520aa;p=dealii.git new preconditioner interface git-svn-id: https://svn.dealii.org/trunk@3031 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/multigrid/mg_base.h b/deal.II/deal.II/include/multigrid/mg_base.h index b53705d6b5..0037d7e679 100644 --- a/deal.II/deal.II/include/multigrid/mg_base.h +++ b/deal.II/deal.II/include/multigrid/mg_base.h @@ -443,8 +443,8 @@ class PreconditionMG * This is the operator used by * LAC iterative solvers. */ - void operator() (VECTOR &dst, - const VECTOR &src) const; + void vmult (VECTOR &dst, + const VECTOR &src) const; private: /** @@ -583,8 +583,8 @@ PreconditionMG::PreconditionMG(MG &mg, template void -PreconditionMG::operator() (VECTOR &dst, - const VECTOR &src) const +PreconditionMG::vmult (VECTOR &dst, + const VECTOR &src) const { multigrid->copy_to_mg(src); multigrid->vcycle(*pre, *post, *coarse); diff --git a/deal.II/deal.II/source/numerics/base.cc b/deal.II/deal.II/source/numerics/base.cc index 523441d82c..b3281538f3 100644 --- a/deal.II/deal.II/source/numerics/base.cc +++ b/deal.II/deal.II/source/numerics/base.cc @@ -154,10 +154,8 @@ void ProblemBase::solve () PrimitiveVectorMemory<> memory; SolverCG<> cg(control,memory); - PreconditionRelaxation<> - prec(system_matrix, - &SparseMatrix::template precondition_SSOR, - 1.2); + PreconditionSSOR<> prec; + prec.initialize (system_matrix, 1.2); // solve cg.solve (system_matrix, solution, right_hand_side, prec); diff --git a/deal.II/deal.II/source/numerics/vectors.cc b/deal.II/deal.II/source/numerics/vectors.cc index f5728bf737..2d2062a377 100644 --- a/deal.II/deal.II/source/numerics/vectors.cc +++ b/deal.II/deal.II/source/numerics/vectors.cc @@ -444,10 +444,8 @@ void VectorTools::project (const DoFHandler &dof, PrimitiveVectorMemory<> memory; SolverCG<> cg(control,memory); - PreconditionRelaxation<> - prec(mass_matrix, - &SparseMatrix::template precondition_SSOR, - 1.2); + PreconditionSSOR<> prec; + prec.initialize(mass_matrix, 1.2); // solve cg.solve (mass_matrix, vec, tmp, prec); @@ -763,10 +761,8 @@ VectorTools::project_boundary_values (const DoFHandler &dof, PrimitiveVectorMemory<> memory; SolverCG<> cg(control,memory); - PreconditionRelaxation<> - prec(mass_matrix, - &SparseMatrix::template precondition_SSOR, - 1.2); + PreconditionSSOR<> prec; + prec.initialize(mass_matrix, 1.2); // solve cg.solve (mass_matrix, boundary_projection, rhs, prec); diff --git a/deal.II/lac/include/lac/eigen.h b/deal.II/lac/include/lac/eigen.h index 6293e80863..5e1ef124ae 100644 --- a/deal.II/lac/include/lac/eigen.h +++ b/deal.II/lac/include/lac/eigen.h @@ -357,6 +357,7 @@ EigenPower::solve (double &value, } memory.free(Vy); + memory.free(Vr); deallog.pop(); // Output @@ -465,7 +466,8 @@ EigenInverse::solve (double &value, } memory.free(Vy); - + memory.free(Vr); + deallog.pop(); // Output diff --git a/deal.II/lac/include/lac/full_matrix.h b/deal.II/lac/include/lac/full_matrix.h index e44678438a..0a4a9bf4b6 100644 --- a/deal.II/lac/include/lac/full_matrix.h +++ b/deal.II/lac/include/lac/full_matrix.h @@ -560,13 +560,6 @@ class FullMatrix : public Subscriptor /** * Exception */ - DeclException2 (ExcInvalidIndex, - int, int, - << "The given index " << arg1 - << " should be less than " << arg2 << "."); - /** - * Exception - */ DeclException2 (ExcDimensionMismatch, int, int, << "The two dimensions " << arg1 << " and " << arg2 @@ -699,8 +692,8 @@ template inline number FullMatrix::operator() (const unsigned int i, const unsigned int j) const { - Assert (i inline number & FullMatrix::operator() (const unsigned int i, const unsigned int j) { - Assert (i #include - template FullMatrix::FullMatrix (const unsigned int n) : val (0), diff --git a/deal.II/lac/include/lac/precondition.h b/deal.II/lac/include/lac/precondition.h index 97e65c8020..410465e54a 100644 --- a/deal.II/lac/include/lac/precondition.h +++ b/deal.II/lac/include/lac/precondition.h @@ -47,60 +47,6 @@ class PreconditionIdentity */ template void Tvmult (VECTOR&, const VECTOR&) const; - /** - * Apply preconditioner. - */ - template - void operator () (VECTOR&, const VECTOR&) const; -}; - - -/** - * Jacobi preconditioner using matrix built-in function. The MATRIX - * class used is required to have a function - * @p{precondition_Jacobi(VECTOR&, const VECTOR&, double} - * - * @author Guido Kanschat, 2000 - */ -template > -class PreconditionJacobi -{ - public: - /** - * 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, const double omega = 1.); - /** - * Apply preconditioner. - */ - template - void vmult (VECTOR&, const VECTOR&) const; - /** - * Apply transpose - * preconditioner. Since this is - * a symmetric preconditioner, - * this function is the same as - * @ref{vmult}. - */ - template - void Tvmult (VECTOR&, const VECTOR&) const; - - private: - /** - * Pointer to the matrix object. - */ - SmartPointer A; - /** - * Relaxation parameter. - */ - double omega; }; @@ -174,8 +120,8 @@ class PreconditionUseMatrix * of this object with the two * arguments given here. */ - void operator() (VECTOR &dst, - const VECTOR &src) const; + void vmult (VECTOR &dst, + const VECTOR &src) const; private: /** @@ -192,94 +138,120 @@ class PreconditionUseMatrix /** - * Preconditioner for builtin relaxation methods. - * Application of this preconditioner includes - * use of the #precondition_...# methods of #SparseMatrix#. + * Base class for other preconditioners. + * Here, only some common features Jacobi, SOR and SSOR preconditioners + * are implemented. For preconditioning, refer to derived classes. * - * \subsection{Use} - * You will usually not want to create a named object of this type, - * although possible. The most common use is like this: - * \begin{verbatim} - * SolverGMRES, - * Vector > gmres(control,memory,500); - * - * gmres.solve (matrix, solution, right_hand_side, - * PreconditionRelaxation,Vector > - * (matrix,&SparseMatrix::template precondition_Jacobi, - * 0.5)); - * \end{verbatim} - * This creates an unnamed object to be passed as the fourth parameter to - * the solver function of the #SolverGMRES# class. It assumes that the - * #SparseMatrix# class has a function #precondition_Jacobi# taking two - * vectors (source and destination) and a relaxation value as parameters. (Unlike - * for the #PreconditionUseMatrix# class, this time it should work, with - * relaxation parameter $0.5$.) - * - * Note that due to the default template parameters, the above example - * could be written shorter as follows: - * \begin{verbatim} - * ... - * gmres.solve (matrix, solution, right_hand_side, - * PreconditionRelaxation<> - * (matrix, - * &SparseMatrix::template precondition_Jacobi, - * 0.5)); - * \end{verbatim} - * - * @author Guido Kanschat, Wolfgang Bangerth, 1999 + * @author Guido Kanschat, 2000 */ -template, class VECTOR = Vector > +template > class PreconditionRelaxation { public: /** - * Type of the preconditioning - * function of the matrix. + * 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. */ - typedef void ( MATRIX::* function_ptr)(VECTOR&, const VECTOR&, - typename MATRIX::value_type) const; - + void initialize (const MATRIX& A, const double omega = 1.); + protected: /** - * Constructor. - * This constructor stores a - * reference to the matrix object - * for later use and selects a - * preconditioning method, which - * must be a member function of - * that matrix. + * Pointer to the matrix object. */ - PreconditionRelaxation(const MATRIX &M, - const function_ptr method, - const double omega = 1.); - + SmartPointer A; /** - * Execute preconditioning. Calls the - * function passed to the constructor - * of this object with the two - * arguments given here, and the - * relaxation parameter passed to the - * constructor. + * Relaxation parameter. */ - void operator() (VECTOR&, const VECTOR&) const; + double omega; +}; - private: +/** + * Jacobi preconditioner using matrix built-in function. The MATRIX + * class used is required to have a function + * @p{precondition_Jacobi(VECTOR&, const VECTOR&, double} + * + * @author Guido Kanschat, 2000 + */ +template > +class PreconditionJacobi : public PreconditionRelaxation +{ + public: /** - * Pointer to the matrix in use. + * Apply preconditioner. */ - const MATRIX& matrix; - + template + void vmult (VECTOR&, const VECTOR&) const; /** - * Pointer to the preconditioning - * function. + * Apply transpose + * preconditioner. Since this is + * a symmetric preconditioner, + * this function is the same as + * @ref{vmult}. */ - const function_ptr precondition; - + template + void Tvmult (VECTOR&, const VECTOR&) const; +}; + + +/** + * SOR preconditioner using matrix built-in function. The MATRIX + * class used is required to have functions + * @p{precondition_SOR(VECTOR&, const VECTOR&, double)} and + * @p{precondition_TSOR(VECTOR&, const VECTOR&, double)}. + * + * @author Guido Kanschat, 2000 + */ +template > +class PreconditionSOR : public PreconditionRelaxation +{ + public: /** - * Relaxation parameter. + * Apply preconditioner. */ - double omega; + template + void vmult (VECTOR&, const VECTOR&) const; + /** + * Apply transpose + * preconditioner. + */ + template + void Tvmult (VECTOR&, const VECTOR&) const; +}; + + +/** + * SSOR preconditioner using matrix built-in function. The MATRIX + * class used is required to have a function + * @p{precondition_SSOR(VECTOR&, const VECTOR&, double} + * + * @author Guido Kanschat, 2000 + */ +template > +class PreconditionSSOR : public PreconditionRelaxation +{ + public: + /** + * Apply preconditioner. + */ + template + void vmult (VECTOR&, const VECTOR&) const; + /** + * Apply transpose + * preconditioner. Since this is + * a symmetric preconditioner, + * this function is the same as + * @ref{vmult}. + */ + template + void Tvmult (VECTOR&, const VECTOR&) const; }; + /** * Preconditioner using an iterative solver. This preconditioner uses * a fully initialized LAC iterative solver for the approximate @@ -307,7 +279,7 @@ class PreconditionLACSolver * Execute preconditioning. */ template - void operator() (VECTOR&, const VECTOR&) const; + void vmult (VECTOR&, const VECTOR&) const; private: /** @@ -351,10 +323,17 @@ class PreconditionedMatrix VectorMemory& mem); /** - * Preconditioned matrix-vector-product. + * Preconditioned + * matrix-vector-product. */ void vmult (VECTOR& dst, const VECTOR& src) const; + /** + * Transposed preconditioned + * matrix-vector-product. + */ + void Tvmult (VECTOR& dst, const VECTOR& src) const; + /** * Residual $b-PAx$. */ @@ -380,13 +359,6 @@ class PreconditionedMatrix /* ---------------------------------- Inline functions ------------------- */ -template -inline void -PreconditionIdentity::operator () (VECTOR& dst, const VECTOR& src) const -{ - dst = src; -} - template inline void PreconditionIdentity::vmult (VECTOR& dst, const VECTOR& src) const @@ -405,56 +377,95 @@ PreconditionIdentity::Tvmult (VECTOR& dst, const VECTOR& src) const template inline void -PreconditionJacobi::initialize (const MATRIX& rA, double o) +PreconditionRelaxation::initialize (const MATRIX& rA, double o) { A = &rA; omega = o; } +//----------------------------------------------------------------------// + template template inline void PreconditionJacobi::vmult (VECTOR& dst, const VECTOR& src) const { -//TODO: Assert object not initialized + Assert (A!=0, ExcNotInitialized()); + A->precondition_Jacobi (dst, src, omega); +} + + +template +template +inline void +PreconditionJacobi::Tvmult (VECTOR& dst, const VECTOR& src) const +{ + Assert (A!=0, ExcNotInitialized()); A->precondition_Jacobi (dst, src, omega); } //----------------------------------------------------------------------// +template +template +inline void +PreconditionSOR::vmult (VECTOR& dst, const VECTOR& src) const +{ + Assert (A!=0, ExcNotInitialized()); + A->precondition_SOR (dst, src, omega); +} -template -PreconditionUseMatrix::PreconditionUseMatrix(const MATRIX& M, - function_ptr method) - : - matrix(M), precondition(method) -{} +template +template +inline void +PreconditionSOR::Tvmult (VECTOR& dst, const VECTOR& src) const +{ + Assert (A!=0, ExcNotInitialized()); + A->precondition_TSOR (dst, src, omega); +} -template -void -PreconditionUseMatrix::operator() (VECTOR& dst, - const VECTOR& src) const + +//----------------------------------------------------------------------// + +template +template +inline void +PreconditionSSOR::vmult (VECTOR& dst, const VECTOR& src) const { - (matrix.*precondition)(dst, src); + Assert (A!=0, ExcNotInitialized()); + A->precondition_SSOR (dst, src, omega); +} + + +template +template +inline void +PreconditionSSOR::Tvmult (VECTOR& dst, const VECTOR& src) const +{ + Assert (A!=0, ExcNotInitialized()); + A->precondition_SSOR (dst, src, omega); } + +//----------------------------------------------------------------------// + + template -PreconditionRelaxation::PreconditionRelaxation(const MATRIX& M, - function_ptr method, - double omega) +PreconditionUseMatrix::PreconditionUseMatrix(const MATRIX& M, + function_ptr method) : - matrix(M), precondition(method), omega(omega) + matrix(M), precondition(method) {} template void -PreconditionRelaxation::operator() (VECTOR& dst, - const VECTOR& src) const +PreconditionUseMatrix::vmult (VECTOR& dst, + const VECTOR& src) const { - (matrix.*precondition)(dst, src, omega); + (matrix.*precondition)(dst, src); } ////////////////////////////////////////////////////////////////////// @@ -471,8 +482,8 @@ PreconditionLACSolver template template void -PreconditionLACSolver::operator() (VECTOR& dst, - const VECTOR& src) const +PreconditionLACSolver::vmult (VECTOR& dst, + const VECTOR& src) const { solver.solve(matrix, dst, src, precondition); } @@ -499,10 +510,25 @@ PreconditionedMatrix VECTOR* h = mem.alloc(); h->reinit(src); A.vmult(*h, src); - P(dst, *h); + P.vmult(dst, *h); mem.free(h); } + +template +inline void +PreconditionedMatrix +::Tvmult (VECTOR& dst, + const VECTOR& src) const +{ + VECTOR* h = mem.alloc(); + h->reinit(src); + A.Tvmult(*h, src); + P.Tvmult(dst, *h); + mem.free(h); +} + + template inline double PreconditionedMatrix @@ -513,7 +539,7 @@ PreconditionedMatrix VECTOR* h = mem.alloc(); h->reinit(src); A.vmult(*h, src); - P(dst, *h); + P.vmult(dst, *h); mem.free(h); dst.sadd(-1.,1.,rhs); return dst.l2_norm (); diff --git a/deal.II/lac/include/lac/precondition_block.h b/deal.II/lac/include/lac/precondition_block.h index bc45f9e907..af9d1be944 100644 --- a/deal.II/lac/include/lac/precondition_block.h +++ b/deal.II/lac/include/lac/precondition_block.h @@ -41,7 +41,7 @@ template class Vector; * For all matrices that are empty above and below the diagonal * blocks (i.e. for all block diagonal matrices) the #BlockJacobi# preconditioner * is a direct solver. For all matrices that are empty only above the diagonal blocks - * (e.g. the matrices one gets by the DG method with downstream numbering) the + * (e.g. the matrices one gets by the DG method with downstream numbering) * #BlockSOR# is a direct solver. * * This first implementation of the #PreconditionBlock# assumes the @@ -84,6 +84,18 @@ class PreconditionBlock: public Subscriptor */ virtual ~PreconditionBlock(); + /** + * Initialize matrix and block + * size. We store the matrix and + * the block size in the + * preconditioner object. In a + * second step, the inverses of + * the diagonal blocks may be + * computed. + */ + void initialize (const SparseMatrix& A, + const unsigned int block_size); + /** * Deletes the inverse diagonal block * matrices if existent, sets the @@ -95,27 +107,20 @@ class PreconditionBlock: public Subscriptor virtual void clear(); /** - * Takes the matrix that should be used - * for the preconditioning. A reference - * to it is stored within this class, - * but ownership of the matrix remains - * with the caller of this function. + * Use only the inverse of the + * first diagonal block to save + * memory and computation time. + * + * Possible applications: + * computing on a cartesian grid, + * all diagonal blocks are the + * same or all diagonal blocks + * are at least similar and + * inversion of one of them still + * yields a preconditioner. */ - void use_matrix(const SparseMatrix &M); + void set_same_diagonal (); - /** - * Set the right block size before calling - * #invert_diagblocks# or calling the - * #operator ()# function - * of #PreconditionBlockSOR# - * or #PreconditionBlockJacobi#. - * If #block_size==1# BlockSOR or - * BlockJacobi are equal to the - * standard SOR or Jacobi - * preconditioner. - */ - void set_block_size (const unsigned int bsize); - /** * Stores the inverse of * the diagonal blocks @@ -147,6 +152,12 @@ class PreconditionBlock: public Subscriptor */ unsigned int block_size () const; + /** + * Determine, whether inverses + * have been computed. + */ + bool inverses_ready () const; + /** * Exception */ @@ -170,29 +181,37 @@ class PreconditionBlock: public Subscriptor */ DeclException0 (ExcInverseMatricesAlreadyExist); - /** - * Exception - */ - DeclException0 (ExcBlockSizeNotSet); - /** * Exception */ DeclException0 (ExcMatrixNotSquare); - - /** - * Exception - */ - DeclException0 (ExcNoMatrixGivenToUse); protected: + /** + * Access to the inverse diagonal blocks. + * + */ + const FullMatrix& inverse (unsigned int i) const; + /** * Size of the blocks. Each diagonal * block is assumed to be of the * same size. */ unsigned int blocksize; - + + /** + * Pointer to the matrix. Make sure that + * the matrix exists as long as this class + * needs it, i.e. until calling #invert_diagblocks#, + * or (if the inverse matrices should not be + * stored) until the last call of the + * preconditoining #vmult# function of the + * derived classes. + */ + SmartPointer > A; + + private: /** * Storage of the inverse matrices of * the diagonal blocks matrices as @@ -201,18 +220,13 @@ class PreconditionBlock: public Subscriptor * using #inverse_type=float# saves memory * in comparison with #inverse_type=double#. */ - vector > inverse; + vector > _inverse; /** - * Pointer to the matrix. Make sure that - * the matrix exists as long as this class - * needs it, i.e. until calling #invert_diagblocks#, - * or (if the inverse matrices should not be - * stored) until the last call of the - * preconditoining #operator()# function of the - * derived classes. + * Flag for diagonal compression. + * @see set_same_diagonal() */ - SmartPointer > A; + bool same_diagonal; }; @@ -226,11 +240,6 @@ class PreconditionBlockJacobi : public PreconditionBlock public: /** * Execute block Jacobi preconditioning. - * Make sure that the right - * block size - * of the matrix is set by - * #set_block_size# - * before calling this function. * * This function will automatically use the * inverse matrices if they exist, if not @@ -239,7 +248,13 @@ class PreconditionBlockJacobi : public PreconditionBlock * matrices in each preconditioning step. */ template - void operator() (Vector&, const Vector&) const; + void vmult (Vector&, const Vector&) const; + + /** + * Same as #vmult#, since Jacobi is symmetric. + */ + template + void Tvmult (Vector&, const Vector&) const; }; @@ -273,11 +288,6 @@ class PreconditionBlockSOR : public PreconditionBlock /** * Execute block SOR preconditioning. - * Make sure that the right - * block size - * of the matrix is set by - * #set_block_size# - * before calling this function. * * This function will automatically use the * inverse matrices if they exist, if not @@ -290,7 +300,7 @@ class PreconditionBlockSOR : public PreconditionBlock * BlockSOR is a direct solver. */ template - void operator() (Vector&, const Vector&) const; + void vmult (Vector&, const Vector&) const; private: /** diff --git a/deal.II/lac/include/lac/precondition_block.templates.h b/deal.II/lac/include/lac/precondition_block.templates.h index ccf1790635..04ac005662 100644 --- a/deal.II/lac/include/lac/precondition_block.templates.h +++ b/deal.II/lac/include/lac/precondition_block.templates.h @@ -14,6 +14,7 @@ #define __deal2__precondition_block_templates_h +#include #include #include #include @@ -23,103 +24,155 @@ template PreconditionBlock::PreconditionBlock (): blocksize(0), - A(0) {}; + A(0), + same_diagonal(false) +{}; template PreconditionBlock::~PreconditionBlock () { - if (inverse.size()!=0) - inverse.erase(inverse.begin(), inverse.end()); + if (_inverse.size()!=0) + _inverse.erase(_inverse.begin(), _inverse.end()); } template void PreconditionBlock::clear () { - if (inverse.size()!=0) - inverse.erase(inverse.begin(), inverse.end()); - blocksize=0; + if (_inverse.size()!=0) + _inverse.erase(_inverse.begin(), _inverse.end()); + blocksize = 0; + same_diagonal = false; } template -void PreconditionBlock::use_matrix( - const SparseMatrix &M) +void PreconditionBlock +::initialize (const SparseMatrix &M, unsigned int bsize) { + Assert (M.m() == M.n(), ExcMatrixNotSquare()); A = &M; + Assert (bsize>0, ExcIndexRange(bsize, 1, M.m())); + Assert (A->m()%bsize==0, ExcWrongBlockSize(bsize, A->m())); + blocksize=bsize; } template -void PreconditionBlock::set_block_size(unsigned int bsize) { - blocksize=bsize; +const FullMatrix& +PreconditionBlock::inverse(unsigned int i) const +{ + if (same_diagonal) + return _inverse[0]; + + Assert (i < _inverse.size(), ExcIndexRange(i,0,_inverse.size())); + return _inverse[i]; } template -unsigned int PreconditionBlock::block_size() const { +unsigned int PreconditionBlock::block_size() const +{ return blocksize; } +template +void +PreconditionBlock::set_same_diagonal() +{ + Assert(_inverse.size()==0, ExcInverseMatricesAlreadyExist()); + same_diagonal = true; +} + + +template +bool +PreconditionBlock::inverses_ready() const +{ + return (_inverse.size() != 0); +} + + template void PreconditionBlock::invert_diagblocks() { - Assert (A!=0, ExcNoMatrixGivenToUse()); - const SparseMatrix &M=*A; - Assert (M.m() == M.n(), ExcMatrixNotSquare()); - Assert (inverse.size()==0, ExcInverseMatricesAlreadyExist()); + Assert (A!=0, ExcNotInitialized()); + Assert (blocksize!=0, ExcNotInitialized()); - Assert (blocksize!=0, ExcBlockSizeNotSet()); - Assert (M.m()%blocksize==0, ExcWrongBlockSize(blocksize, M.m())); + const SparseMatrix &M=*A; + Assert (_inverse.size()==0, ExcInverseMatricesAlreadyExist()); const unsigned int n_cells = M.m()/blocksize; - // cell_row, cell_column are the - // numbering of the blocks (cells). - // row_cell, column_cell are the local - // numbering of the unknowns in the - // blocks. - // row, column are the global numbering - // of the unkowns. - - // set the #inverse# array to the right - // size. we could do it like this: - // inverse = vector<>(n_cells,FullMatrix<>()) - // but this would involve copying many - // FullMatrix objects. - // - // the following is a neat trick which - // avoids copying - if (true) - { - vector > tmp(n_cells, - FullMatrix(blocksize)); - inverse.swap (tmp); - }; - FullMatrix M_cell(blocksize); - - for (unsigned int cell=0, row=0; cell(n_cells,FullMatrix<>()) + // but this would involve copying many + // FullMatrix objects. + // + // the following is a neat trick which + // avoids copying + vector > tmp(n_cells, + FullMatrix(blocksize)); + _inverse.swap (tmp); + + M_cell.clear (); + + for (unsigned int cell=0, row=0; cell::invert_diagblocks() template template void PreconditionBlockJacobi -::operator() (Vector &dst, - const Vector &src) const +::vmult (Vector &dst, + const Vector &src) const { // introduce the following typedef // since in the use of exceptions, @@ -150,20 +203,11 @@ void PreconditionBlockJacobi // exceptions that do not take // args... typedef PreconditionBlock BaseClass; - Assert(A!=0, typename BaseClass::ExcNoMatrixGivenToUse()); + Assert(A!=0, ExcNotInitialized()); const SparseMatrix &M=*A; const unsigned int n_cells=M.m()/blocksize; - Assert (M.m() == M.n(), - typename BaseClass::ExcMatrixNotSquare()); - Assert (blocksize!=0, - typename BaseClass::ExcBlockSizeNotSet()); - Assert (M.m()%blocksize==0, - typename BaseClass::ExcWrongBlockSize(blocksize, M.m())); - Assert (inverse.size()==0 || inverse.size()==n_cells, - typename BaseClass::ExcWrongNumberOfInverses(inverse.size(), n_cells)); - Vector b_cell(blocksize), x_cell(blocksize); // cell_row, cell_column are the @@ -175,7 +219,7 @@ void PreconditionBlockJacobi // of the unkowns. unsigned int row, row_cell, begin_diag_block=0; - if (inverse.size()==0) + if (!inverses_ready()) { FullMatrix M_cell(blocksize); for (unsigned int cell=0; cell { b_cell(row_cell)=src(row); } - inverse[cell].vmult(x_cell, b_cell); + inverse(cell).vmult(x_cell, b_cell); // distribute x_cell to dst for (row=cell*blocksize, row_cell=0; row_cell begin_diag_block+=blocksize; } } + + +template +template +void PreconditionBlockJacobi +::Tvmult (Vector &dst, + const Vector &src) const +{ + vmult(dst, src); +} + + /*--------------------- PreconditionBlockSOR -----------------------*/ template PreconditionBlockSOR::PreconditionBlockSOR(const number omega): - omega(omega) {} + omega(omega) +{} template -PreconditionBlockSOR::~PreconditionBlockSOR(){} +PreconditionBlockSOR::~PreconditionBlockSOR() +{} @@ -231,12 +289,12 @@ PreconditionBlockSOR::set_omega(number om) } - +//TODO: implement Tvmult template template -void PreconditionBlockSOR::operator() (Vector &dst, - const Vector &src) const +void PreconditionBlockSOR::vmult (Vector &dst, + const Vector &src) const { // introduce the following typedef // since in the use of exceptions, @@ -256,20 +314,11 @@ void PreconditionBlockSOR::operator() (Vector // args... typedef PreconditionBlock BaseClass; - Assert (A!=0, typename BaseClass::ExcNoMatrixGivenToUse()); + Assert (A!=0, ExcNotInitialized()); const SparseMatrix &M=*A; const unsigned int n_cells=M.m()/blocksize; - Assert (M.m() == M.n(), - typename BaseClass::ExcMatrixNotSquare()); - Assert (blocksize!=0, - typename BaseClass::ExcBlockSizeNotSet()); - Assert (M.m()%blocksize==0, - typename BaseClass::ExcWrongBlockSize(blocksize, M.m())); - Assert (inverse.size()==0 || inverse.size()==n_cells, - typename BaseClass::ExcWrongNumberOfInverses(inverse.size(), n_cells)); - const SparsityPattern &spars = M.get_sparsity_pattern(); const unsigned int *rowstart = spars.get_rowstart_indices(); const unsigned int *columns = spars.get_column_numbers(); @@ -286,7 +335,7 @@ void PreconditionBlockSOR::operator() (Vector unsigned int row, column, row_cell, begin_diag_block=0; number2 b_cell_row; - if (inverse.size()==0) + if (!inverses_ready()) { FullMatrix M_cell(blocksize); for (unsigned int cell=0; cell::operator() (Vector } b_cell(row_cell)=b_cell_row; } - inverse[cell].vmult(x_cell, b_cell); + inverse(cell).vmult(x_cell, b_cell); // distribute x_cell to dst for (row=cell*blocksize, row_cell=0; row_cell /** * Solve primal problem only. */ - template + template typename Solver::ReturnState solve (const MATRIX &A, VECTOR &x, const VECTOR &b, - const Preconditioner& precondition); + const PRECONDITIONER& precondition); protected: /** @@ -257,7 +257,7 @@ SolverBicgstab::iterate(const MATRIX& A, beta = rhobar * alpha / (rho * omega); rho = rhobar; p.sadd(beta, 1., r, -beta*omega, v); - precondition(y,p); + precondition.vmult(y,p); A.vmult(v,y); rhobar = rbar * v; @@ -269,7 +269,7 @@ SolverBicgstab::iterate(const MATRIX& A, return typename Solver::ReturnState(breakdown); s.equ(1., r, -alpha, v); - precondition(z,s); + precondition.vmult(z,s); A.vmult(t,z); rhobar = t*s; omega = rhobar/(t*t); diff --git a/deal.II/lac/include/lac/solver_cg.h b/deal.II/lac/include/lac/solver_cg.h index 1b7f41e5e6..0b2d344ce1 100644 --- a/deal.II/lac/include/lac/solver_cg.h +++ b/deal.II/lac/include/lac/solver_cg.h @@ -70,12 +70,12 @@ class SolverCG : private Solver /** * Solver method. */ - template + template typename Solver::ReturnState solve (const MATRIX &A, VECTOR &x, const VECTOR &b, - const Preconditioner& precondition); + const PRECONDITIONER& precondition); protected: /** @@ -164,12 +164,12 @@ SolverCG::print_vectors(const unsigned int, template -template +template typename Solver::ReturnState SolverCG::solve (const MATRIX &A, VECTOR &x, const VECTOR &b, - const Preconditioner& precondition) + const PRECONDITIONER& precondition) { SolverControl::State conv=SolverControl::iterate; @@ -210,7 +210,7 @@ SolverCG::solve (const MATRIX &A, }; g.scale(-1.); - precondition(h,g); + precondition.vmult(h,g); d.equ(-1.,h); @@ -234,7 +234,7 @@ SolverCG::solve (const MATRIX &A, if (conv) break; - precondition(h,g); + precondition.vmult(h,g); beta = gh; gh = g*h; diff --git a/deal.II/lac/include/lac/solver_gmres.h b/deal.II/lac/include/lac/solver_gmres.h index 2dbf439fe8..d264fdb8ad 100644 --- a/deal.II/lac/include/lac/solver_gmres.h +++ b/deal.II/lac/include/lac/solver_gmres.h @@ -271,7 +271,7 @@ unsigned int dim = 0; if (left_precondition) { A.residual(p,x,b); - precondition(v,p); + precondition.vmult(v,p); } else { A.residual(v,x,b); }; @@ -311,9 +311,9 @@ unsigned int dim = 0; if (left_precondition) { A.vmult(p, *tmp_vectors[inner_iteration]); - precondition(vv,p); + precondition.vmult(vv,p); } else { - precondition(p,*tmp_vectors[inner_iteration]); + precondition.vmult(p,*tmp_vectors[inner_iteration]); A.vmult(vv,p); }; @@ -376,7 +376,7 @@ unsigned int dim = 0; p = 0.; for (unsigned int i=0; i /** * Solver method. */ - template + template typename Solver::ReturnState solve (const MATRIX &A, VECTOR &x, const VECTOR &b, - const Preconditioner& precondition); + const PRECONDITIONER& precondition); /** * Exception @@ -168,12 +168,12 @@ SolverMinRes::print_vectors(const unsigned int, template -template +template typename Solver::ReturnState SolverMinRes::solve (const MATRIX &A, VECTOR &x, const VECTOR &b, - const Preconditioner& precondition) + const PRECONDITIONER& precondition) { SolverControl::State conv=SolverControl::iterate; @@ -234,7 +234,7 @@ SolverMinRes::solve (const MATRIX &A, // positiv definite and symmetric // M v = u[1] - precondition (v,u[1]); + precondition.vmult (v,u[1]); delta[1] = v * u[1]; // Preconditioner positive @@ -270,7 +270,7 @@ SolverMinRes::solve (const MATRIX &A, // precondition: solve M v = u[2] // Preconditioner has to be positiv // definite and symmetric. - precondition(v,u[2]); + precondition.vmult(v,u[2]); delta[2] = v * u[2]; diff --git a/deal.II/lac/include/lac/solver_qmrs.h b/deal.II/lac/include/lac/solver_qmrs.h index 20ce79596f..31946803a8 100644 --- a/deal.II/lac/include/lac/solver_qmrs.h +++ b/deal.II/lac/include/lac/solver_qmrs.h @@ -102,12 +102,12 @@ class SolverQMRS : private Solver /** * Solver method. */ - template + template typename Solver::ReturnState solve (const MATRIX &A, VECTOR &x, const VECTOR &b, - const Preconditioner& precondition); + const PRECONDITIONER& precondition); /** * Interface for derived class. @@ -168,9 +168,9 @@ class SolverQMRS : private Solver /** * The iteration loop itself. */ - template + template typename Solver::ReturnState - iterate(const MATRIX& A, const Preconditioner& precondition); + iterate(const MATRIX& A, const PRECONDITIONER& precondition); /** * The current iteration step. */ @@ -210,12 +210,12 @@ SolverQMRS::print_vectors(const unsigned int, template -template +template typename Solver::ReturnState SolverQMRS::solve (const MATRIX &A, VECTOR &x, const VECTOR &b, - const Preconditioner& precondition) + const PRECONDITIONER& precondition) { deallog.push("QMRS"); @@ -266,10 +266,10 @@ SolverQMRS::solve (const MATRIX &A, template -template +template typename Solver::ReturnState SolverQMRS::iterate(const MATRIX& A, - const Preconditioner& precondition) + const PRECONDITIONER& precondition) { /* Remark: the matrix A in the article is the preconditioned matrix. * Therefore, we have to precondition x before we compute the first residual. @@ -296,7 +296,7 @@ SolverQMRS::iterate(const MATRIX& A, d.reinit(x); // Apply right preconditioning to x - precondition(q,x); + precondition.vmult(q,x); // Preconditioned residual res = A.residual(v, q, b); @@ -305,7 +305,7 @@ SolverQMRS::iterate(const MATRIX& A, p = v; - precondition(q,p); + precondition.vmult(q,p); tau = v.norm_sqr(); //deallog << "tau:" << tau << endl; @@ -359,12 +359,12 @@ while (state == SolverControl::iterate) return breakdown; // Step 7 rho_old = rho; - precondition(q,v); + precondition.vmult(q,v); rho = q*v; beta = rho/rho_old; p.sadd(beta,v); - precondition(q,p); + precondition.vmult(q,p); } return exceeded; } diff --git a/deal.II/lac/include/lac/solver_richardson.h b/deal.II/lac/include/lac/solver_richardson.h index eb5a6332e9..fc6ad64a47 100644 --- a/deal.II/lac/include/lac/solver_richardson.h +++ b/deal.II/lac/include/lac/solver_richardson.h @@ -73,11 +73,11 @@ class SolverRichardson : private Solver /** * Solve $Ax=b$ for $x$. */ - template + template typename Solver::ReturnState solve (const MATRIX &A, VECTOR &x, const VECTOR &b, - const Preconditioner& precondition); + const PRECONDITIONER& precondition); /** * Set the damping-coefficient. @@ -146,12 +146,12 @@ SolverRichardson::SolverRichardson(SolverControl &cn, template -template +template typename Solver::ReturnState SolverRichardson::solve (const MATRIX &A, VECTOR &x, const VECTOR &b, - const Preconditioner& precondition) + const PRECONDITIONER& precondition) { SolverControl::State conv=SolverControl::iterate; @@ -170,7 +170,7 @@ SolverRichardson::solve (const MATRIX &A, if (conv != SolverControl::iterate) break; - precondition(d,r); + precondition.vmult(d,r); x.add(additional_data.omega,d); print_vectors(iter,x,r,d); } diff --git a/deal.II/lac/include/lac/sparse_ilu.h b/deal.II/lac/include/lac/sparse_ilu.h index e092baa40a..2dae4fe823 100644 --- a/deal.II/lac/include/lac/sparse_ilu.h +++ b/deal.II/lac/include/lac/sparse_ilu.h @@ -192,8 +192,8 @@ class SparseILU : protected SparseMatrix * Same as #apply_decomposition#, format for LAC. */ template - void operator() (Vector &dst, - const Vector &src) const; + void vmult (Vector &dst, + const Vector &src) const; /** * Exception @@ -218,8 +218,8 @@ class SparseILU : protected SparseMatrix template template void -SparseILU::operator() (Vector &dst, - const Vector &src) const +SparseILU::vmult (Vector &dst, + const Vector &src) const { apply_decomposition(dst, src); } diff --git a/deal.II/lac/include/lac/sparse_vanka.h b/deal.II/lac/include/lac/sparse_vanka.h index cb318d9b94..bffe5898ac 100644 --- a/deal.II/lac/include/lac/sparse_vanka.h +++ b/deal.II/lac/include/lac/sparse_vanka.h @@ -203,8 +203,8 @@ class SparseVanka * update vector in #dst#. */ template - void operator() (Vector &dst, - const Vector &src) const; + void vmult (Vector &dst, + const Vector &src) const; /** * Exception @@ -256,7 +256,7 @@ class SparseVanka * function with a * #vector(n_dofs,true)#. * - * The #operator()# of this class + * The #vmult# of this class * of course calls this function * with a null pointer */ @@ -516,7 +516,7 @@ class SparseBlockVanka : public SparseVanka * Apply the preconditioner. */ template - void operator() (Vector &dst, + void vmult (Vector &dst, const Vector &src) const; private: diff --git a/deal.II/lac/include/lac/sparse_vanka.templates.h b/deal.II/lac/include/lac/sparse_vanka.templates.h index bb77852398..b3d0df796d 100644 --- a/deal.II/lac/include/lac/sparse_vanka.templates.h +++ b/deal.II/lac/include/lac/sparse_vanka.templates.h @@ -186,13 +186,13 @@ SparseVanka::compute_inverse (const unsigned int row, inverses[row]->gauss_jordan(); }; - +//TODO: implement Tvmult template template void -SparseVanka::operator ()(Vector &dst, - const Vector &src) const +SparseVanka::vmult (Vector &dst, + const Vector &src) const { // first set output vector to zero dst.clear (); @@ -543,11 +543,12 @@ SparseBlockVanka::compute_dof_masks (const SparseMatrix &M, }; }; +//TODO: implement Tvmult template template -void SparseBlockVanka::operator() (Vector &dst, - const Vector &src) const +void SparseBlockVanka::vmult (Vector &dst, + const Vector &src) const { dst.clear (); diff --git a/deal.II/lac/source/precondition_block.cc b/deal.II/lac/source/precondition_block.cc index d94b57ac28..426c93ce76 100644 --- a/deal.II/lac/source/precondition_block.cc +++ b/deal.II/lac/source/precondition_block.cc @@ -35,25 +35,36 @@ template class PreconditionBlock; // explicit instantiations for "float" PreconditionBlock template class PreconditionBlockJacobi; -template void PreconditionBlockJacobi::operator() ( - Vector &, const Vector &) const; -template void PreconditionBlockJacobi::operator() ( - Vector &, const Vector &) const; - +template void PreconditionBlockJacobi::vmult +(Vector &, const Vector &) const; +template void PreconditionBlockJacobi::vmult +(Vector &, const Vector &) const; +template void PreconditionBlockJacobi::Tvmult +(Vector &, const Vector &) const; +template void PreconditionBlockJacobi::Tvmult +(Vector &, const Vector &) const; template class PreconditionBlockJacobi; -template void PreconditionBlockJacobi::operator() ( - Vector &, const Vector &) const; -template void PreconditionBlockJacobi::operator() ( - Vector &, const Vector &) const; +template void PreconditionBlockJacobi::vmult +(Vector &, const Vector &) const; +template void PreconditionBlockJacobi::vmult +(Vector &, const Vector &) const; +template void PreconditionBlockJacobi::Tvmult +(Vector &, const Vector &) const; +template void PreconditionBlockJacobi::Tvmult +(Vector &, const Vector &) const; template class PreconditionBlockJacobi; -template void PreconditionBlockJacobi::operator() ( - Vector &, const Vector &) const; -template void PreconditionBlockJacobi::operator() ( - Vector &, const Vector &) const; +template void PreconditionBlockJacobi::vmult +(Vector &, const Vector &) const; +template void PreconditionBlockJacobi::vmult +(Vector &, const Vector &) const; +template void PreconditionBlockJacobi::Tvmult +(Vector &, const Vector &) const; +template void PreconditionBlockJacobi::Tvmult +(Vector &, const Vector &) const; /*--------------------- PreconditionBlockGaussSeidel -----------------------*/ @@ -61,9 +72,9 @@ template void PreconditionBlockJacobi::operator() ( // explicit instantiations for "float" PreconditionBlock template class PreconditionBlockSOR; -template void PreconditionBlockSOR::operator() ( +template void PreconditionBlockSOR::vmult ( Vector &, const Vector &) const; -template void PreconditionBlockSOR::operator() ( +template void PreconditionBlockSOR::vmult ( Vector &, const Vector &) const; @@ -76,16 +87,16 @@ template void PreconditionBlockSOR::operator() ( template class PreconditionBlockSOR; -template void PreconditionBlockSOR::operator() ( +template void PreconditionBlockSOR::vmult ( Vector &, const Vector &) const; -template void PreconditionBlockSOR::operator() ( +template void PreconditionBlockSOR::vmult ( Vector &, const Vector &) const; template class PreconditionBlockSOR; -template void PreconditionBlockSOR::operator() ( +template void PreconditionBlockSOR::vmult ( Vector &, const Vector &) const; -template void PreconditionBlockSOR::operator() ( +template void PreconditionBlockSOR::vmult ( Vector &, const Vector &) const; diff --git a/deal.II/lac/source/sparse_vanka.cc b/deal.II/lac/source/sparse_vanka.cc index b23b455b4d..9828eea8dd 100644 --- a/deal.II/lac/source/sparse_vanka.cc +++ b/deal.II/lac/source/sparse_vanka.cc @@ -19,16 +19,16 @@ template class SparseVanka; template class SparseVanka; -template void SparseVanka::operator () (Vector &dst, - const Vector &src) const; -template void SparseVanka::operator () (Vector &dst, - const Vector &src) const; +template void SparseVanka::vmult (Vector &dst, + const Vector &src) const; +template void SparseVanka::vmult (Vector &dst, + const Vector &src) const; template class SparseBlockVanka; template class SparseBlockVanka; -template void SparseBlockVanka::operator () (Vector &dst, - const Vector &src) const; -template void SparseBlockVanka::operator () (Vector &dst, - const Vector &src) const; +template void SparseBlockVanka::vmult (Vector &dst, + const Vector &src) const; +template void SparseBlockVanka::vmult (Vector &dst, + const Vector &src) const; diff --git a/tests/big-tests/multigrid/multigrid.cc b/tests/big-tests/multigrid/multigrid.cc index 891b8cbe58..d28e210a3e 100644 --- a/tests/big-tests/multigrid/multigrid.cc +++ b/tests/big-tests/multigrid/multigrid.cc @@ -154,8 +154,8 @@ MGSmootherLAC::smooth (const unsigned int level, SolverControl control(2,1.e-300,false,false); PrimitiveVectorMemory<> mem; SolverRichardson<> rich(control, mem); - PreconditionRelaxation<> - prec((*matrices)[level], &SparseMatrix ::template precondition_SSOR, 1.); + PreconditionSSOR<> prec; + prec.initialize((*matrices)[level], 1.); rich.solve((*matrices)[level], u, rhs, prec); } @@ -377,10 +377,8 @@ void LaplaceProblem::solve () PrimitiveVectorMemory<> vector_memory; SolverCG<> cg (solver_control, vector_memory); - PreconditionRelaxation<> - preconditioner(global_system_matrix, - &SparseMatrix::template precondition_SSOR, - 1.2); + PreconditionSSOR<> preconditioner; + preconditioner.initialize(global_system_matrix, 1.2); solution.clear (); cg.solve (global_system_matrix, solution, system_rhs, diff --git a/tests/deal.II/block_matrices.cc b/tests/deal.II/block_matrices.cc index 7eb3c1c726..d3e75ea572 100644 --- a/tests/deal.II/block_matrices.cc +++ b/tests/deal.II/block_matrices.cc @@ -266,8 +266,8 @@ void LaplaceProblem::solve () PrimitiveVectorMemory vector_memory; SolverCG cg (solver_control, vector_memory); - PreconditionRelaxation preconditioner - (system_matrix, &Matrix::precondition_Jacobi, 0.8); + PreconditionJacobi preconditioner; + preconditioner.initialize (system_matrix, 0.8); cg.solve (system_matrix, solution, system_rhs, preconditioner); diff --git a/tests/deal.II/mg.cc b/tests/deal.II/mg.cc index 84cfb6f38e..52241020bf 100644 --- a/tests/deal.II/mg.cc +++ b/tests/deal.II/mg.cc @@ -180,8 +180,8 @@ MGSmootherLAC::smooth (const unsigned int level, SolverControl control(2,1.e-300,false,false); PrimitiveVectorMemory<> mem; SolverRichardson<> rich(control, mem); - PreconditionRelaxation<> - prec((*matrices)[level], &SparseMatrix ::template precondition_SSOR, 1.); + PreconditionSSOR<> prec; + prec.initialize((*matrices)[level], 1.); rich.solve((*matrices)[level], u, rhs, prec); }