]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Ammended native deal.II preconditioners and solvers for use in LinearOperators.
authorJean-Paul Pelteret <jppelteret@gmail.com>
Thu, 24 Sep 2015 14:30:05 +0000 (16:30 +0200)
committerMatthias Maier <tamiko@43-1.org>
Fri, 23 Oct 2015 17:57:35 +0000 (12:57 -0500)
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.

17 files changed:
doc/news/changes.h
include/deal.II/lac/iterative_inverse.h
include/deal.II/lac/linear_operator.h
include/deal.II/lac/pointer_matrix.h
include/deal.II/lac/precondition.h
include/deal.II/lac/precondition_selector.h
include/deal.II/lac/sparse_decomposition.h
include/deal.II/lac/sparse_direct.h
include/deal.II/lac/sparse_ilu.h
include/deal.II/lac/sparse_mic.h
include/deal.II/lac/sparse_mic.templates.h
include/deal.II/lac/sparse_vanka.h
include/deal.II/lac/sparse_vanka.templates.h
source/lac/sparse_direct.cc
source/lac/sparse_mic.cc
tests/lac/linear_operator_08.cc [new file with mode: 0644]
tests/lac/linear_operator_08.with_cxx11=on.output [new file with mode: 0644]

index 2449e7eebb5f90f644f27de0e34902640b344606..3edcdc60fb0c358288da66ea21c924f81cf61de5 100644 (file)
@@ -170,6 +170,17 @@ inconvenience this causes.
   (Timo Heister, 2015/09/30)
   </li>
 
+  <li> 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 .
+  <br>
+  (Jean-Paul Pelteret, 2015/09/24 - 2015/10/19)
+  </li>
+
   <li> New: MGTransferPrebuilt can now be used with parallel::distributed::Vector
   and TrilinosWrappers::SparseMatrix as a transfer matrix.
   <br>
@@ -372,6 +383,13 @@ inconvenience this causes.
   (Alberto Sartori, 2015/10/22)
   </li>
 
+  <li> 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.
+  <br>
+  (Jean-Paul Pelteret, 2015/10/19)
+  </li>
+
   <li> New: Ghost cells for the multigrid levels in
   parallel::distributed::Triangulation are now correctly created also for
   periodic boundary conditions.
index dbbf17ccfc01de97251bf7f58e6d9e24b979d1c2..84375351cfdd21d65bfcba180abc62d3c71ea13b 100644 (file)
@@ -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 VECTOR>
 class IterativeInverse : public Subscriptor
 {
 public:
+  /**
+   * Declare type for container size.
+   */
+  typedef typename PointerMatrixBase<VECTOR>::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 <tt>src</tt>, 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 <class VECTOR2>
   void vmult (VECTOR2 &dst, const VECTOR2 &src) const;
 
+  /**
+   * Same as before, but uses the transpose of the matrix.
+   */
+  template <class VECTOR2>
+  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<VECTOR>::vmult (VECTOR &dst, const VECTOR &src) const
 }
 
 
+template <class VECTOR>
+inline void
+IterativeInverse<VECTOR>::Tvmult (VECTOR &dst, const VECTOR &src) const
+{
+  AssertThrow(false, ExcNotImplemented());
+}
+
+
 template <class VECTOR>
 template <class VECTOR2>
 inline void
@@ -174,6 +211,32 @@ IterativeInverse<VECTOR>::vmult (VECTOR2 &dst, const VECTOR2 &src) const
 }
 
 
+template <class VECTOR>
+template <class VECTOR2>
+inline void
+IterativeInverse<VECTOR>::Tvmult (VECTOR2 &dst, const VECTOR2 &src) const
+{
+  AssertThrow(false, ExcNotImplemented());
+}
+
+
+template <class VECTOR>
+inline typename IterativeInverse<VECTOR>::size_type
+IterativeInverse<VECTOR>::m () const
+{
+  Assert (matrix.get()!=0, ExcNotInitialized());
+  return matrix->m();
+}
+
+template <class VECTOR>
+inline typename IterativeInverse<VECTOR>::size_type
+IterativeInverse<VECTOR>::n () const
+{
+  Assert (matrix.get()!=0, ExcNotInitialized());
+  return matrix->n();
+}
+
+
 
 DEAL_II_NAMESPACE_CLOSE
 
index 2b4808d403f15af0795e783bd4da41b737b440c4..eae3309d82bd879b5c21f412125c5b9321531f8e 100644 (file)
@@ -641,6 +641,7 @@ inverse_operator(const LinearOperator<typename Solver::vector_type, typename Sol
 
   return_op.vmult = [op, &solver, &preconditioner](Vector &v, const Vector &u)
   {
+    op.reinit_range_vector(v, /*bool fast =*/ false);
     solver.solve(op, v, u, preconditioner);
   };
 
@@ -650,7 +651,7 @@ inverse_operator(const LinearOperator<typename Solver::vector_type, typename Sol
     static GrowingVectorMemory<typename Solver::vector_type> 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<typename Solver::vector_type, typename Sol
   return_op.Tvmult = [op, &solver, &preconditioner]( Vector &v, const
                                                      Vector &u)
   {
+    op.reinit_range_vector(v, /*bool fast =*/ false);
     solver.solve(transpose_operator(op), v, u, preconditioner);
   };
 
@@ -668,7 +670,7 @@ inverse_operator(const LinearOperator<typename Solver::vector_type, typename Sol
     static GrowingVectorMemory<typename Solver::vector_type> 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);
index 708f595e1cf167afeb4beac554c0638e101f426d..4b27ee6bcb5af63cbf4ecfd8b3533529974213d0 100644 (file)
@@ -51,6 +51,11 @@ template<class VECTOR>
 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 MATRIX, class VECTOR>
 class PointerMatrix : public PointerMatrixBase<VECTOR>
 {
 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 <tt>*M</tt> 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<const MATRIX,PointerMatrix<MATRIX,VECTOR> > m;
+  SmartPointer<const MATRIX,PointerMatrix<MATRIX,VECTOR> > 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 MATRIX, class VECTOR>
 class PointerMatrixAux : public PointerMatrixBase<VECTOR>
 {
 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 <tt>*M</tt> 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<const MATRIX,PointerMatrixAux<MATRIX,VECTOR> > m;
+  SmartPointer<const MATRIX,PointerMatrixAux<MATRIX,VECTOR> > 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 <typename number>
 class PointerMatrixVector : public PointerMatrixBase<Vector<number> >
 {
 public:
+  /**
+   * Declare type for container size.
+   */
+  typedef typename PointerMatrixBase< Vector<number> >::size_type size_type;
+
   /**
    * Constructor.  The pointer in the argument is stored in this class. As
    * usual, the lifetime of <tt>*M</tt> must be longer than the one of the
@@ -365,6 +436,12 @@ public:
   PointerMatrixVector (const Vector<number> *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<number> &dst,
                            const Vector<number> &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<const Vector<number>,PointerMatrixVector<number> > m;
+  SmartPointer<const Vector<number>,PointerMatrixVector<number> > mtrx;
 };
 
 
@@ -585,13 +674,13 @@ PointerMatrixBase<VECTOR>::~PointerMatrixBase ()
 
 template<class MATRIX, class VECTOR>
 PointerMatrix<MATRIX, VECTOR>::PointerMatrix (const MATRIX *M)
-  : m(M, typeid(*this).name())
+  : mtrx(M, typeid(*this).name())
 {}
 
 
 template<class MATRIX, class VECTOR>
 PointerMatrix<MATRIX, VECTOR>::PointerMatrix (const char *name)
-  : m(0, name)
+  : mtrx(0, name)
 {}
 
 
@@ -599,7 +688,7 @@ template<class MATRIX, class VECTOR>
 PointerMatrix<MATRIX, VECTOR>::PointerMatrix (
   const MATRIX *M,
   const char *name)
-  : m(M, name)
+  : mtrx(M, name)
 {}
 
 
@@ -607,7 +696,7 @@ template<class MATRIX, class VECTOR>
 inline void
 PointerMatrix<MATRIX, VECTOR>::clear ()
 {
-  m = 0;
+  mtrx = 0;
 }
 
 
@@ -615,7 +704,7 @@ template<class MATRIX, class VECTOR>
 inline const PointerMatrix<MATRIX, VECTOR> &
 PointerMatrix<MATRIX, VECTOR>::operator= (const MATRIX *M)
 {
-  m = M;
+  mtrx = M;
   return *this;
 }
 
@@ -624,9 +713,9 @@ template<class MATRIX, class VECTOR>
 inline bool
 PointerMatrix<MATRIX, VECTOR>::empty () const
 {
-  if (m == 0)
+  if (mtrx == 0)
     return true;
-  return m->empty();
+  return mtrx->empty();
 }
 
 template<class MATRIX, class VECTOR>
@@ -634,8 +723,8 @@ inline void
 PointerMatrix<MATRIX, VECTOR>::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<MATRIX, VECTOR>::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<MATRIX, VECTOR>::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<MATRIX, VECTOR>::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<class MATRIX, class VECTOR>
+inline typename PointerMatrix<MATRIX, VECTOR>::size_type
+PointerMatrix<MATRIX, VECTOR>::m () const
+{
+  Assert (mtrx != 0, ExcNotInitialized());
+  return mtrx->m();
+}
+
+
+template<class MATRIX, class VECTOR>
+inline typename PointerMatrix<MATRIX, VECTOR>::size_type
+PointerMatrix<MATRIX, VECTOR>::n () const
+{
+  Assert (mtrx != 0, ExcNotInitialized());
+  return mtrx->n();
+}
+
 
 //----------------------------------------------------------------------//
 
@@ -678,7 +784,7 @@ PointerMatrixAux<MATRIX, VECTOR>::PointerMatrixAux (
   VectorMemory<VECTOR> *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<MATRIX, VECTOR>::PointerMatrixAux (
   VectorMemory<VECTOR> *mem,
   const char *name)
   : mem(mem, name),
-    m(0, name)
+    mtrx(0, name)
 {
   if (mem == 0) mem = &my_memory;
 }
@@ -701,7 +807,7 @@ PointerMatrixAux<MATRIX, VECTOR>::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<class MATRIX, class VECTOR>
 inline void
 PointerMatrixAux<MATRIX, VECTOR>::clear ()
 {
-  m = 0;
+  mtrx = 0;
 }
 
 
@@ -719,7 +825,7 @@ template<class MATRIX, class VECTOR>
 inline const PointerMatrixAux<MATRIX, VECTOR> &
 PointerMatrixAux<MATRIX, VECTOR>::operator= (const MATRIX *M)
 {
-  m = M;
+  mtrx = M;
   return *this;
 }
 
@@ -738,9 +844,9 @@ template<class MATRIX, class VECTOR>
 inline bool
 PointerMatrixAux<MATRIX, VECTOR>::empty () const
 {
-  if (m == 0)
+  if (mtrx == 0)
     return true;
-  return m->empty();
+  return mtrx->empty();
 }
 
 template<class MATRIX, class VECTOR>
@@ -751,8 +857,8 @@ PointerMatrixAux<MATRIX, VECTOR>::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<MATRIX, VECTOR>::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<MATRIX, VECTOR>::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<MATRIX, VECTOR>::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<class MATRIX, class VECTOR>
+inline typename PointerMatrixAux<MATRIX, VECTOR>::size_type
+PointerMatrixAux<MATRIX, VECTOR>::m () const
+{
+  Assert (mtrx != 0, ExcNotInitialized());
+  return mtrx->m();
+}
+
+
+template<class MATRIX, class VECTOR>
+inline typename PointerMatrixAux<MATRIX, VECTOR>::size_type
+PointerMatrixAux<MATRIX, VECTOR>::n () const
+{
+  Assert (mtrx != 0, ExcNotInitialized());
+  return mtrx->n();
+}
+
+
 //----------------------------------------------------------------------//
 
 
 template<typename number>
 PointerMatrixVector<number>::PointerMatrixVector (const Vector<number> *M)
-  : m(M, typeid(*this).name())
+  : mtrx(M, typeid(*this).name())
 {}
 
 
 template<typename number>
 PointerMatrixVector<number>::PointerMatrixVector (const char *name)
-  : m(0, name)
+  : mtrx(0, name)
 {}
 
 
@@ -822,7 +946,7 @@ template<typename number>
 PointerMatrixVector<number>::PointerMatrixVector (
   const Vector<number> *M,
   const char *name)
-  : m(M, name)
+  : mtrx(M, name)
 {}
 
 
@@ -830,7 +954,7 @@ template<typename number>
 inline void
 PointerMatrixVector<number>::clear ()
 {
-  m = 0;
+  mtrx = 0;
 }
 
 
@@ -838,7 +962,7 @@ template<typename number>
 inline const PointerMatrixVector<number> &
 PointerMatrixVector<number>::operator= (const Vector<number> *M)
 {
-  m = M;
+  mtrx = M;
   return *this;
 }
 
@@ -847,9 +971,9 @@ template<typename number>
 inline bool
 PointerMatrixVector<number>::empty () const
 {
-  if (m == 0)
+  if (mtrx == 0)
     return true;
-  return m->empty();
+  return mtrx->empty();
 }
 
 template<typename number>
@@ -858,10 +982,10 @@ PointerMatrixVector<number>::vmult (
   Vector<number> &dst,
   const Vector<number> &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<number>::Tvmult (
   Vector<number> &dst,
   const Vector<number> &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<number>::vmult_add (
   Vector<number> &dst,
   const Vector<number> &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<number>::Tvmult_add (
   Vector<number> &dst,
   const Vector<number> &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<typename number>
+inline typename PointerMatrixVector<number>::size_type
+PointerMatrixVector<number>::m () const
+{
+  Assert (mtrx != 0, ExcNotInitialized());
+  return mtrx->m();
+}
+
+
+template<typename number>
+inline typename PointerMatrixVector<number>::size_type
+PointerMatrixVector<number>::n () const
+{
+  Assert (mtrx != 0, ExcNotInitialized());
+  return mtrx->n();
 }
 
 
index c21ead9852bbaecf1db69178c84a338e30c2203e..246af14d9cd928e9c4c906ed6a6735e7b0f5b805 100644 (file)
@@ -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 <class MATRIX>
-  void initialize (const MATRIX &,
+  void initialize (const MATRIX &matrix,
                    const AdditionalData &parameters);
 
   /**
@@ -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 MATRIX = SparseMatrix<double> >
 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 MATRIX = SparseMatrix<double> >
 class PreconditionPSOR : public PreconditionRelaxation<MATRIX>
@@ -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<size_type> &permutation,
+                    const std::vector<size_type> &inverse_permutation,
+                    const typename PreconditionRelaxation<MATRIX>::AdditionalData
+                    &parameters = typename PreconditionRelaxation<MATRIX>::AdditionalData());
+
+    /**
+     * Storage for the permutation vector.
+     */
+    const std::vector<size_type> &permutation;
+    /**
+     * Storage for the inverse permutation vector.
+     */
+    const std::vector<size_type> &inverse_permutation;
+    /**
+     * Relaxation parameters
+     */
+    typename PreconditionRelaxation<MATRIX>::AdditionalData parameters;
+  };
 
   /**
    * Initialize matrix and relaxation parameter. The matrix is just stored in
@@ -622,6 +751,19 @@ public:
                    const typename PreconditionRelaxation<MATRIX>::AdditionalData &
                    parameters = typename PreconditionRelaxation<MATRIX>::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 <code>vmult</code> 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 MATRIX=SparseMatrix<double>, class VECTOR=Vector<double> >
 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 <class MATRIX>
 inline void
 PreconditionIdentity::initialize (
-  const MATRIX &,
+  const MATRIX &matrix,
   const PreconditionIdentity::AdditionalData &)
-{}
+{
+  _m = matrix.m();
+  _n = matrix.n();
+}
 
 
 template<class VECTOR>
@@ -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 <class MATRIX>
 inline void
 PreconditionRichardson::initialize (
-  const MATRIX &,
+  const MATRIX &matrix,
   const PreconditionRichardson::AdditionalData &parameters)
 {
   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 <class MATRIX>
@@ -956,6 +1153,21 @@ PreconditionRelaxation<MATRIX>::clear ()
   A = 0;
 }
 
+template <class MATRIX>
+inline typename PreconditionRelaxation<MATRIX>::size_type
+PreconditionRelaxation<MATRIX>::m () const
+{
+  Assert (A!=0, ExcNotInitialized());
+  return A->m();
+}
+
+template <class MATRIX>
+inline typename PreconditionRelaxation<MATRIX>::size_type
+PreconditionRelaxation<MATRIX>::n () const
+{
+  Assert (A!=0, ExcNotInitialized());
+  return A->n();
+}
 
 //---------------------------------------------------------------------------
 
@@ -1143,6 +1355,19 @@ PreconditionPSOR<MATRIX>::initialize (
 }
 
 
+template <class MATRIX>
+inline void
+PreconditionPSOR<MATRIX>::initialize (
+  const MATRIX         &A,
+  const AdditionalData &additional_data)
+{
+  initialize(A,
+             additional_data.permutation,
+             additional_data.inverse_permutation,
+             additional_data.parameters);
+}
+
+
 template <class MATRIX>
 template<class VECTOR>
 inline void
@@ -1165,6 +1390,19 @@ PreconditionPSOR<MATRIX>::Tvmult (VECTOR &dst, const VECTOR &src) const
   this->A->TPSOR (dst, *permutation, *inverse_permutation, this->relaxation);
 }
 
+template <class MATRIX>
+PreconditionPSOR<MATRIX>::AdditionalData::AdditionalData (
+  const std::vector<size_type> &permutation,
+  const std::vector<size_type> &inverse_permutation,
+  const typename PreconditionRelaxation<MATRIX>::AdditionalData &parameters)
+  :
+  permutation(permutation),
+  inverse_permutation(inverse_permutation),
+  parameters(parameters)
+{
+
+}
+
 
 //---------------------------------------------------------------------------
 
@@ -1620,6 +1858,24 @@ void PreconditionChebyshev<MATRIX,VECTOR>::clear ()
 }
 
 
+template <class MATRIX, class VECTOR>
+inline
+typename PreconditionChebyshev<MATRIX,VECTOR>::size_type
+PreconditionChebyshev<MATRIX,VECTOR>::m () const
+{
+  Assert (matrix_ptr!=0, ExcNotInitialized());
+  return matrix_ptr->m();
+}
+
+
+template <class MATRIX, class VECTOR>
+inline
+typename PreconditionChebyshev<MATRIX,VECTOR>::size_type
+PreconditionChebyshev<MATRIX,VECTOR>::n () const
+{
+  Assert (matrix_ptr!=0, ExcNotInitialized());
+  return matrix_ptr->n();
+}
 
 #endif // DOXYGEN
 
index 3357f2e795ab94cad5bbb33c0d249c7e7025b5b2..3cdcfe857b5257196f6a93cf2679fda67ccf59d6 100644 (file)
@@ -87,13 +87,18 @@ template <class number> 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 MATRIX = SparseMatrix<double>,
           class VECTOR = dealii::Vector<double> >
 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<MATRIX,VECTOR>::use_matrix(const MATRIX &M)
   A=&M;
 }
 
+
+template <class MATRIX, class VECTOR>
+inline typename PreconditionSelector<MATRIX,VECTOR>::size_type
+PreconditionSelector<MATRIX,VECTOR>::m () const
+{
+  Assert(A!=0, ExcNoMatrixGivenToUse());
+  return A->m();
+}
+
+
+template <class MATRIX, class VECTOR>
+inline typename PreconditionSelector<MATRIX,VECTOR>::size_type
+PreconditionSelector<MATRIX,VECTOR>::n () const
+{
+  Assert(A!=0, ExcNoMatrixGivenToUse());
+  return A->n();
+}
+
+
 template <class MATRIX, class VECTOR>
 void PreconditionSelector<MATRIX,VECTOR>::vmult (VECTOR &dst,
                                                  const VECTOR &src) const
@@ -212,6 +254,36 @@ void PreconditionSelector<MATRIX,VECTOR>::vmult (VECTOR &dst,
 }
 
 
+template <class MATRIX, class VECTOR>
+void PreconditionSelector<MATRIX,VECTOR>::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 <class MATRIX, class VECTOR>
 std::string PreconditionSelector<MATRIX,VECTOR>::get_precondition_names()
 {
index 330bc3306cf4869e5a058d86659b78c2c5fb5e65..2edfe9b99a70ab9c658091c9413a5f80f7bcc426 100644 (file)
@@ -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 <typename number>
 class SparseLUDecomposition : protected SparseMatrix<number>,
@@ -121,7 +122,7 @@ public:
   /**
    * Declare type for container size.
    */
-  typedef types::global_dof_index size_type;
+  typedef typename SparseMatrix<number>::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 <i>M*src</i> on <i>dst</i> with
+   * <i>M</i> being this matrix.
+   *
+   * Source and destination must not be the same vector.
+   *
+   */
+  template <class OutVector, class InVector>
+  void vmult_add (OutVector &dst,
+                  const InVector &src) const;
+
+  /**
+   * Adding Matrix-vector multiplication. Add <i>M<sup>T</sup>*src</i> to
+   * <i>dst</i> with <i>M</i> 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 <class OutVector, class InVector>
+  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<number>::empty () const
 }
 
 
+template <typename number>
+inline typename SparseLUDecomposition<number>::size_type
+SparseLUDecomposition<number>::m () const
+{
+  return SparseMatrix<number>::m();
+}
+
+
+template <typename number>
+inline typename SparseLUDecomposition<number>::size_type
+SparseLUDecomposition<number>::n () const
+{
+  return SparseMatrix<number>::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 <typename number>
+template <class OutVector, class InVector>
+inline void
+SparseLUDecomposition<number>::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 <typename number>
+template <class OutVector, class InVector>
+inline void
+SparseLUDecomposition<number>::Tvmult_add (OutVector &dst,
+                                           const InVector &src) const
+{
+  OutVector tmp;
+  tmp.reinit(dst);
+  this->Tvmult(tmp, src);
+  dst += tmp;
+}
 
 //---------------------------------------------------------------------------
 
index 837e9260d06bb22c22361db01bafbdb529bcc8ac..d01450f5e977fbc3521de9160f61b85344214917 100644 (file)
@@ -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<double> &dst,
                const BlockVector<double> &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
index 9c7b68a6e8980fabda5a0c540c8d3ae97f39a5ef..53c342b9fc3e57a34f9a06ba6c57e687981be86e 100644 (file)
@@ -63,7 +63,7 @@ public:
   /**
    * Declare type for container size.
    */
-  typedef types::global_dof_index size_type;
+  typedef typename SparseLUDecomposition<number>::size_type size_type;
 
   /**
    * Constructor. Does nothing.
index 4ecd6095e7b22ff3c5aac2b05b20eca24a9acb2b..91b4d8c1b87e08c6dae2ab17a9decde4354b156a 100644 (file)
@@ -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 <typename number>
 class SparseMIC : public SparseLUDecomposition<number>
@@ -107,6 +108,19 @@ public:
   void vmult (Vector<somenumber>       &dst,
               const Vector<somenumber> &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 <typename somenumber>
+  void Tvmult (Vector<somenumber>       &dst,
+               const Vector<somenumber> &src) const;
+
   /**
    * Determine an estimate for the memory consumption (in bytes) of this
    * object.
index 59b939405c01d9a430ea0c7fb1a278f2e7a0c425..46df280943e6fea0f2329e629c846a3a7e912bf1 100644 (file)
@@ -187,6 +187,17 @@ SparseMIC<number>::vmult (Vector<somenumber>       &dst,
 }
 
 
+// Exists for full compatibility with the LinearOperator class
+template <typename number>
+template <typename somenumber>
+void
+SparseMIC<number>::Tvmult (Vector<somenumber>       &dst,
+                           const Vector<somenumber> &src) const
+{
+  AssertThrow(false, ExcNotImplemented());
+}
+
+
 
 template <typename number>
 std::size_t
index 41799e942da6dcc54e923dd19ceb8b46fb4275b5..b9eb2323475f677250e53a6df91c080ff8ec7ea6 100644 (file)
@@ -124,7 +124,8 @@ template <typename number> 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<typename number>
 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<bool> &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<bool> &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 <code>vmult</code> function of derived classes).
+   */
+  void initialize (const SparseMatrix<number> &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<number2>       &dst,
               const Vector<number2> &src) const;
 
+  /**
+   * Apply transpose preconditioner. This function takes the residual in
+   * @p src  and returns the resulting update vector in @p dst.
+   */
+  template<typename number2>
+  void Tvmult (Vector<number2>       &dst,
+               const Vector<number2> &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<bool> &selected;
+  const std::vector<bool> *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<SmartPointer<FullMatrix<float>,SparseVanka<number> > > 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<typename number>
+inline typename SparseVanka<number>::size_type
+SparseVanka<number>::m () const
+{
+  Assert(_m != 0, ExcNotInitialized());
+  return _m;
+}
+
+template<typename number>
+inline typename SparseVanka<number>::size_type
+SparseVanka<number>::n () const
+{
+  Assert(_n != 0, ExcNotInitialized());
+  return _n;
+}
+
+template<typename number>
+template<typename number2>
+inline void
+SparseVanka<number>::Tvmult (Vector<number2>       &dst,
+                             const Vector<number2> &src) const
+{
+  AssertThrow(false, ExcNotImplemented());
+}
+
+#endif // DOXYGEN
+
 DEAL_II_NAMESPACE_CLOSE
 
 #endif
index 141980b5df42510a7943eca66f7e49d8d428e3ae..99994602cd30849df655858a37e332c8ae6a080b 100644 (file)
 
 DEAL_II_NAMESPACE_OPEN
 
+template<typename number>
+SparseVanka<number>::SparseVanka()
+  :
+  matrix (),
+  conserve_mem (false),
+  selected (),
+  n_threads (0),
+  inverses (),
+  _m (0),
+  _n (0)
+{
+
+}
+
 template<typename number>
 SparseVanka<number>::SparseVanka(const SparseMatrix<number> &M,
-                                 const std::vector<bool>    &selected,
+                                 const std::vector<bool>    &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<number>::~SparseVanka()
 }
 
 
+template<typename number>
+void
+SparseVanka<number>::initialize(const SparseMatrix<number> &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 <typename number>
 void
 SparseVanka<number>::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<number>::compute_inverses ()
 
   for (size_type i=0; (i<matrix->m()) && (thread+1<n_threads); ++i)
     {
-      if (selected[i] == true)
+      if ((*selected)[i] == true)
         ++c;
       if (c == n_inverses_per_thread)
         {
@@ -141,7 +180,7 @@ SparseVanka<number>::compute_inverses (const size_type begin,
   // traverse all rows of the matrix
   // which are selected
   for (size_type row=begin; row<end; ++row)
-    if (selected[row] == true)
+    if ((*selected)[row] == true)
       compute_inverse (row, local_indices);
 }
 
@@ -152,6 +191,9 @@ void
 SparseVanka<number>::compute_inverse (const size_type         row,
                                       std::vector<size_type> &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<number>::vmult (Vector<number2>       &dst,
                             const Vector<number2> &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<number>::vmult (Vector<number2>       &dst,
 }
 
 
-
 template<typename number>
 template<typename number2>
 void
@@ -236,7 +280,7 @@ SparseVanka<number>::apply_preconditioner (Vector<number2>         &dst,
   // which are selected
   const size_type n = matrix->m();
   for (size_type row=0; row<n; ++row)
-    if ((selected[row] == true) &&
+    if (((*selected)[row] == true) &&
         ((range_is_restricted == false) || ((*dof_mask)[row] == true)))
       {
         const size_type row_length = structure.row_length(row);
@@ -356,7 +400,7 @@ std::size_t
 SparseVanka<number>::memory_consumption () const
 {
   std::size_t mem = (sizeof(*this) +
-                     MemoryConsumption::memory_consumption (selected));
+                     MemoryConsumption::memory_consumption (*selected));
   for (size_type i=0; i<inverses.size(); ++i)
     mem += MemoryConsumption::memory_consumption (*inverses[i]);
 
@@ -364,6 +408,19 @@ SparseVanka<number>::memory_consumption () const
 }
 
 
+template <typename number>
+SparseVanka<number>::AdditionalData::AdditionalData (
+  const std::vector<bool> &selected,
+  const bool               conserve_mem,
+  const unsigned int       n_threads)
+  :
+  selected(selected),
+  conserve_mem (conserve_mem),
+  n_threads (n_threads)
+{}
+
+
+//---------------------------------------------------------------------------
 
 
 template <typename number>
index 18a9ddf8134530ecb7cab71991506223f0625c2f..637eea45505721de345a543b37caea06e0b92bb8 100644 (file)
@@ -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)                        \
index eb3e2e6d13a1d364fcc06bb1e8cba0d459af0b0a..0fb5ca5d6d4fcc61796870774ccb50d2408b0cc7 100644 (file)
@@ -25,20 +25,29 @@ template void SparseMIC<double>::initialize<double> (const SparseMatrix<double>
                                                      const AdditionalData &data);
 template void SparseMIC<double>::vmult<double> (Vector<double> &,
                                                 const Vector<double> &) const;
+template void SparseMIC<double>::Tvmult<double> (Vector<double> &,
+                                                 const Vector<double> &) const;
 template void SparseMIC<double>::initialize<float> (const SparseMatrix<float> &,
                                                     const AdditionalData &data);
 template void SparseMIC<double>::vmult<float> (Vector<float> &,
                                                const Vector<float> &) const;
+template void SparseMIC<double>::Tvmult<float> (Vector<float> &,
+                                                const Vector<float> &) const;
 
 template class SparseMIC<float>;
 template void SparseMIC<float>::initialize<double> (const SparseMatrix<double> &,
                                                     const AdditionalData &data);
 template void SparseMIC<float>::vmult<double> (Vector<double> &,
                                                const Vector<double> &) const;
+template void SparseMIC<float>::Tvmult<double> (Vector<double> &,
+                                                const Vector<double> &) const;
 template void SparseMIC<float>::initialize<float> (const SparseMatrix<float> &,
                                                    const AdditionalData &data);
 template void SparseMIC<float>::vmult<float> (Vector<float> &,
                                               const Vector<float> &) const;
+template void SparseMIC<float>::Tvmult<float> (Vector<float> &,
+                                               const Vector<float> &) 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 (file)
index 0000000..965fda2
--- /dev/null
@@ -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 <deal.II/lac/block_sparse_matrix.h>
+#include <deal.II/lac/block_vector.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/linear_operator.h>
+#include <deal.II/lac/packaged_operation.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/precondition_block.h>
+#include <deal.II/lac/precondition_selector.h>
+#include <deal.II/lac/sparse_ilu.h>
+#include <deal.II/lac/sparse_decomposition.h>
+#include <deal.II/lac/sparse_mic.h>
+#include <deal.II/lac/sparse_vanka.h>
+
+#include <deal.II/lac/solver_control.h>
+#include <deal.II/lac/solver_bicgstab.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/solver_gmres.h>
+#include <deal.II/lac/solver_minres.h>
+#include <deal.II/lac/solver_qmrs.h>
+#include <deal.II/lac/solver_relaxation.h>
+#include <deal.II/lac/solver_richardson.h>
+#include <deal.II/lac/solver_selector.h>
+#include <deal.II/lac/sparse_direct.h>
+#include <deal.II/lac/iterative_inverse.h>
+
+
+using namespace dealii;
+
+
+template<class PRECONDITIONER, class MATRIX, class VECTOR,
+         class ADDITIONAL_DATA = typename PRECONDITIONER::AdditionalData>
+void
+test_preconditioner (const MATRIX &A,
+                     const VECTOR &b,
+                     const ADDITIONAL_DATA &data = ADDITIONAL_DATA())
+{
+  const auto lo_A = linear_operator<VECTOR>(A);
+
+  PRECONDITIONER preconditioner;
+  preconditioner.initialize(A, data);
+
+  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
+  {
+    // Using exemplar matrix
+    const auto lo_A_inv_approx = linear_operator<VECTOR>(A, preconditioner);
+    const VECTOR x_approx = lo_A_inv_approx*b;
+  }
+  {
+    // Stand-alone
+    const auto lo_A_inv_approx = linear_operator<VECTOR>(preconditioner);
+    const VECTOR x_approx = lo_A_inv_approx*b;
+  }
+}
+
+// For Vector <double>
+// Cannot use more generic function as Vector <double>
+// does not define vector_type
+template<class PRECONDITIONER>
+void
+test_preconditioner (const SparseMatrix<double> &A,
+                     const Vector<double> &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<double> > solver (solver_control);
+    const auto lo_A_inv = inverse_operator(lo_A,
+                                           solver,
+                                           preconditioner);
+
+    const Vector<double> 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<double> 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<double> x_approx = lo_A_inv_approx*b;
+    deallog.pop();
+  }
+}
+
+template<class SOLVER>
+void
+test_solver (const SparseMatrix<double> &A,
+             const Vector<double> &b)
+{
+  // Standard solver
+  {
+    deallog.push("Standard solver");
+    SolverControl solver_control (100, 1.0e-10);
+    SOLVER solver (solver_control);
+
+    PreconditionJacobi< SparseMatrix<double> > preconditioner;
+    preconditioner.initialize(A);
+
+    Vector<double> 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<double> > preconditioner;
+    preconditioner.initialize(A);
+
+    const auto lo_A_inv = inverse_operator(lo_A,
+                                           solver,
+                                           preconditioner);
+    const Vector<double> x = lo_A_inv*b;
+    deallog.pop();
+  }
+}
+
+template <typename MatrixType>
+class PreconditionBlockIdentity : public BlockMatrixBase<MatrixType>
+{
+public:
+
+  struct AdditionalData
+  {
+    AdditionalData () {}
+  };
+
+  virtual ~PreconditionBlockIdentity ()
+  {}
+
+  void
+  initialize (const BlockMatrixBase<MatrixType> &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<typename VECTOR>
+  void
+  vmult(VECTOR &dst, const VECTOR &src) const
+  {
+    dst = src;
+  }
+
+  template<class VECTOR>
+  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<double> A (sparsity_pattern);
+    Vector<double> 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<double> > PREC;
+      test_preconditioner<PREC>(A, b);
+    }
+    {
+      deallog << "PreconditionIdentity" << std::endl;
+      typedef PreconditionIdentity PREC;
+      test_preconditioner<PREC>(A, b);
+    }
+    {
+      deallog << "PreconditionJacobi" << std::endl;
+      typedef PreconditionJacobi< SparseMatrix<double> > PREC;
+      test_preconditioner<PREC>(A, b);
+    }
+    {
+      deallog << "PreconditionPSOR" << std::endl;
+      typedef PreconditionPSOR< SparseMatrix<double> > PREC;
+      std::vector<unsigned int> permutation(b.size());
+      std::vector<unsigned int> inverse_permutation(b.size());
+      test_preconditioner<PREC>(A, b,
+                                typename PREC::AdditionalData(permutation,
+                                                              inverse_permutation));
+    }
+    {
+      deallog << "PreconditionRichardson" << std::endl;
+      typedef PreconditionRichardson PREC;
+      test_preconditioner<PREC>(A, b);
+    }
+    {
+      deallog << "PreconditionSelector" << std::endl;
+      const auto lo_A = linear_operator(A);
+
+      PreconditionSelector< SparseMatrix<double>, Vector<double> >
+      preconditioner ("jacobi");
+      preconditioner.use_matrix(A);
+
+      SolverControl solver_control (100, 1.0e-10);
+      SolverCG< Vector<double> > solver (solver_control);
+
+      // Exact inverse
+      const auto lo_A_inv = inverse_operator(lo_A,
+                                             solver,
+                                             preconditioner);
+
+      const Vector<double> x = lo_A_inv*b;
+
+      // Approximate inverse
+      const auto lo_A_inv_approx = linear_operator(preconditioner);
+      const Vector<double> x_approx = lo_A_inv_approx*b;
+    }
+    {
+      deallog << "PreconditionSOR" << std::endl;
+      typedef PreconditionSOR< SparseMatrix<double> > PREC;
+      test_preconditioner<PREC>(A, b);
+    }
+    {
+      deallog << "PreconditionSSOR" << std::endl;
+      typedef PreconditionSSOR< SparseMatrix<double> > PREC;
+      test_preconditioner<PREC>(A, b);
+    }
+    {
+      deallog << "SparseILU" << std::endl;
+      typedef SparseILU<double> PREC;
+      test_preconditioner<PREC>(A, b);
+    }
+    {
+      deallog << "SparseMIC" << std::endl;
+      typedef SparseMIC<double> PREC;
+      test_preconditioner<PREC>(A, b);
+    }
+    {
+      deallog << "SparseVanka" << std::endl;
+      typedef SparseVanka<double> PREC;
+      typedef typename PREC::AdditionalData PREC_AD;
+      test_preconditioner<PREC>(A, b,
+                                PREC_AD(std::vector<bool>(rc, true)));
+    }
+    deallog.pop();
+
+    // === SOLVERS ===
+    deallog << std::endl;
+    deallog << "Solvers" << std::endl;
+    deallog.push("Solvers");
+
+    {
+      deallog << "SolverBicgstab" << std::endl;
+      typedef SolverBicgstab< Vector<double> > SLVR;
+      test_solver<SLVR> (A, b);
+    }
+    {
+      deallog << "SolverCG" << std::endl;
+      typedef SolverCG< Vector<double> > SLVR;
+      test_solver<SLVR> (A, b);
+    }
+    {
+      deallog << "SolverFGMRES" << std::endl;
+      typedef SolverFGMRES< Vector<double> > SLVR;
+      test_solver<SLVR> (A, b);
+    }
+    {
+      deallog << "SolverGMRES" << std::endl;
+      typedef SolverGMRES< Vector<double> > SLVR;
+      test_solver<SLVR> (A, b);
+    }
+    {
+      deallog << "SolverMinRes" << std::endl;
+      typedef SolverMinRes< Vector<double> > SLVR;
+      test_solver<SLVR> (A, b);
+    }
+    {
+      deallog << "SolverQMRS" << std::endl;
+      typedef SolverQMRS< Vector<double> > SLVR;
+      test_solver<SLVR> (A, b);
+    }
+    {
+      deallog << "SolverRelaxation" << std::endl;
+      typedef SolverRelaxation< Vector<double> > SLVR;
+      test_solver<SLVR> (A, b);
+    }
+    {
+      deallog << "SolverRichardson" << std::endl;
+      typedef SolverRichardson< Vector<double> > SLVR;
+      test_solver<SLVR> (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<double> > solver;
+      solver.select("cg");
+      solver.set_control(solver_control);
+
+      PreconditionJacobi< SparseMatrix<double> > preconditioner;
+      preconditioner.initialize(A);
+
+      const auto lo_A_inv = inverse_operator(lo_A,
+                                             solver,
+                                             preconditioner);
+
+      const Vector<double> b (rc);
+      const Vector<double> 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<double> x = lo_A_inv*b;
+    }
+    {
+      deallog << "IterativeInverse" << std::endl;
+
+      PreconditionJacobi< SparseMatrix<double> > preconditioner;
+      preconditioner.initialize(A);
+
+      ReductionControl solver_control (10, 1.e-30, 1.e-2);
+      IterativeInverse< Vector<double> > 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<double> 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<blks; ++bi)
+        for (unsigned int bj=0; bj<blks; ++bj)
+          csp.block(bi,bj).reinit(rc,rc);
+
+      csp.collect_sizes();
+      sparsity_pattern.copy_from(csp);
+    }
+
+    BlockSparseMatrix<double> A (sparsity_pattern);
+    BlockVector<double> b (blks,rc);
+    for (unsigned int bi=0; bi<blks; ++bi)
+      for (unsigned int i=0; i<rc; ++i)
+        {
+          A.block(bi,bi).diag_element(i) = 2.0;
+          b.block(bi)(i) = bi*rc + i;
+        }
+
+    // === PRECONDITIONERS ===
+    {
+      deallog << "PreconditionBlockIdentity" << std::endl;
+      typedef PreconditionBlockIdentity< SparseMatrix<double> > PREC;
+      test_preconditioner<PREC>(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 (file)
index 0000000..0aeb146
--- /dev/null
@@ -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

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.