/**
- * Base class for iterative solvers.
- *
-//TODO:[?] * HAS TO BE UPDATED!
- *
- * This class defines possible
+ * Base class for iterative solvers. This class defines possible
* return states of linear solvers and provides interfaces to a memory
* pool and the control object.
*
*
* @sect3{Requirements for template classes}
*
- * The class is templated to allow for different matrix and vector
- * classes, since iterative solvers do not rely on any special structure
- * of matrices or the format of storage. However, there are some common
- * requirements a matrix or vector type must fulfil to qualify as an
- * applicable type for the solvers in this hierarchy. These requirements
- * are listed following. The classes do not declare any concrete
- * class, they are rather intended to form a `signature' which a concrete
- * class has to conform to.
+ * Since iterative solvers do not rely on any special structure of
+ * matrices or the format of storage, but only require that matrices
+ * and vector define certain operations such as matrix-vector
+ * products, or scalar products between vectors, this class as well as
+ * the derived classes implementing concrete linear solvers are
+ * templated on the types of matrices and vectors. However, there are
+ * some common requirements a matrix or vector type must fulfill to
+ * qualify as an applicable type for the solvers in this
+ * hierarchy. These requirements are listed following. The listed
+ * classes are not any concrete class, they are rather intended to
+ * form a `signature' which a concrete class has to conform to. Note
+ * that the matrix and vector classes within this library of course
+ * conform to this interface.
*
* @begin{verbatim}
* class Matrix
* {
* public:
- * // Application to a Vector
- * void vmult (Vector& dst, const Vector& src) const;
- *
- * // Application of a preconditioner to
- * // a Vector, i.e. $dst=\tilde A^(-1) src$,
- * // where $\tilde A^(-1)$ is an approximation
- * // to the inverse if the matrix stored in
- * // this object.
- * void precondition (Vector& dst, const Vector& src) const;
+ * // Application of matrix to vector src.
+ * // write result into dst
+ * void vmult (Vector &dst, const Vector &src) const;
*
* // Application of transpose to a Vector.
- * // Only used by special iterative methods.
- * void T_vmult (Vector& dst, const Vector& src) const;
- *
- * // Application of a transposed preconditioner
- * // to a Vector. Only used by special
- * // iterative methods
- *
- * void T_precondition (Vector& dst, const Vector& src) const;
+ * // Only used by certain iterative methods.
+ * void Tvmult (Vector &dst, const Vector &src) const;
* };
*
*
* class Vector
* {
* public:
+ * // resize and/or clear vector. note
+ * // that the second argument must have
+ * // a default value equal to false
+ * void reinit (const unsigned int size,
+ * bool leave_elements_uninitialized = false);
+ *
* // scalar product
- * double operator * (const Vector& v) const;
+ * double operator * (const Vector &v) const;
*
* // addition of vectors
* // $y = y + x$.
- * void add (const Vector& x);
+ * void add (const Vector &x);
+ *
* // $y = y + ax$.
- * void add (double a, const Vector& x);
+ * void add (const double a,
+ * const Vector &x);
*
- * // scaled addition of vectors
- * // $y = ay + x$.
- * void sadd (double a,
- * const Vector& x);
* // $y = ay + bx$.
- * void sadd (double a,
- * double b, const Vector& x);
- * // $y = ay + bx + cz$.
- * void sadd (double a,
- * double b, const Vector& x,
- * double c, const Vector& z);
+ * void sadd (const double a,
+ * const double b,
+ * const Vector &x);
*
* // $y = ax$.
- * void equ (double a, const Vector& x);
- * // $y = ax + bz$.
- * void equ (double a, const Vector& x,
- * double b, const Vector& z);
+ * void equ (const double a,
+ * const Vector &x);
+ *
+ * // scale the elements of the vector
+ * // by a fixed value
+ * void scale (const double a);
+ *
+ * // return the l2 norm of the vector
+ * double l2_norm () const;
* };
* @end{verbatim}
*
+ * In addition, for some solvers there has to be a global function
+ * @p{swap(vector &a, vector &b)} that exchanges the values of the two vectors.
+ *
+ * The preconditioners used must have the same interface as matrices,
+ * i.e. in particular they have to provide a member function @p{vmult}
+ * which denotes the application of the preconditioner.
+ *
*
* @sect3{AdditionalData}
*
/**
* Bicgstab algorithm by van der Vorst.
*
+ * For the requirements on matrices and vectors in order to work with
+ * this class, see the documentation of the @ref{Solver} base class.
+ *
* Like all other solver classes, this class has a local structure called
* @p{AdditionalData} which is used to pass additional parameters to the
* solver, like damping parameters or the number of temporary vectors. We
/**
* Preconditioned cg method.
*
+ * For the requirements on matrices and vectors in order to work with
+ * this class, see the documentation of the @ref{Solver} base class.
+ *
* Like all other solver classes, this class has a local structure called
* @p{AdditionalData} which is used to pass additional parameters to the
* solver, like damping parameters or the number of temporary vectors. We
virtual ~SolverCG ();
/**
- * Solver method.
+ * Solve the linear system $Ax=b$
+ * for x.
*/
template<class MATRIX, class PRECONDITIONER>
void
- solve (const MATRIX &A,
- VECTOR &x,
- const VECTOR &b,
- const PRECONDITIONER& precondition);
+ solve (const MATRIX &A,
+ VECTOR &x,
+ const VECTOR &b,
+ const PRECONDITIONER &precondition);
protected:
/**
template<class VECTOR>
template<class MATRIX, class PRECONDITIONER>
void
-SolverCG<VECTOR>::solve (const MATRIX &A,
- VECTOR &x,
- const VECTOR &b,
- const PRECONDITIONER& precondition)
+SolverCG<VECTOR>::solve (const MATRIX &A,
+ VECTOR &x,
+ const VECTOR &b,
+ const PRECONDITIONER &precondition)
{
SolverControl::State conv=SolverControl::iterate;
* of temporary vectors as commented upon above. By default, the number
* of these vectors is set to 30.
*
+ * For the requirements on matrices and vectors in order to work with
+ * this class, see the documentation of the @ref{Solver} base class.
+ *
* @author Wolfgang Bangerth
*/
template <class VECTOR = Vector<double> >
const AdditionalData &data=AdditionalData());
/**
- * Solver method.
+ * Solve the linear system $Ax=b$
+ * for x.
*/
template<class MATRIX, class PRECONDITIONER>
- void solve (const MATRIX &A,
- VECTOR &x,
- const VECTOR &b,
- const PRECONDITIONER& precondition);
+ void
+ solve (const MATRIX &A,
+ VECTOR &x,
+ const VECTOR &b,
+ const PRECONDITIONER &precondition);
DeclException1 (ExcTooFewTmpVectors,
int,
template<class VECTOR>
template<class MATRIX, class PRECONDITIONER>
void
-SolverGMRES<VECTOR>::solve (const MATRIX& A,
- VECTOR & x,
- const VECTOR& b,
- const PRECONDITIONER& precondition)
+SolverGMRES<VECTOR>::solve (const MATRIX &A,
+ VECTOR &x,
+ const VECTOR &b,
+ const PRECONDITIONER &precondition)
{
// this code was written a very
// long time ago by people not
/**
* Preconditioned MinRes method.
*
+ * For the requirements on matrices and vectors in order to work with
+ * this class, see the documentation of the @ref{Solver} base class.
+ *
* Like all other solver classes, this class has a local structure called
* @p{AdditionalData} which is used to pass additional parameters to the
* solver, like damping parameters or the number of temporary vectors. We
virtual ~SolverMinRes ();
/**
- * Solver method.
+ * Solve the linear system $Ax=b$
+ * for x.
*/
template<class MATRIX, class PRECONDITIONER>
void
- solve (const MATRIX &A,
- VECTOR &x,
- const VECTOR &b,
- const PRECONDITIONER& precondition);
+ solve (const MATRIX &A,
+ VECTOR &x,
+ const VECTOR &b,
+ const PRECONDITIONER &precondition);
/**
* Exception
template<class VECTOR>
template<class MATRIX, class PRECONDITIONER>
typename Solver<VECTOR>::ReturnState
-SolverMinRes<VECTOR>::solve (const MATRIX &A,
- VECTOR &x,
- const VECTOR &b,
- const PRECONDITIONER& precondition)
+SolverMinRes<VECTOR>::solve (const MATRIX &A,
+ VECTOR &x,
+ const VECTOR &b,
+ const PRECONDITIONER &precondition)
{
SolverControl::State conv=SolverControl::iterate;
r_l2 = r0;
- u[0].reinit(VS,0);
+ u[0].reinit(VS);
delta[0] = 1.;
- m[0].reinit(VS,0);
- m[1].reinit(VS,0);
- m[2].reinit(VS,0);
+ m[0].reinit(VS);
+ m[1].reinit(VS);
+ m[2].reinit(VS);
conv = control().check(0,r_l2);
while (conv==SolverControl::iterate)
- {
-
+ {
if (delta[1]!=0)
v.scale(1./sqrt(delta[1]));
else
- v.reinit(VS,0);
+ v.reinit(VS);
A.vmult(u[2],v);
u[2].add (-sqrt(delta[1]/delta[0]), u[0]);
* preconditioner is used: left preconditioning seems to require the
* inverse.
*
+ * For the requirements on matrices and vectors in order to work with
+ * this class, see the documentation of the @ref{Solver} base class.
+ *
* Like all other solver classes, this class has a local structure called
* @p{AdditionalData} which is used to pass additional parameters to the
* solver, like damping parameters or the number of temporary vectors. We
const AdditionalData &data=AdditionalData());
/**
- * Solver method.
+ * Solve the linear system $Ax=b$
+ * for x.
*/
template<class MATRIX, class PRECONDITIONER>
void
- solve (const MATRIX &A,
- VECTOR &x,
- const VECTOR &b,
- const PRECONDITIONER& precondition);
+ solve (const MATRIX &A,
+ VECTOR &x,
+ const VECTOR &b,
+ const PRECONDITIONER &precondition);
/**
* Interface for derived class.
template<class VECTOR>
template<class MATRIX, class PRECONDITIONER>
void
-SolverQMRS<VECTOR>::solve (const MATRIX &A,
- VECTOR &x,
- const VECTOR &b,
- const PRECONDITIONER& precondition)
+SolverQMRS<VECTOR>::solve (const MATRIX &A,
+ VECTOR &x,
+ const VECTOR &b,
+ const PRECONDITIONER &precondition)
{
deallog.push("QMRS");
template<class VECTOR>
template<class MATRIX, class PRECONDITIONER>
bool
-SolverQMRS<VECTOR>::iterate(const MATRIX& A,
- const PRECONDITIONER& precondition)
+SolverQMRS<VECTOR>::iterate(const MATRIX &A,
+ const PRECONDITIONER &precondition)
{
/* Remark: the matrix A in the article is the preconditioned matrix.
* Therefore, we have to precondition x before we compute the first residual.
precondition.vmult(q,p);
tau = v.norm_sqr();
- //deallog << "tau:" << tau << std::endl;
rho = q*v;
- //deallog << "rho:" << rho << std::endl;
-
-
-while (state == SolverControl::iterate)
+
+ while (state == SolverControl::iterate)
{
step++; it++;
// Step 1
return true;
// Step 3
alpha = rho/sigma;
- //deallog << "alpha:" << alpha << std::endl;
v.add(-alpha,t);
// Step 4
psi = 1./(1.+theta);
tau *= theta*psi;
- //deallog << "psi:" << psi << std::endl;
- //deallog << "theta:" << theta << std::endl;
- //deallog << "tau:" << tau << std::endl;
-
d.sadd(psi*theta_old, psi*alpha, p);
x.add(d);
* Implementation of the richardson iteration method. The stopping criterion
* is the norm of the residual.
*
+ * For the requirements on matrices and vectors in order to work with
+ * this class, see the documentation of the @ref{Solver} base class.
+ *
* Like all other solver classes, this class has a local structure called
* @p{AdditionalData} which is used to pass additional parameters to the
* solver, like damping parameters or the number of temporary vectors. We
virtual ~SolverRichardson ();
/**
- * Solve $Ax=b$ for $x$.
+ * Solve the linear system $Ax=b$
+ * for x.
*/
template<class MATRIX, class PRECONDITIONER>
- void solve (const MATRIX &A,
- VECTOR &x,
- const VECTOR &b,
- const PRECONDITIONER& precondition);
+ void
+ solve (const MATRIX &A,
+ VECTOR &x,
+ const VECTOR &b,
+ const PRECONDITIONER &precondition);
/**
* Solve $A^Tx=b$ for $x$.
*/
template<class MATRIX, class PRECONDITIONER>
- void Tsolve (const MATRIX &A,
- VECTOR &x,
- const VECTOR &b,
- const PRECONDITIONER& precondition);
+ void
+ Tsolve (const MATRIX &A,
+ VECTOR &x,
+ const VECTOR &b,
+ const PRECONDITIONER &precondition);
/**
* Set the damping-coefficient.
template<class VECTOR>
template<class MATRIX, class PRECONDITIONER>
void
-SolverRichardson<VECTOR>::solve (const MATRIX &A,
- VECTOR &x,
- const VECTOR &b,
- const PRECONDITIONER& precondition)
+SolverRichardson<VECTOR>::solve (const MATRIX &A,
+ VECTOR &x,
+ const VECTOR &b,
+ const PRECONDITIONER &precondition)
{
SolverControl::State conv=SolverControl::iterate;
template<class VECTOR>
template<class MATRIX, class PRECONDITIONER>
void
-SolverRichardson<VECTOR>::Tsolve (const MATRIX &A,
- VECTOR &x,
- const VECTOR &b,
- const PRECONDITIONER& precondition)
+SolverRichardson<VECTOR>::Tsolve (const MATRIX &A,
+ VECTOR &x,
+ const VECTOR &b,
+ const PRECONDITIONER &precondition)
{
SolverControl::State conv=SolverControl::iterate;
projects / dealii-svn.git / commitdiff