From: Jean-Paul Pelteret Date: Thu, 24 Sep 2015 14:30:05 +0000 (+0200) Subject: Ammended native deal.II preconditioners and solvers for use in LinearOperators. X-Git-Tag: v8.4.0-rc2~296^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=67c8e925d4107969c786033b5db096ed92c46de3;p=dealii.git Ammended native deal.II preconditioners and solvers for use in LinearOperators. Improved: Native deal.II preconditioners can now be directly wrapped as LinearOperators. This allows them to be used in the construction of more complex linear solvers, such as where the approximate inverse to a matrix is to be defined as a preconditioner to said matrix. New: Added a LinearOperator test for deal.II SparseMatrix solvers and preconditioners. Added a test to check that the various combinations of deal.II solvers and preconditioners (for the SparseMatrix type) work with the linear_operator and inverse_operator functions. By 'work' one means that they both compile and also produce a sensible answer. Note that the preconditioners are expected to act as per usual (i.e. in the preconditioner slot in inverse_operator) as well as a standard matrix type in linear_operator. This is necessary as they will be used in this context later when a schur_matrix_operator and approximate_schur_matrix_operator is developed. Fixed: Reinitialised (zeroed) range vectors in the construction of an inverse_operator before the solve calls. This corrects a bug wherein the starting point of the solver was garbage leading to potential divergence in the linear solver. --- diff --git a/doc/news/changes.h b/doc/news/changes.h index 2449e7eebb..3edcdc60fb 100644 --- a/doc/news/changes.h +++ b/doc/news/changes.h @@ -170,6 +170,17 @@ inconvenience this causes. (Timo Heister, 2015/09/30) +
  • Improved: The interface to all deal.II type solvers and preconditioners have + been updated such that they function as expected with the LinearOperator class and + its associated functions (i.e. linear_operator, transpose_operator and + inverse_operator). These preconditioners can now be wrapped as a LinearOperator, + facilitating the construction of approximate matrix inverses such as in the + development of a block matrix preconditioner. + An example of this functionality can be found in tests/lac/linear_operator_08.cc . +
    + (Jean-Paul Pelteret, 2015/09/24 - 2015/10/19) +
  • +
  • New: MGTransferPrebuilt can now be used with parallel::distributed::Vector and TrilinosWrappers::SparseMatrix as a transfer matrix.
    @@ -372,6 +383,13 @@ inconvenience this causes. (Alberto Sartori, 2015/10/22)
  • +
  • Fixed: The range vectors in the construction of an inverse_operator + is now reinitialised before solve calls. This ensures a consistent starting + point for the solver. +
    + (Jean-Paul Pelteret, 2015/10/19) +
  • +
  • New: Ghost cells for the multigrid levels in parallel::distributed::Triangulation are now correctly created also for periodic boundary conditions. diff --git a/include/deal.II/lac/iterative_inverse.h b/include/deal.II/lac/iterative_inverse.h index dbbf17ccfc..84375351cf 100644 --- a/include/deal.II/lac/iterative_inverse.h +++ b/include/deal.II/lac/iterative_inverse.h @@ -72,13 +72,19 @@ DEAL_II_NAMESPACE_OPEN * for further details. * * @ingroup Matrix2 - * @author Guido Kanschat - * @date 2010 + * @author Guido Kanschat; extension for full compatibility with + * LinearOperator class: Jean-Paul Pelteret + * @date 2010, 2015 */ template class IterativeInverse : public Subscriptor { public: + /** + * Declare type for container size. + */ + typedef typename PointerMatrixBase::size_type size_type; + /** * Initialization function. Provide a matrix and preconditioner for the * solve in vmult(). @@ -96,6 +102,11 @@ public: */ void vmult (VECTOR &dst, const VECTOR &src) const; + /** + * Same as before, but uses the transpose of the matrix. + */ + void Tvmult (VECTOR &dst, const VECTOR &src) const; + /** * Solve for right hand side src, but allow for the fact that the * vectors given to this function have different type from the vectors used @@ -104,6 +115,24 @@ public: template void vmult (VECTOR2 &dst, const VECTOR2 &src) const; + /** + * Same as before, but uses the transpose of the matrix. + */ + template + void Tvmult (VECTOR2 &dst, const VECTOR2 &src) const; + + /** + * Return the dimension of the codomain (or range) space. To remember: the + * matrix is of dimension $m \times n$. + */ + size_type m () const; + + /** + * Return the dimension of the domain space. To remember: the matrix is of + * dimension $m \times n$. + */ + size_type n () const; + /** * The solver, which allows selection of the actual solver as well as * adjustment of parameters. @@ -157,6 +186,14 @@ IterativeInverse::vmult (VECTOR &dst, const VECTOR &src) const } +template +inline void +IterativeInverse::Tvmult (VECTOR &dst, const VECTOR &src) const +{ + AssertThrow(false, ExcNotImplemented()); +} + + template template inline void @@ -174,6 +211,32 @@ IterativeInverse::vmult (VECTOR2 &dst, const VECTOR2 &src) const } +template +template +inline void +IterativeInverse::Tvmult (VECTOR2 &dst, const VECTOR2 &src) const +{ + AssertThrow(false, ExcNotImplemented()); +} + + +template +inline typename IterativeInverse::size_type +IterativeInverse::m () const +{ + Assert (matrix.get()!=0, ExcNotInitialized()); + return matrix->m(); +} + +template +inline typename IterativeInverse::size_type +IterativeInverse::n () const +{ + Assert (matrix.get()!=0, ExcNotInitialized()); + return matrix->n(); +} + + DEAL_II_NAMESPACE_CLOSE diff --git a/include/deal.II/lac/linear_operator.h b/include/deal.II/lac/linear_operator.h index 2b4808d403..eae3309d82 100644 --- a/include/deal.II/lac/linear_operator.h +++ b/include/deal.II/lac/linear_operator.h @@ -641,6 +641,7 @@ inverse_operator(const LinearOperator vector_memory; Vector *v2 = vector_memory.alloc(); - op.reinit_range_vector(*v2, /*bool fast =*/ true); + op.reinit_range_vector(*v2, /*bool fast =*/ false); solver.solve(op, *v2, u, preconditioner); v += *v2; vector_memory.free(v2); @@ -659,6 +660,7 @@ inverse_operator(const LinearOperator vector_memory; Vector *v2 = vector_memory.alloc(); - op.reinit_range_vector(*v2, /*bool fast =*/ true); + op.reinit_range_vector(*v2, /*bool fast =*/ false); solver.solve(transpose_operator(op), *v2, u, preconditioner); v += *v2; vector_memory.free(v2); diff --git a/include/deal.II/lac/pointer_matrix.h b/include/deal.II/lac/pointer_matrix.h index 708f595e1c..4b27ee6bcb 100644 --- a/include/deal.II/lac/pointer_matrix.h +++ b/include/deal.II/lac/pointer_matrix.h @@ -51,6 +51,11 @@ template class PointerMatrixBase : public Subscriptor { public: + /** + * Declare type for container size. + */ + typedef types::global_dof_index size_type; + /** * Value type of this matrix. since the matrix itself is unknown, we take * the value type of the vector. Therefore, matrix entries must be @@ -96,6 +101,18 @@ public: */ virtual void Tvmult_add (VECTOR &dst, const VECTOR &src) const = 0; + + /** + * Return the dimension of the codomain (or range) space. To remember: the + * matrix is of dimension $m \times n$. + */ + virtual size_type m () const = 0; + + /** + * Return the dimension of the domain space. To remember: the matrix is of + * dimension $m \times n$. + */ + virtual size_type n () const = 0; }; @@ -107,12 +124,18 @@ public: * only has a vector as template argument. Therefore, this interface provides * an abstract base class for matrices. * - * @author Guido Kanschat 2000, 2001, 2002 + * @author Guido Kanschat 2000, 2001, 2002; extension for full compatibility + * with LinearOperator class: Jean-Paul Pelteret, 2015 */ template class PointerMatrix : public PointerMatrixBase { public: + /** + * Declare type for container size. + */ + typedef typename PointerMatrixBase::size_type size_type; + /** * Constructor. The pointer in the argument is stored in this class. As * usual, the lifetime of *M must be longer than the one of the @@ -148,6 +171,12 @@ public: PointerMatrix(const MATRIX *M, const char *name); + /** + * Destructor + */ + virtual ~PointerMatrix () + {} + // Use doc from base class virtual void clear(); @@ -187,11 +216,23 @@ public: virtual void Tvmult_add (VECTOR &dst, const VECTOR &src) const; + /** + * Return the dimension of the codomain (or range) space. To remember: the + * matrix is of dimension $m \times n$. + */ + virtual size_type m () const; + + /** + * Return the dimension of the domain space. To remember: the matrix is of + * dimension $m \times n$. + */ + virtual size_type n () const; + private: /** * The pointer to the actual matrix. */ - SmartPointer > m; + SmartPointer > mtrx; }; @@ -207,12 +248,18 @@ private: * and by the fact that it implements the functions vmult_add() and * Tvmult_add() only using vmult() and Tvmult() of the MATRIX. * - * @author Guido Kanschat 2006 + * @author Guido Kanschat 2006; extension for full compatibility with + * LinearOperator class: Jean-Paul Pelteret, 2015 */ template class PointerMatrixAux : public PointerMatrixBase { public: + /** + * Declare type for container size. + */ + typedef typename PointerMatrixBase::size_type size_type; + /** * Constructor. The pointer in the argument is stored in this class. As * usual, the lifetime of *M must be longer than the one of the @@ -253,6 +300,12 @@ public: const MATRIX *M, const char *name); + /** + * Destructor + */ + virtual ~PointerMatrixAux () + {} + // Use doc from base class virtual void clear(); @@ -297,6 +350,18 @@ public: virtual void Tvmult_add (VECTOR &dst, const VECTOR &src) const; + /** + * Return the dimension of the codomain (or range) space. To remember: the + * matrix is of dimension $m \times n$. + */ + virtual size_type m () const; + + /** + * Return the dimension of the domain space. To remember: the matrix is of + * dimension $m \times n$. + */ + virtual size_type n () const; + private: /** * The backup memory if none was provided. @@ -311,7 +376,7 @@ private: /** * The pointer to the actual matrix. */ - SmartPointer > m; + SmartPointer > mtrx; }; @@ -324,12 +389,18 @@ private: * (#vmult()) and scalar multiplication (#Tvmult()) functions of the Vector * class. * - * @author Guido Kanschat, 2006 + * @author Guido Kanschat, 2006; extension for full compatibility with + * LinearOperator class: Jean-Paul Pelteret, 2015 */ template class PointerMatrixVector : public PointerMatrixBase > { public: + /** + * Declare type for container size. + */ + typedef typename PointerMatrixBase< Vector >::size_type size_type; + /** * Constructor. The pointer in the argument is stored in this class. As * usual, the lifetime of *M must be longer than the one of the @@ -365,6 +436,12 @@ public: PointerMatrixVector (const Vector *M, const char *name); + /** + * Destructor + */ + virtual ~PointerMatrixVector () + {} + // Use doc from base class virtual void clear(); @@ -418,11 +495,23 @@ public: virtual void Tvmult_add (Vector &dst, const Vector &src) const; + /** + * Return the dimension of the codomain (or range) space. To remember: the + * matrix is of dimension $m \times n$. + */ + virtual size_type m () const; + + /** + * Return the dimension of the domain space. To remember: the matrix is of + * dimension $m \times n$. + */ + virtual size_type n () const; + private: /** * The pointer to the actual matrix. */ - SmartPointer,PointerMatrixVector > m; + SmartPointer,PointerMatrixVector > mtrx; }; @@ -585,13 +674,13 @@ PointerMatrixBase::~PointerMatrixBase () template PointerMatrix::PointerMatrix (const MATRIX *M) - : m(M, typeid(*this).name()) + : mtrx(M, typeid(*this).name()) {} template PointerMatrix::PointerMatrix (const char *name) - : m(0, name) + : mtrx(0, name) {} @@ -599,7 +688,7 @@ template PointerMatrix::PointerMatrix ( const MATRIX *M, const char *name) - : m(M, name) + : mtrx(M, name) {} @@ -607,7 +696,7 @@ template inline void PointerMatrix::clear () { - m = 0; + mtrx = 0; } @@ -615,7 +704,7 @@ template inline const PointerMatrix & PointerMatrix::operator= (const MATRIX *M) { - m = M; + mtrx = M; return *this; } @@ -624,9 +713,9 @@ template inline bool PointerMatrix::empty () const { - if (m == 0) + if (mtrx == 0) return true; - return m->empty(); + return mtrx->empty(); } template @@ -634,8 +723,8 @@ inline void PointerMatrix::vmult (VECTOR &dst, const VECTOR &src) const { - Assert (m != 0, ExcNotInitialized()); - m->vmult (dst, src); + Assert (mtrx != 0, ExcNotInitialized()); + mtrx->vmult (dst, src); } @@ -644,8 +733,8 @@ inline void PointerMatrix::Tvmult (VECTOR &dst, const VECTOR &src) const { - Assert (m != 0, ExcNotInitialized()); - m->Tvmult (dst, src); + Assert (mtrx != 0, ExcNotInitialized()); + mtrx->Tvmult (dst, src); } @@ -654,8 +743,8 @@ inline void PointerMatrix::vmult_add (VECTOR &dst, const VECTOR &src) const { - Assert (m != 0, ExcNotInitialized()); - m->vmult_add (dst, src); + Assert (mtrx != 0, ExcNotInitialized()); + mtrx->vmult_add (dst, src); } @@ -664,11 +753,28 @@ inline void PointerMatrix::Tvmult_add (VECTOR &dst, const VECTOR &src) const { - Assert (m != 0, ExcNotInitialized()); - m->Tvmult_add (dst, src); + Assert (mtrx != 0, ExcNotInitialized()); + mtrx->Tvmult_add (dst, src); } +template +inline typename PointerMatrix::size_type +PointerMatrix::m () const +{ + Assert (mtrx != 0, ExcNotInitialized()); + return mtrx->m(); +} + + +template +inline typename PointerMatrix::size_type +PointerMatrix::n () const +{ + Assert (mtrx != 0, ExcNotInitialized()); + return mtrx->n(); +} + //----------------------------------------------------------------------// @@ -678,7 +784,7 @@ PointerMatrixAux::PointerMatrixAux ( VectorMemory *mem, const MATRIX *M) : mem(mem, typeid(*this).name()), - m(M, typeid(*this).name()) + mtrx(M, typeid(*this).name()) { if (mem == 0) mem = &my_memory; } @@ -689,7 +795,7 @@ PointerMatrixAux::PointerMatrixAux ( VectorMemory *mem, const char *name) : mem(mem, name), - m(0, name) + mtrx(0, name) { if (mem == 0) mem = &my_memory; } @@ -701,7 +807,7 @@ PointerMatrixAux::PointerMatrixAux ( const MATRIX *M, const char *name) : mem(mem, name), - m(M, name) + mtrx(M, name) { if (mem == 0) mem = &my_memory; } @@ -711,7 +817,7 @@ template inline void PointerMatrixAux::clear () { - m = 0; + mtrx = 0; } @@ -719,7 +825,7 @@ template inline const PointerMatrixAux & PointerMatrixAux::operator= (const MATRIX *M) { - m = M; + mtrx = M; return *this; } @@ -738,9 +844,9 @@ template inline bool PointerMatrixAux::empty () const { - if (m == 0) + if (mtrx == 0) return true; - return m->empty(); + return mtrx->empty(); } template @@ -751,8 +857,8 @@ PointerMatrixAux::vmult (VECTOR &dst, if (mem == 0) mem = &my_memory; Assert (mem != 0, ExcNotInitialized()); - Assert (m != 0, ExcNotInitialized()); - m->vmult (dst, src); + Assert (mtrx != 0, ExcNotInitialized()); + mtrx->vmult (dst, src); } @@ -764,8 +870,8 @@ PointerMatrixAux::Tvmult (VECTOR &dst, if (mem == 0) mem = &my_memory; Assert (mem != 0, ExcNotInitialized()); - Assert (m != 0, ExcNotInitialized()); - m->Tvmult (dst, src); + Assert (mtrx != 0, ExcNotInitialized()); + mtrx->Tvmult (dst, src); } @@ -777,10 +883,10 @@ PointerMatrixAux::vmult_add (VECTOR &dst, if (mem == 0) mem = &my_memory; Assert (mem != 0, ExcNotInitialized()); - Assert (m != 0, ExcNotInitialized()); + Assert (mtrx != 0, ExcNotInitialized()); VECTOR *v = mem->alloc(); v->reinit(dst); - m->vmult (*v, src); + mtrx->vmult (*v, src); dst += *v; mem->free(v); } @@ -794,27 +900,45 @@ PointerMatrixAux::Tvmult_add (VECTOR &dst, if (mem == 0) mem = &my_memory; Assert (mem != 0, ExcNotInitialized()); - Assert (m != 0, ExcNotInitialized()); + Assert (mtrx != 0, ExcNotInitialized()); VECTOR *v = mem->alloc(); v->reinit(dst); - m->Tvmult (*v, src); + mtrx->Tvmult (*v, src); dst += *v; mem->free(v); } +template +inline typename PointerMatrixAux::size_type +PointerMatrixAux::m () const +{ + Assert (mtrx != 0, ExcNotInitialized()); + return mtrx->m(); +} + + +template +inline typename PointerMatrixAux::size_type +PointerMatrixAux::n () const +{ + Assert (mtrx != 0, ExcNotInitialized()); + return mtrx->n(); +} + + //----------------------------------------------------------------------// template PointerMatrixVector::PointerMatrixVector (const Vector *M) - : m(M, typeid(*this).name()) + : mtrx(M, typeid(*this).name()) {} template PointerMatrixVector::PointerMatrixVector (const char *name) - : m(0, name) + : mtrx(0, name) {} @@ -822,7 +946,7 @@ template PointerMatrixVector::PointerMatrixVector ( const Vector *M, const char *name) - : m(M, name) + : mtrx(M, name) {} @@ -830,7 +954,7 @@ template inline void PointerMatrixVector::clear () { - m = 0; + mtrx = 0; } @@ -838,7 +962,7 @@ template inline const PointerMatrixVector & PointerMatrixVector::operator= (const Vector *M) { - m = M; + mtrx = M; return *this; } @@ -847,9 +971,9 @@ template inline bool PointerMatrixVector::empty () const { - if (m == 0) + if (mtrx == 0) return true; - return m->empty(); + return mtrx->empty(); } template @@ -858,10 +982,10 @@ PointerMatrixVector::vmult ( Vector &dst, const Vector &src) const { - Assert (m != 0, ExcNotInitialized()); + Assert (mtrx != 0, ExcNotInitialized()); Assert (dst.size() == 1, ExcDimensionMismatch(dst.size(), 1)); - dst(0) = *m * src; + dst(0) = *mtrx * src; } @@ -871,10 +995,10 @@ PointerMatrixVector::Tvmult ( Vector &dst, const Vector &src) const { - Assert (m != 0, ExcNotInitialized()); + Assert (mtrx != 0, ExcNotInitialized()); Assert(src.size() == 1, ExcDimensionMismatch(src.size(), 1)); - dst.equ (src(0), *m); + dst.equ (src(0), *mtrx); } @@ -884,10 +1008,10 @@ PointerMatrixVector::vmult_add ( Vector &dst, const Vector &src) const { - Assert (m != 0, ExcNotInitialized()); + Assert (mtrx != 0, ExcNotInitialized()); Assert (dst.size() == 1, ExcDimensionMismatch(dst.size(), 1)); - dst(0) += *m * src; + dst(0) += *mtrx * src; } @@ -897,10 +1021,28 @@ PointerMatrixVector::Tvmult_add ( Vector &dst, const Vector &src) const { - Assert (m != 0, ExcNotInitialized()); + Assert (mtrx != 0, ExcNotInitialized()); Assert(src.size() == 1, ExcDimensionMismatch(src.size(), 1)); - dst.add (src(0), *m); + dst.add (src(0), *mtrx); +} + + +template +inline typename PointerMatrixVector::size_type +PointerMatrixVector::m () const +{ + Assert (mtrx != 0, ExcNotInitialized()); + return mtrx->m(); +} + + +template +inline typename PointerMatrixVector::size_type +PointerMatrixVector::n () const +{ + Assert (mtrx != 0, ExcNotInitialized()); + return mtrx->n(); } diff --git a/include/deal.II/lac/precondition.h b/include/deal.II/lac/precondition.h index c21ead9852..246af14d9c 100644 --- a/include/deal.II/lac/precondition.h +++ b/include/deal.II/lac/precondition.h @@ -66,11 +66,17 @@ namespace parallel * Alternatively, the IdentityMatrix class can be used to precondition in this * way. * - * @author Guido Kanschat, 1999 + * @author Guido Kanschat, 1999; extension for full compatibility with + * LinearOperator class: Jean-Paul Pelteret, 2015 */ class PreconditionIdentity : public Subscriptor { public: + /** + * Declare type for container size. + */ + typedef types::global_dof_index size_type; + /** * This function is only present to provide the interface of a * preconditioner to be handed to a smoother. This does nothing. @@ -83,6 +89,11 @@ public: AdditionalData () {} }; + /** + * Constructor, sets the domain and range sizes to their defaults. + */ + PreconditionIdentity(); + /** * The matrix argument is ignored and here just for compatibility with more * complex preconditioners. @@ -122,6 +133,35 @@ public: * preconditioner to be handed to a smoother. This does nothing. */ void clear () {} + + /** + * Return the dimension of the codomain (or range) space. To remember: the + * matrix is of dimension $m \times n$. + * + * @note This function should only be called if the preconditioner has been + * initialized. + */ + size_type m () const; + + /** + * Return the dimension of the domain space. To remember: the matrix is of + * dimension $m \times n$. + * + * @note This function should only be called if the preconditioner has been + * initialized. + */ + size_type n () const; + +private: + /** + * The dimension of the range space. + */ + size_type _m; + + /** + * The dimension of the domain space. + */ + size_type _n; }; @@ -136,11 +176,17 @@ public: * multiplied. Still, this class is useful in multigrid smoother objects * (MGSmootherRelaxation). * - * @author Guido Kanschat, 2005 + * @author Guido Kanschat, 2005; extension for full compatibility with + * LinearOperator class: Jean-Paul Pelteret, 2015 */ class PreconditionRichardson : public Subscriptor { public: + /** + * Declare type for container size. + */ + typedef types::global_dof_index size_type; + /** * Parameters for Richardson preconditioner. */ @@ -160,7 +206,8 @@ public: }; /** - * Constructor, sets the relaxation parameter to its default. + * Constructor, sets the relaxation parameter, domain and range sizes + * to their default. */ PreconditionRichardson(); @@ -175,7 +222,7 @@ public: * compatibility with more complex preconditioners. */ template - void initialize (const MATRIX &, + void initialize (const MATRIX &matrix, const AdditionalData ¶meters); /** @@ -209,11 +256,39 @@ public: */ void clear () {} + /** + * Return the dimension of the codomain (or range) space. To remember: the + * matrix is of dimension $m \times n$. + * + * @note This function should only be called if the preconditioner has been + * initialized. + */ + size_type m () const; + + /** + * Return the dimension of the domain space. To remember: the matrix is of + * dimension $m \times n$. + * + * @note This function should only be called if the preconditioner has been + * initialized. + */ + size_type n () const; + private: /** * The relaxation parameter multiplied with the vectors. */ double relaxation; + + /** + * The dimension of the range space. + */ + size_type _m; + + /** + * The dimension of the domain space. + */ + size_type _n; }; @@ -301,12 +376,18 @@ private: * Jacobi, SOR and SSOR preconditioners are implemented. For preconditioning, * refer to derived classes. * - * @author Guido Kanschat, 2000 + * @author Guido Kanschat, 2000; extension for full compatibility with + * LinearOperator class: Jean-Paul Pelteret, 2015 */ template > class PreconditionRelaxation : public Subscriptor { public: + /** + * Declare type for container size. + */ + typedef typename MATRIX::size_type size_type; + /** * Class for parameters. */ @@ -337,6 +418,18 @@ public: */ void clear(); + /** + * Return the dimension of the codomain (or range) space. To remember: the + * matrix is of dimension $m \times n$. + */ + size_type m () const; + + /** + * Return the dimension of the domain space. To remember: the matrix is of + * dimension $m \times n$. + */ + size_type n () const; + protected: /** * Pointer to the matrix object. @@ -514,7 +607,7 @@ public: /** * Declare type for container size. */ - typedef types::global_dof_index size_type; + typedef typename MATRIX::size_type size_type; /** * A typedef to the base class. @@ -594,7 +687,8 @@ private: * solver.solve (A, x, b, precondition); * @endcode * - * @author Guido Kanschat, 2003 + * @author Guido Kanschat, 2003; extension for full compatibility with + * LinearOperator class: Jean-Paul Pelteret, 2015 */ template > class PreconditionPSOR : public PreconditionRelaxation @@ -603,7 +697,42 @@ public: /** * Declare type for container size. */ - typedef types::global_dof_index size_type; + typedef typename MATRIX::size_type size_type; + + /** + * Parameters for PreconditionPSOR. + */ + class AdditionalData + { + public: + /** + * Constructor. For the parameters' description, see below. + * + * The permutation vectors are stored as a reference. 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. + */ + AdditionalData (const std::vector &permutation, + const std::vector &inverse_permutation, + const typename PreconditionRelaxation::AdditionalData + ¶meters = typename PreconditionRelaxation::AdditionalData()); + + /** + * Storage for the permutation vector. + */ + const std::vector &permutation; + /** + * Storage for the inverse permutation vector. + */ + const std::vector &inverse_permutation; + /** + * Relaxation parameters + */ + typename PreconditionRelaxation::AdditionalData parameters; + }; /** * Initialize matrix and relaxation parameter. The matrix is just stored in @@ -622,6 +751,19 @@ public: const typename PreconditionRelaxation::AdditionalData & parameters = typename PreconditionRelaxation::AdditionalData()); + /** + * Initialize matrix and relaxation parameter. The matrix is just stored in + * the preconditioner object. + * + * For more detail about possible parameters, see the class documentation + * and the documentation of the PreconditionPSOR::AdditionalData class. + * + * After this function is called the preconditioner is ready to be used + * (using the vmult function of derived classes). + */ + void initialize (const MATRIX &A, + const AdditionalData &additional_data); + /** * Apply preconditioner. */ @@ -658,7 +800,8 @@ private: * This class is useful e.g. in multigrid smoother objects, since it is * trivially %parallel (assuming that matrix-vector products are %parallel). * - * @author Martin Kronbichler, 2009 + * @author Martin Kronbichler, 2009; extension for full compatibility with + * LinearOperator class: Jean-Paul Pelteret, 2015 */ template , class VECTOR=Vector > class PreconditionChebyshev : public Subscriptor @@ -667,7 +810,7 @@ public: /** * Declare type for container size. */ - typedef types::global_dof_index size_type; + typedef typename MATRIX::size_type size_type; /** * Standardized data struct to pipe additional parameters to the @@ -778,6 +921,18 @@ public: */ void clear (); + /** + * Return the dimension of the codomain (or range) space. To remember: the + * matrix is of dimension $m \times n$. + */ + size_type m () const; + + /** + * Return the dimension of the domain space. To remember: the matrix is of + * dimension $m \times n$. + */ + size_type n () const; + private: /** @@ -824,12 +979,22 @@ private: #ifndef DOXYGEN +inline +PreconditionIdentity::PreconditionIdentity () + : + _m (0), + _n (0) +{} + template inline void PreconditionIdentity::initialize ( - const MATRIX &, + const MATRIX &matrix, const PreconditionIdentity::AdditionalData &) -{} +{ + _m = matrix.m(); + _n = matrix.n(); +} template @@ -864,6 +1029,20 @@ PreconditionIdentity::Tvmult_add (VECTOR &dst, const VECTOR &src) const dst.add(src); } +inline typename PreconditionIdentity::size_type +PreconditionIdentity::m () const +{ + Assert(_m != 0, ExcNotInitialized()); + return _m; +} + +inline typename PreconditionIdentity::size_type +PreconditionIdentity::n () const +{ + Assert(_n != 0, ExcNotInitialized()); + return _n; +} + //--------------------------------------------------------------------------- inline @@ -877,7 +1056,9 @@ PreconditionRichardson::AdditionalData::AdditionalData ( inline PreconditionRichardson::PreconditionRichardson () : - relaxation(0) + relaxation(0), + _m (0), + _n (0) { AdditionalData add_data; relaxation=add_data.relaxation; @@ -897,10 +1078,12 @@ PreconditionRichardson::initialize ( template inline void PreconditionRichardson::initialize ( - const MATRIX &, + const MATRIX &matrix, const PreconditionRichardson::AdditionalData ¶meters) { relaxation = parameters.relaxation; + _m = matrix.m(); + _n = matrix.n(); } @@ -937,6 +1120,20 @@ PreconditionRichardson::Tvmult_add (VECTOR &dst, const VECTOR &src) const dst.add(relaxation,src); } +inline typename PreconditionRichardson::size_type +PreconditionRichardson::m () const +{ + Assert(_m != 0, ExcNotInitialized()); + return _m; +} + +inline typename PreconditionRichardson::size_type +PreconditionRichardson::n () const +{ + Assert(_n != 0, ExcNotInitialized()); + return _n; +} + //--------------------------------------------------------------------------- template @@ -956,6 +1153,21 @@ PreconditionRelaxation::clear () A = 0; } +template +inline typename PreconditionRelaxation::size_type +PreconditionRelaxation::m () const +{ + Assert (A!=0, ExcNotInitialized()); + return A->m(); +} + +template +inline typename PreconditionRelaxation::size_type +PreconditionRelaxation::n () const +{ + Assert (A!=0, ExcNotInitialized()); + return A->n(); +} //--------------------------------------------------------------------------- @@ -1143,6 +1355,19 @@ PreconditionPSOR::initialize ( } +template +inline void +PreconditionPSOR::initialize ( + const MATRIX &A, + const AdditionalData &additional_data) +{ + initialize(A, + additional_data.permutation, + additional_data.inverse_permutation, + additional_data.parameters); +} + + template template inline void @@ -1165,6 +1390,19 @@ PreconditionPSOR::Tvmult (VECTOR &dst, const VECTOR &src) const this->A->TPSOR (dst, *permutation, *inverse_permutation, this->relaxation); } +template +PreconditionPSOR::AdditionalData::AdditionalData ( + const std::vector &permutation, + const std::vector &inverse_permutation, + const typename PreconditionRelaxation::AdditionalData ¶meters) + : + permutation(permutation), + inverse_permutation(inverse_permutation), + parameters(parameters) +{ + +} + //--------------------------------------------------------------------------- @@ -1620,6 +1858,24 @@ void PreconditionChebyshev::clear () } +template +inline +typename PreconditionChebyshev::size_type +PreconditionChebyshev::m () const +{ + Assert (matrix_ptr!=0, ExcNotInitialized()); + return matrix_ptr->m(); +} + + +template +inline +typename PreconditionChebyshev::size_type +PreconditionChebyshev::n () const +{ + Assert (matrix_ptr!=0, ExcNotInitialized()); + return matrix_ptr->n(); +} #endif // DOXYGEN diff --git a/include/deal.II/lac/precondition_selector.h b/include/deal.II/lac/precondition_selector.h index 3357f2e795..3cdcfe857b 100644 --- a/include/deal.II/lac/precondition_selector.h +++ b/include/deal.II/lac/precondition_selector.h @@ -87,13 +87,18 @@ template class SparseMatrix; * started (that is several times e.g. in a nonlinear iteration) this * preselected solver and preconditioner is called. * - * @author Ralf Hartmann, 1999 + * @author Ralf Hartmann, 1999; extension for full compatibility with + * LinearOperator class: Jean-Paul Pelteret, 2015 */ template , class VECTOR = dealii::Vector > class PreconditionSelector : public Subscriptor { public: + /** + * Declare type for container size. + */ + typedef typename MATRIX::size_type size_type; /** * Constructor. @p omega denotes the damping parameter of the @@ -113,12 +118,30 @@ public: */ void use_matrix(const MATRIX &M); + /** + * Return the dimension of the codomain (or range) space. To remember: the + * matrix is of dimension $m \times n$. + */ + size_type m () const; + + /** + * Return the dimension of the domain space. To remember: the matrix is of + * dimension $m \times n$. + */ + size_type n () const; + /** * Precondition procedure. Calls the preconditioning that was specified in * the constructor. */ virtual void vmult (VECTOR &dst, const VECTOR &src) const; + /** + * Transpose precondition procedure. Calls the preconditioning that was + * specified in the constructor. + */ + virtual void Tvmult (VECTOR &dst, const VECTOR &src) const; + /** * Get the names of all implemented preconditionings. */ @@ -182,6 +205,25 @@ void PreconditionSelector::use_matrix(const MATRIX &M) A=&M; } + +template +inline typename PreconditionSelector::size_type +PreconditionSelector::m () const +{ + Assert(A!=0, ExcNoMatrixGivenToUse()); + return A->m(); +} + + +template +inline typename PreconditionSelector::size_type +PreconditionSelector::n () const +{ + Assert(A!=0, ExcNoMatrixGivenToUse()); + return A->n(); +} + + template void PreconditionSelector::vmult (VECTOR &dst, const VECTOR &src) const @@ -212,6 +254,36 @@ void PreconditionSelector::vmult (VECTOR &dst, } +template +void PreconditionSelector::Tvmult (VECTOR &dst, + const VECTOR &src) const +{ + if (preconditioning=="none") + { + dst=src; + } + else + { + Assert(A!=0, ExcNoMatrixGivenToUse()); + + if (preconditioning=="jacobi") + { + A->precondition_Jacobi(dst,src,omega); // Symmetric operation + } + else if (preconditioning=="sor") + { + A->precondition_TSOR(dst,src,omega); + } + else if (preconditioning=="ssor") + { + A->precondition_SSOR(dst,src,omega); // Symmetric operation + } + else + Assert(false,ExcNotImplemented()); + } +} + + template std::string PreconditionSelector::get_precondition_names() { diff --git a/include/deal.II/lac/sparse_decomposition.h b/include/deal.II/lac/sparse_decomposition.h index 330bc3306c..2edfe9b99a 100644 --- a/include/deal.II/lac/sparse_decomposition.h +++ b/include/deal.II/lac/sparse_decomposition.h @@ -102,7 +102,8 @@ DEAL_II_NAMESPACE_OPEN * get_strengthen_diagonal() method. * * @author Stephen "Cheffo" Kolaroff, 2002, based on SparseILU implementation - * by Wolfgang Bangerth; unified interface: Ralf Hartmann, 2003 + * by Wolfgang Bangerth; unified interface: Ralf Hartmann, 2003; extension for + * full compatibility with LinearOperator class: Jean-Paul Pelteret, 2015 */ template class SparseLUDecomposition : protected SparseMatrix, @@ -121,7 +122,7 @@ public: /** * Declare type for container size. */ - typedef types::global_dof_index size_type; + typedef typename SparseMatrix::size_type size_type; /** * Destruction. Mark the destructor pure to ensure that this class isn't @@ -217,6 +218,42 @@ public: */ bool empty () const; + /** + * Return the dimension of the codomain (or range) space. It calls the + * inherited SparseMatrix::m() function. To remember: the matrix is + * of dimension $m \times n$. + */ + size_type m () const; + + /** + * Return the dimension of the domain space. It calls the inherited + * SparseMatrix::n() function. To remember: the matrix is of dimension + * $m \times n$. + */ + size_type n () const; + + /** + * Adding Matrix-vector multiplication. Add M*src on dst with + * M being this matrix. + * + * Source and destination must not be the same vector. + * + */ + template + void vmult_add (OutVector &dst, + const InVector &src) const; + + /** + * Adding Matrix-vector multiplication. Add MT*src to + * dst with M being this matrix. This function does the same + * as vmult_add() but takes the transposed matrix. + * + * Source and destination must not be the same vector. + */ + template + void Tvmult_add (OutVector &dst, + const InVector &src) const; + /** * Determine an estimate for the memory consumption (in bytes) of this * object. @@ -317,6 +354,52 @@ SparseLUDecomposition::empty () const } +template +inline typename SparseLUDecomposition::size_type +SparseLUDecomposition::m () const +{ + return SparseMatrix::m(); +} + + +template +inline typename SparseLUDecomposition::size_type +SparseLUDecomposition::n () const +{ + return SparseMatrix::n(); +} + +// Note: This function is required for full compatibility with +// the LinearOperator class. ::MatrixInterfaceWithVmultAdd +// picks up the vmult_add function in the protected SparseMatrix +// base class. +template +template +inline void +SparseLUDecomposition::vmult_add (OutVector &dst, + const InVector &src) const +{ + OutVector tmp; + tmp.reinit(dst); + this->vmult(tmp, src); + dst += tmp; +} + +// Note: This function is required for full compatibility with +// the LinearOperator class. ::MatrixInterfaceWithVmultAdd +// picks up the vmult_add function in the protected SparseMatrix +// base class. +template +template +inline void +SparseLUDecomposition::Tvmult_add (OutVector &dst, + const InVector &src) const +{ + OutVector tmp; + tmp.reinit(dst); + this->Tvmult(tmp, src); + dst += tmp; +} //--------------------------------------------------------------------------- diff --git a/include/deal.II/lac/sparse_direct.h b/include/deal.II/lac/sparse_direct.h index 837e9260d0..d01450f5e9 100644 --- a/include/deal.II/lac/sparse_direct.h +++ b/include/deal.II/lac/sparse_direct.h @@ -74,7 +74,8 @@ DEAL_II_NAMESPACE_OPEN * * @ingroup Solvers Preconditioners * - * @author Wolfgang Bangerth, 2004 + * @author Wolfgang Bangerth, 2004; extension for full compatibility + * with LinearOperator class: Jean-Paul Pelteret, 2015 */ class SparseDirectUMFPACK : public Subscriptor { @@ -185,6 +186,18 @@ public: void Tvmult (BlockVector &dst, const BlockVector &src) const; + /** + * Return the dimension of the codomain (or range) space. To remember: the + * matrix is of dimension $m \times n$. + */ + size_type m () const; + + /** + * Return the dimension of the domain space. To remember: the matrix is of + * dimension $m \times n$. + */ + size_type n () const; + /** * @} */ @@ -281,6 +294,16 @@ public: "UMFPACK.")); private: + /** + * The dimension of the range space. + */ + size_type _m; + + /** + * The dimension of the domain space. + */ + size_type _n; + /** * The UMFPACK routines allocate objects in which they store information * about symbolic and numeric values of the decomposition. The actual data diff --git a/include/deal.II/lac/sparse_ilu.h b/include/deal.II/lac/sparse_ilu.h index 9c7b68a6e8..53c342b9fc 100644 --- a/include/deal.II/lac/sparse_ilu.h +++ b/include/deal.II/lac/sparse_ilu.h @@ -63,7 +63,7 @@ public: /** * Declare type for container size. */ - typedef types::global_dof_index size_type; + typedef typename SparseLUDecomposition::size_type size_type; /** * Constructor. Does nothing. diff --git a/include/deal.II/lac/sparse_mic.h b/include/deal.II/lac/sparse_mic.h index 4ecd6095e7..91b4d8c1b8 100644 --- a/include/deal.II/lac/sparse_mic.h +++ b/include/deal.II/lac/sparse_mic.h @@ -40,7 +40,8 @@ DEAL_II_NAMESPACE_OPEN * by the condition $\text{rowsum}(A) = \text{rowsum}(B)$. * * @author Stephen "Cheffo" Kolaroff, 2002, unified interface: Ralf Hartmann - * 2003. + * 2003; extension for full compatibility with LinearOperator class: Jean-Paul + * Pelteret, 2015. */ template class SparseMIC : public SparseLUDecomposition @@ -107,6 +108,19 @@ public: void vmult (Vector &dst, const Vector &src) const; + /** + * Apply the transpose of the incomplete decomposition, i.e. do one forward- + * backward step $dst=(LU)^{-1}src$. + * + * Call @p initialize before calling this function. + * + * @note This function has not yet been implemented + * + */ + template + void Tvmult (Vector &dst, + const Vector &src) const; + /** * Determine an estimate for the memory consumption (in bytes) of this * object. diff --git a/include/deal.II/lac/sparse_mic.templates.h b/include/deal.II/lac/sparse_mic.templates.h index 59b939405c..46df280943 100644 --- a/include/deal.II/lac/sparse_mic.templates.h +++ b/include/deal.II/lac/sparse_mic.templates.h @@ -187,6 +187,17 @@ SparseMIC::vmult (Vector &dst, } +// Exists for full compatibility with the LinearOperator class +template +template +void +SparseMIC::Tvmult (Vector &dst, + const Vector &src) const +{ + AssertThrow(false, ExcNotImplemented()); +} + + template std::size_t diff --git a/include/deal.II/lac/sparse_vanka.h b/include/deal.II/lac/sparse_vanka.h index 41799e942d..b9eb232347 100644 --- a/include/deal.II/lac/sparse_vanka.h +++ b/include/deal.II/lac/sparse_vanka.h @@ -124,7 +124,8 @@ template class SparseBlockVanka; * @ref Instantiations * in the manual). * - * @author Guido Kanschat, Wolfgang Bangerth; 1999, 2000 + * @author Guido Kanschat, Wolfgang Bangerth; 1999, 2000; extension for full + * compatibility with LinearOperator class: Jean-Paul Pelteret, 2015 */ template class SparseVanka @@ -135,6 +136,14 @@ public: */ typedef types::global_dof_index size_type; + /** + * Constructor. Does nothing. + * + * Call the initialize() function before using this object as preconditioner + * (vmult()). + */ + SparseVanka (); + /** * Constructor. Gets the matrix for preconditioning and a bit vector with * entries @p true for all rows to be updated. A reference to this vector @@ -170,6 +179,50 @@ public: */ ~SparseVanka(); + /** + * Parameters for SparseVanka. + */ + class AdditionalData + { + public: + /** + * Constructor. For the parameters' description, see below. + */ + AdditionalData (const std::vector &selected, + const bool conserve_memory = false, + const unsigned int n_threads = MultithreadInfo::n_threads()); + + /** + * Indices of those degrees of freedom that we shall work on. + */ + const std::vector &selected; + + /** + * Conserve memory flag. + */ + const bool conserve_mem; + + /** + * Number of threads to be used when building the inverses. Only relevant in + * multithreaded mode. + */ + const unsigned int n_threads; + }; + + + /** + * If the default constructor is used then this function needs to be called + * before an object of this class is used as preconditioner. + * + * For more detail about possible parameters, see the class documentation + * and the documentation of the SparseVanka::AdditionalData class. + * + * After this function is called the preconditioner is ready to be used + * (using the vmult function of derived classes). + */ + void initialize (const SparseMatrix &M, + const AdditionalData &additional_data); + /** * Do the preconditioning. This function takes the residual in @p src and * returns the resulting update vector in @p dst. @@ -178,6 +231,32 @@ public: void vmult (Vector &dst, const Vector &src) const; + /** + * Apply transpose preconditioner. This function takes the residual in + * @p src and returns the resulting update vector in @p dst. + */ + template + void Tvmult (Vector &dst, + const Vector &src) const; + + /** + * Return the dimension of the codomain (or range) space. To remember: the + * matrix is of dimension $m \times n$. + * + * @note This function should only be called if the preconditioner has been + * initialized. + */ + size_type m () const; + + /** + * Return the dimension of the domain space. To remember: the matrix is of + * dimension $m \times n$. + * + * @note This function should only be called if the preconditioner has been + * initialized. + */ + size_type n () const; + protected: /** * Apply the inverses corresponding to those degrees of freedom that have a @@ -220,18 +299,18 @@ private: /** * Conserve memory flag. */ - const bool conserve_mem; + bool conserve_mem; /** * Indices of those degrees of freedom that we shall work on. */ - const std::vector &selected; + const std::vector *selected; /** * Number of threads to be used when building the inverses. Only relevant in * multithreaded mode. */ - const unsigned int n_threads; + unsigned int n_threads; /** * Array of inverse matrices, one for each degree of freedom. Only those @@ -239,6 +318,16 @@ private: */ mutable std::vector,SparseVanka > > inverses; + /** + * The dimension of the range space. + */ + size_type _m; + + /** + * The dimension of the domain space. + */ + size_type _n; + /** * Compute the inverses of all selected diagonal elements. */ @@ -479,6 +568,37 @@ private: }; /*@}*/ +/* ---------------------------------- Inline functions ------------------- */ + +#ifndef DOXYGEN + +template +inline typename SparseVanka::size_type +SparseVanka::m () const +{ + Assert(_m != 0, ExcNotInitialized()); + return _m; +} + +template +inline typename SparseVanka::size_type +SparseVanka::n () const +{ + Assert(_n != 0, ExcNotInitialized()); + return _n; +} + +template +template +inline void +SparseVanka::Tvmult (Vector &dst, + const Vector &src) const +{ + AssertThrow(false, ExcNotImplemented()); +} + +#endif // DOXYGEN + DEAL_II_NAMESPACE_CLOSE #endif diff --git a/include/deal.II/lac/sparse_vanka.templates.h b/include/deal.II/lac/sparse_vanka.templates.h index 141980b5df..99994602cd 100644 --- a/include/deal.II/lac/sparse_vanka.templates.h +++ b/include/deal.II/lac/sparse_vanka.templates.h @@ -29,20 +29,36 @@ DEAL_II_NAMESPACE_OPEN +template +SparseVanka::SparseVanka() + : + matrix (), + conserve_mem (false), + selected (), + n_threads (0), + inverses (), + _m (0), + _n (0) +{ + +} + template SparseVanka::SparseVanka(const SparseMatrix &M, - const std::vector &selected, + const std::vector &selected_dofs, const bool conserve_mem, const unsigned int n_threads) : matrix (&M, typeid(*this).name()), conserve_mem (conserve_mem), - selected (selected), + selected (&selected_dofs), n_threads (n_threads), - inverses (M.m(), 0) + inverses (M.m(), 0), + _m (M.m()), + _n (M.n()) { Assert (M.m() == M.n(), ExcNotQuadratic ()); - Assert (M.m() == selected.size(), ExcDimensionMismatch(M.m(), selected.size())); + Assert (M.m() == selected->size(), ExcDimensionMismatch(M.m(), selected->size())); if (conserve_mem == false) compute_inverses (); @@ -62,15 +78,38 @@ SparseVanka::~SparseVanka() } +template +void +SparseVanka::initialize(const SparseMatrix &M, + const AdditionalData &additional_data) +{ + matrix = &M; + conserve_mem = additional_data.conserve_mem; + selected = &(additional_data.selected); + n_threads = additional_data.n_threads; + inverses.resize(M.m()); + _m = M.m(); + _n = M.n(); + + Assert (M.m() == M.n(), ExcNotQuadratic ()); + Assert (M.m() == selected->size(), ExcDimensionMismatch(M.m(), selected->size())); + + if (conserve_mem == false) + compute_inverses (); +} + template void SparseVanka::compute_inverses () { + Assert(matrix != 0, ExcNotInitialized()); + Assert(selected != 0, ExcNotInitialized()); + #ifndef DEAL_II_WITH_THREADS compute_inverses (0, matrix->m()); #else - const size_type n_inverses = std::count (selected.begin(), - selected.end(), + const size_type n_inverses = std::count (selected->begin(), + selected->end(), true); const size_type n_inverses_per_thread = std::max(n_inverses / n_threads, @@ -100,7 +139,7 @@ SparseVanka::compute_inverses () for (size_type i=0; (im()) && (thread+1::compute_inverses (const size_type begin, // traverse all rows of the matrix // which are selected for (size_type row=begin; row::compute_inverse (const size_type row, std::vector &local_indices) { + Assert(matrix != 0, ExcNotInitialized()); + Assert(selected != 0, ExcNotInitialized()); + // first define an alias to the sparsity // pattern of the matrix, since this // will be used quite often @@ -184,6 +226,9 @@ void SparseVanka::vmult (Vector &dst, const Vector &src) const { + Assert(matrix != 0, ExcNotInitialized()); + Assert(selected != 0, ExcNotInitialized()); + // first set output vector to zero dst = 0; // then pass on to the function @@ -192,7 +237,6 @@ SparseVanka::vmult (Vector &dst, } - template template void @@ -236,7 +280,7 @@ SparseVanka::apply_preconditioner (Vector &dst, // which are selected const size_type n = matrix->m(); for (size_type row=0; row::memory_consumption () const { std::size_t mem = (sizeof(*this) + - MemoryConsumption::memory_consumption (selected)); + MemoryConsumption::memory_consumption (*selected)); for (size_type i=0; i::memory_consumption () const } +template +SparseVanka::AdditionalData::AdditionalData ( + const std::vector &selected, + const bool conserve_mem, + const unsigned int n_threads) + : + selected(selected), + conserve_mem (conserve_mem), + n_threads (n_threads) +{} + + +//--------------------------------------------------------------------------- template diff --git a/source/lac/sparse_direct.cc b/source/lac/sparse_direct.cc index 18a9ddf813..637eea4550 100644 --- a/source/lac/sparse_direct.cc +++ b/source/lac/sparse_direct.cc @@ -53,6 +53,8 @@ initialize (const SparsityPattern &) SparseDirectUMFPACK::SparseDirectUMFPACK () : + _m (0), + _n (0), symbolic_decomposition (0), numeric_decomposition (0), control (UMFPACK_CONTROL) @@ -206,6 +208,9 @@ factorize (const Matrix &matrix) clear (); + _m = matrix.m(); + _n = matrix.n(); + const size_type N = matrix.m(); // copy over the data from the matrix to the data structures UMFPACK @@ -468,6 +473,20 @@ SparseDirectUMFPACK::Tvmult ( this->solve(dst, /*transpose=*/ true); } +typename SparseDirectUMFPACK::size_type +SparseDirectUMFPACK::m () const +{ + Assert (_m!=0, ExcNotInitialized()); + return _m; +} + +typename SparseDirectUMFPACK::size_type +SparseDirectUMFPACK::n () const +{ + Assert (_n!=0, ExcNotInitialized()); + return _n; +} + // explicit instantiations for SparseMatrixUMFPACK #define InstantiateUMFPACK(MATRIX) \ diff --git a/source/lac/sparse_mic.cc b/source/lac/sparse_mic.cc index eb3e2e6d13..0fb5ca5d6d 100644 --- a/source/lac/sparse_mic.cc +++ b/source/lac/sparse_mic.cc @@ -25,20 +25,29 @@ template void SparseMIC::initialize (const SparseMatrix const AdditionalData &data); template void SparseMIC::vmult (Vector &, const Vector &) const; +template void SparseMIC::Tvmult (Vector &, + const Vector &) const; template void SparseMIC::initialize (const SparseMatrix &, const AdditionalData &data); template void SparseMIC::vmult (Vector &, const Vector &) const; +template void SparseMIC::Tvmult (Vector &, + const Vector &) const; template class SparseMIC; template void SparseMIC::initialize (const SparseMatrix &, const AdditionalData &data); template void SparseMIC::vmult (Vector &, const Vector &) const; +template void SparseMIC::Tvmult (Vector &, + const Vector &) const; template void SparseMIC::initialize (const SparseMatrix &, const AdditionalData &data); template void SparseMIC::vmult (Vector &, const Vector &) const; +template void SparseMIC::Tvmult (Vector &, + const Vector &) const; + DEAL_II_NAMESPACE_CLOSE diff --git a/tests/lac/linear_operator_08.cc b/tests/lac/linear_operator_08.cc new file mode 100644 index 0000000000..965fda2f21 --- /dev/null +++ b/tests/lac/linear_operator_08.cc @@ -0,0 +1,454 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2015 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +// Test internal preconditioner and solver options + +#include "../tests.h" + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +using namespace dealii; + + +template +void +test_preconditioner (const MATRIX &A, + const VECTOR &b, + const ADDITIONAL_DATA &data = ADDITIONAL_DATA()) +{ + const auto lo_A = linear_operator(A); + + PRECONDITIONER preconditioner; + preconditioner.initialize(A, data); + + SolverControl solver_control (100, 1.0e-10); + SolverCG solver (solver_control); + + // Exact inverse + const auto lo_A_inv = inverse_operator(lo_A, + solver, + preconditioner); + + const VECTOR x = lo_A_inv*b; + + // Approximate inverse + { + // Using exemplar matrix + const auto lo_A_inv_approx = linear_operator(A, preconditioner); + const VECTOR x_approx = lo_A_inv_approx*b; + } + { + // Stand-alone + const auto lo_A_inv_approx = linear_operator(preconditioner); + const VECTOR x_approx = lo_A_inv_approx*b; + } +} + +// For Vector +// Cannot use more generic function as Vector +// does not define vector_type +template +void +test_preconditioner (const SparseMatrix &A, + const Vector &b, + const typename PRECONDITIONER::AdditionalData &data + = typename PRECONDITIONER::AdditionalData ()) +{ + const auto lo_A = linear_operator(A); + + PRECONDITIONER preconditioner; + preconditioner.initialize(A, data); + + // Exact inverse + { + deallog.push("Exact inverse"); + SolverControl solver_control (100, 1.0e-10); + SolverCG< Vector > solver (solver_control); + const auto lo_A_inv = inverse_operator(lo_A, + solver, + preconditioner); + + const Vector x = lo_A_inv*b; + deallog.pop(); + } + + // Approximate inverses: + // Using an exemplar matrix + { + deallog.push("Exemplar matrix"); + const auto lo_A_inv_approx = linear_operator(A, preconditioner); + const Vector x_approx = lo_A_inv_approx*b; + deallog.pop(); + } + // Stand-alone + { + deallog.push("Stand-alone"); + const auto lo_A_inv_approx = linear_operator(preconditioner); + const Vector x_approx = lo_A_inv_approx*b; + deallog.pop(); + } +} + +template +void +test_solver (const SparseMatrix &A, + const Vector &b) +{ + // Standard solver + { + deallog.push("Standard solver"); + SolverControl solver_control (100, 1.0e-10); + SOLVER solver (solver_control); + + PreconditionJacobi< SparseMatrix > preconditioner; + preconditioner.initialize(A); + + Vector x; + x.reinit(b); + + solver.solve(A,x,b,preconditioner); + deallog.pop(); + } + + // Linear operator + { + deallog.push("Linear operator"); + const auto lo_A = linear_operator(A); + + SolverControl solver_control (100, 1.0e-10); + SOLVER solver (solver_control); + + PreconditionJacobi< SparseMatrix > preconditioner; + preconditioner.initialize(A); + + const auto lo_A_inv = inverse_operator(lo_A, + solver, + preconditioner); + const Vector x = lo_A_inv*b; + deallog.pop(); + } +} + +template +class PreconditionBlockIdentity : public BlockMatrixBase +{ +public: + + struct AdditionalData + { + AdditionalData () {} + }; + + virtual ~PreconditionBlockIdentity () + {} + + void + initialize (const BlockMatrixBase &matrix, + const AdditionalData &additional_data = AdditionalData()) + { + this->row_block_indices = matrix.get_row_indices(); + this->column_block_indices = matrix.get_column_indices(); + + this->sub_objects.reinit (matrix.n_block_rows(), + matrix.n_block_cols()); + } + + template + void + vmult(VECTOR &dst, const VECTOR &src) const + { + dst = src; + } + + template + void + Tvmult(VECTOR &dst, const VECTOR &src) const + { + dst = src; + } +}; + +// Not tested: +// PreconditionLU +// PreconditionMG +// PreconditionUseMatrix +// The following don't work as expected: vmult acts on Vectors, not BlockVectors +// PreconditionBlockJacobi +// PreconditionBlockSOR +// PreconditionBlockSSOR +// SparseBlockVanka + +int main() +{ + initlog(); + deallog.depth_console(0); + deallog << std::setprecision(10); + + // deal.II SparseMatrix + { + const unsigned int rc=10; + SparsityPattern sparsity_pattern (rc, rc, 0); + sparsity_pattern.compress(); + + SparseMatrix A (sparsity_pattern); + Vector b (rc); + for (unsigned int i=0; i < rc; ++i) + { + A.diag_element(i) = 2.0; + b(i) = i; + } + + // === PRECONDITIONERS === + deallog << "Preconditioners" << std::endl; + deallog.push("Preconditioners"); + + { + deallog << "PreconditionChebyshev" << std::endl; + typedef PreconditionChebyshev< SparseMatrix > PREC; + test_preconditioner(A, b); + } + { + deallog << "PreconditionIdentity" << std::endl; + typedef PreconditionIdentity PREC; + test_preconditioner(A, b); + } + { + deallog << "PreconditionJacobi" << std::endl; + typedef PreconditionJacobi< SparseMatrix > PREC; + test_preconditioner(A, b); + } + { + deallog << "PreconditionPSOR" << std::endl; + typedef PreconditionPSOR< SparseMatrix > PREC; + std::vector permutation(b.size()); + std::vector inverse_permutation(b.size()); + test_preconditioner(A, b, + typename PREC::AdditionalData(permutation, + inverse_permutation)); + } + { + deallog << "PreconditionRichardson" << std::endl; + typedef PreconditionRichardson PREC; + test_preconditioner(A, b); + } + { + deallog << "PreconditionSelector" << std::endl; + const auto lo_A = linear_operator(A); + + PreconditionSelector< SparseMatrix, Vector > + preconditioner ("jacobi"); + preconditioner.use_matrix(A); + + SolverControl solver_control (100, 1.0e-10); + SolverCG< Vector > solver (solver_control); + + // Exact inverse + const auto lo_A_inv = inverse_operator(lo_A, + solver, + preconditioner); + + const Vector x = lo_A_inv*b; + + // Approximate inverse + const auto lo_A_inv_approx = linear_operator(preconditioner); + const Vector x_approx = lo_A_inv_approx*b; + } + { + deallog << "PreconditionSOR" << std::endl; + typedef PreconditionSOR< SparseMatrix > PREC; + test_preconditioner(A, b); + } + { + deallog << "PreconditionSSOR" << std::endl; + typedef PreconditionSSOR< SparseMatrix > PREC; + test_preconditioner(A, b); + } + { + deallog << "SparseILU" << std::endl; + typedef SparseILU PREC; + test_preconditioner(A, b); + } + { + deallog << "SparseMIC" << std::endl; + typedef SparseMIC PREC; + test_preconditioner(A, b); + } + { + deallog << "SparseVanka" << std::endl; + typedef SparseVanka PREC; + typedef typename PREC::AdditionalData PREC_AD; + test_preconditioner(A, b, + PREC_AD(std::vector(rc, true))); + } + deallog.pop(); + + // === SOLVERS === + deallog << std::endl; + deallog << "Solvers" << std::endl; + deallog.push("Solvers"); + + { + deallog << "SolverBicgstab" << std::endl; + typedef SolverBicgstab< Vector > SLVR; + test_solver (A, b); + } + { + deallog << "SolverCG" << std::endl; + typedef SolverCG< Vector > SLVR; + test_solver (A, b); + } + { + deallog << "SolverFGMRES" << std::endl; + typedef SolverFGMRES< Vector > SLVR; + test_solver (A, b); + } + { + deallog << "SolverGMRES" << std::endl; + typedef SolverGMRES< Vector > SLVR; + test_solver (A, b); + } + { + deallog << "SolverMinRes" << std::endl; + typedef SolverMinRes< Vector > SLVR; + test_solver (A, b); + } + { + deallog << "SolverQMRS" << std::endl; + typedef SolverQMRS< Vector > SLVR; + test_solver (A, b); + } + { + deallog << "SolverRelaxation" << std::endl; + typedef SolverRelaxation< Vector > SLVR; + test_solver (A, b); + } + { + deallog << "SolverRichardson" << std::endl; + typedef SolverRichardson< Vector > SLVR; + test_solver (A, b); + } + { + deallog << "SolverSelector" << std::endl; + const auto lo_A = linear_operator(A); + + ReductionControl solver_control (10, 1.e-30, 1.e-2); + SolverSelector< Vector > solver; + solver.select("cg"); + solver.set_control(solver_control); + + PreconditionJacobi< SparseMatrix > preconditioner; + preconditioner.initialize(A); + + const auto lo_A_inv = inverse_operator(lo_A, + solver, + preconditioner); + + const Vector b (rc); + const Vector x = lo_A_inv*b; + } + { + deallog << "SparseDirectUMFPACK" << std::endl; + const auto lo_A = linear_operator(A); + + SparseDirectUMFPACK solver; + solver.initialize(A); + + const auto lo_A_inv = linear_operator(solver); + const Vector x = lo_A_inv*b; + } + { + deallog << "IterativeInverse" << std::endl; + + PreconditionJacobi< SparseMatrix > preconditioner; + preconditioner.initialize(A); + + ReductionControl solver_control (10, 1.e-30, 1.e-2); + IterativeInverse< Vector > A_inv; + A_inv.initialize(A,preconditioner); + A_inv.solver.select("cg"); + A_inv.solver.set_control(solver_control); + + const auto lo_A_inv = linear_operator(A_inv); + const Vector x = lo_A_inv*b; + } + deallog.pop(); + + + deallog << "SparseMatrix OK" << std::endl; + } + + // deal.II BlockSparseMatrix + { + const unsigned int blks=2; + const unsigned int rc=10; + BlockSparsityPattern sparsity_pattern; + { + BlockCompressedSimpleSparsityPattern csp(blks, blks); + for (unsigned int bi=0; bi A (sparsity_pattern); + BlockVector b (blks,rc); + for (unsigned int bi=0; bi > PREC; + test_preconditioner(A, b); + } + + deallog << "BlockSparseMatrix OK" << std::endl; + } + +} diff --git a/tests/lac/linear_operator_08.with_cxx11=on.output b/tests/lac/linear_operator_08.with_cxx11=on.output new file mode 100644 index 0000000000..0aeb146b82 --- /dev/null +++ b/tests/lac/linear_operator_08.with_cxx11=on.output @@ -0,0 +1,90 @@ + +DEAL::Preconditioners +DEAL:Preconditioners::PreconditionChebyshev +DEAL:Preconditioners:cg::Starting value 0.9486832981 +DEAL:Preconditioners:cg::Convergence step 1 value 0.000000000 +DEAL:Preconditioners:Exact inverse:cg::Starting value 16.88194302 +DEAL:Preconditioners:Exact inverse:cg::Convergence step 1 value 2.601335705e-15 +DEAL:Preconditioners::PreconditionIdentity +DEAL:Preconditioners:Exact inverse:cg::Starting value 16.88194302 +DEAL:Preconditioners:Exact inverse:cg::Convergence step 1 value 3.903187958e-15 +DEAL:Preconditioners::PreconditionJacobi +DEAL:Preconditioners:Exact inverse:cg::Starting value 16.88194302 +DEAL:Preconditioners:Exact inverse:cg::Convergence step 1 value 0.000000000 +DEAL:Preconditioners::PreconditionPSOR +DEAL:Preconditioners:Exact inverse:cg::Starting value 16.88194302 +DEAL:Preconditioners:Exact inverse:cg::Convergence step 1 value 0.000000000 +DEAL:Preconditioners::PreconditionRichardson +DEAL:Preconditioners:Exact inverse:cg::Starting value 16.88194302 +DEAL:Preconditioners:Exact inverse:cg::Convergence step 1 value 0.000000000 +DEAL:Preconditioners::PreconditionSelector +DEAL:Preconditioners:cg::Starting value 16.88194302 +DEAL:Preconditioners:cg::Convergence step 1 value 0.000000000 +DEAL:Preconditioners::PreconditionSOR +DEAL:Preconditioners:Exact inverse:cg::Starting value 16.88194302 +DEAL:Preconditioners:Exact inverse:cg::Convergence step 1 value 0.000000000 +DEAL:Preconditioners::PreconditionSSOR +DEAL:Preconditioners:Exact inverse:cg::Starting value 16.88194302 +DEAL:Preconditioners:Exact inverse:cg::Convergence step 1 value 0.000000000 +DEAL:Preconditioners::SparseILU +DEAL:Preconditioners:Exact inverse:cg::Starting value 16.88194302 +DEAL:Preconditioners:Exact inverse:cg::Convergence step 1 value 0.000000000 +DEAL:Preconditioners::SparseMIC +DEAL:Preconditioners:Exact inverse:cg::Starting value 16.88194302 +DEAL:Preconditioners:Exact inverse:cg::Convergence step 1 value 0.000000000 +DEAL:Preconditioners::SparseVanka +DEAL:Preconditioners:Exact inverse:cg::Starting value 16.88194302 +DEAL:Preconditioners:Exact inverse:cg::Convergence step 1 value 0.000000000 +DEAL:: +DEAL::Solvers +DEAL:Solvers::SolverBicgstab +DEAL:Solvers:Standard solver:Bicgstab::Starting value 16.88194302 +DEAL:Solvers:Standard solver:Bicgstab::Convergence step 1 value 0.000000000 +DEAL:Solvers:Linear operator:Bicgstab::Starting value 16.88194302 +DEAL:Solvers:Linear operator:Bicgstab::Convergence step 1 value 0.000000000 +DEAL:Solvers::SolverCG +DEAL:Solvers:Standard solver:cg::Starting value 16.88194302 +DEAL:Solvers:Standard solver:cg::Convergence step 1 value 0.000000000 +DEAL:Solvers:Linear operator:cg::Starting value 16.88194302 +DEAL:Solvers:Linear operator:cg::Convergence step 1 value 0.000000000 +DEAL:Solvers::SolverFGMRES +DEAL:Solvers:Standard solver:FGMRES::Starting value 16.88194302 +DEAL:Solvers:Standard solver:FGMRES::Convergence step 1 value 2.744725072e-15 +DEAL:Solvers:Linear operator:FGMRES::Starting value 16.88194302 +DEAL:Solvers:Linear operator:FGMRES::Convergence step 1 value 2.744725072e-15 +DEAL:Solvers::SolverGMRES +DEAL:Solvers:Standard solver:GMRES::Starting value 8.440971508 +DEAL:Solvers:Standard solver:GMRES::Convergence step 1 value 1.372362536e-15 +DEAL:Solvers:Linear operator:GMRES::Starting value 8.440971508 +DEAL:Solvers:Linear operator:GMRES::Convergence step 1 value 1.372362536e-15 +DEAL:Solvers::SolverMinRes +DEAL:Solvers:Standard solver:minres::Starting value 11.93733639 +DEAL:Solvers:Standard solver:minres::Convergence step 1 value 2.059168646e-15 +DEAL:Solvers:Linear operator:minres::Starting value 11.93733639 +DEAL:Solvers:Linear operator:minres::Convergence step 1 value 2.059168646e-15 +DEAL:Solvers::SolverQMRS +DEAL:Solvers:Standard solver:QMRS::Starting value 16.88194302 +DEAL:Solvers:Standard solver:QMRS::Convergence step 1 value 0.000000000 +DEAL:Solvers:Linear operator:QMRS::Starting value 16.88194302 +DEAL:Solvers:Linear operator:QMRS::Convergence step 1 value 0.000000000 +DEAL:Solvers::SolverRelaxation +DEAL:Solvers:Standard solver:Relaxation::Starting value 16.88194302 +DEAL:Solvers:Standard solver:Relaxation::Convergence step 1 value 0.000000000 +DEAL:Solvers:Linear operator:Relaxation::Starting value 16.88194302 +DEAL:Solvers:Linear operator:Relaxation::Convergence step 1 value 0.000000000 +DEAL:Solvers::SolverRichardson +DEAL:Solvers:Standard solver:Richardson::Starting value 16.88194302 +DEAL:Solvers:Standard solver:Richardson::Convergence step 1 value 0.000000000 +DEAL:Solvers:Linear operator:Richardson::Starting value 16.88194302 +DEAL:Solvers:Linear operator:Richardson::Convergence step 1 value 0.000000000 +DEAL:Solvers::SolverSelector +DEAL:Solvers:cg::Convergence step 0 value 0.000000000 +DEAL:Solvers::SparseDirectUMFPACK +DEAL:Solvers::IterativeInverse +DEAL:Solvers:cg::Starting value 16.88194302 +DEAL:Solvers:cg::Convergence step 1 value 0.000000000 +DEAL::SparseMatrix OK +DEAL::PreconditionBlockIdentity +DEAL:cg::Starting value 49.69909456 +DEAL:cg::Convergence step 1 value 0.000000000 +DEAL::BlockSparseMatrix OK