(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>
(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.
* 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().
*/
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
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.
}
+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
}
+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
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);
};
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);
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);
};
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);
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
*/
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;
};
* 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
PointerMatrix(const MATRIX *M,
const char *name);
+ /**
+ * Destructor
+ */
+ virtual ~PointerMatrix ()
+ {}
+
// Use doc from base class
virtual void clear();
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;
};
* 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
const MATRIX *M,
const char *name);
+ /**
+ * Destructor
+ */
+ virtual ~PointerMatrixAux ()
+ {}
+
// Use doc from base class
virtual void clear();
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.
/**
* The pointer to the actual matrix.
*/
- SmartPointer<const MATRIX,PointerMatrixAux<MATRIX,VECTOR> > m;
+ SmartPointer<const MATRIX,PointerMatrixAux<MATRIX,VECTOR> > mtrx;
};
* (#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
PointerMatrixVector (const Vector<number> *M,
const char *name);
+ /**
+ * Destructor
+ */
+ virtual ~PointerMatrixVector ()
+ {}
+
// Use doc from base class
virtual void clear();
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;
};
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)
{}
PointerMatrix<MATRIX, VECTOR>::PointerMatrix (
const MATRIX *M,
const char *name)
- : m(M, name)
+ : mtrx(M, name)
{}
inline void
PointerMatrix<MATRIX, VECTOR>::clear ()
{
- m = 0;
+ mtrx = 0;
}
inline const PointerMatrix<MATRIX, VECTOR> &
PointerMatrix<MATRIX, VECTOR>::operator= (const MATRIX *M)
{
- m = M;
+ mtrx = M;
return *this;
}
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>
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);
}
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);
}
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);
}
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();
+}
+
//----------------------------------------------------------------------//
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;
}
VectorMemory<VECTOR> *mem,
const char *name)
: mem(mem, name),
- m(0, name)
+ mtrx(0, name)
{
if (mem == 0) mem = &my_memory;
}
const MATRIX *M,
const char *name)
: mem(mem, name),
- m(M, name)
+ mtrx(M, name)
{
if (mem == 0) mem = &my_memory;
}
inline void
PointerMatrixAux<MATRIX, VECTOR>::clear ()
{
- m = 0;
+ mtrx = 0;
}
inline const PointerMatrixAux<MATRIX, VECTOR> &
PointerMatrixAux<MATRIX, VECTOR>::operator= (const MATRIX *M)
{
- m = M;
+ mtrx = M;
return *this;
}
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>
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);
}
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);
}
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);
}
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)
{}
PointerMatrixVector<number>::PointerMatrixVector (
const Vector<number> *M,
const char *name)
- : m(M, name)
+ : mtrx(M, name)
{}
inline void
PointerMatrixVector<number>::clear ()
{
- m = 0;
+ mtrx = 0;
}
inline const PointerMatrixVector<number> &
PointerMatrixVector<number>::operator= (const Vector<number> *M)
{
- m = M;
+ mtrx = M;
return *this;
}
inline bool
PointerMatrixVector<number>::empty () const
{
- if (m == 0)
+ if (mtrx == 0)
return true;
- return m->empty();
+ return mtrx->empty();
}
template<typename number>
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;
}
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);
}
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;
}
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();
}
* 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.
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.
* 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;
};
* 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.
*/
};
/**
- * Constructor, sets the relaxation parameter to its default.
+ * Constructor, sets the relaxation parameter, domain and range sizes
+ * to their default.
*/
PreconditionRichardson();
* compatibility with more complex preconditioners.
*/
template <class MATRIX>
- void initialize (const MATRIX &,
+ void initialize (const MATRIX &matrix,
const AdditionalData ¶meters);
/**
*/
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;
};
* 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.
*/
*/
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.
/**
* 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.
* 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>
/**
* 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
+ ¶meters = 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
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.
*/
* 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
/**
* 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
*/
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:
/**
#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>
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
inline
PreconditionRichardson::PreconditionRichardson ()
:
- relaxation(0)
+ relaxation(0),
+ _m (0),
+ _n (0)
{
AdditionalData add_data;
relaxation=add_data.relaxation;
template <class MATRIX>
inline void
PreconditionRichardson::initialize (
- const MATRIX &,
+ const MATRIX &matrix,
const PreconditionRichardson::AdditionalData ¶meters)
{
relaxation = parameters.relaxation;
+ _m = matrix.m();
+ _n = matrix.n();
}
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>
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();
+}
//---------------------------------------------------------------------------
}
+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
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 ¶meters)
+ :
+ permutation(permutation),
+ inverse_permutation(inverse_permutation),
+ parameters(parameters)
+{
+
+}
+
//---------------------------------------------------------------------------
}
+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
* 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
*/
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.
*/
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
}
+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()
{
* 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>,
/**
* 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
*/
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.
}
+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;
+}
//---------------------------------------------------------------------------
*
* @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
{
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;
+
/**
* @}
*/
"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
/**
* Declare type for container size.
*/
- typedef types::global_dof_index size_type;
+ typedef typename SparseLUDecomposition<number>::size_type size_type;
/**
* Constructor. Does nothing.
* 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>
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.
}
+// 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
* @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
*/
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
*/
~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.
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
/**
* 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
*/
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.
*/
};
/*@}*/
+/* ---------------------------------- 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
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 ();
}
+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,
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)
{
// 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);
}
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
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
}
-
template<typename number>
template<typename number2>
void
// 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);
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]);
}
+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>
SparseDirectUMFPACK::SparseDirectUMFPACK ()
:
+ _m (0),
+ _n (0),
symbolic_decomposition (0),
numeric_decomposition (0),
control (UMFPACK_CONTROL)
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
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) \
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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+ }
+
+}
--- /dev/null
+
+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