--- /dev/null
+/*-------------------- precondition.h --------------------*/
+//$Id$
+// Guido Kanschat, University of Heidelberg, 1999
+
+#ifndef __lac_precondition_h
+#define __lac_precondition_h
+
+/**
+ * No preconditioning.
+ * @author Guido Kanschat, 1999
+ */
+template<class VECTOR>
+class PreconditionIdentity
+{
+ public:
+ /**
+ * Execute preconditioning.
+ */
+ void operator() (VECTOR&, const VECTOR&) const;
+};
+
+
+/**
+ * Preconditioner using a matrix-builtin function.
+ * This class forms a preconditioner suitable for the LAC solver
+ * classes. Since many preconditioning methods are based on matrix
+ * entries, these have to be implemented as member functions of the
+ * underlying matrix implementation. This class now is intended to
+ * allow easy access to these member functions from LAC solver
+ * classes.
+ *
+ * It seems that all builtin preconditioners have a relaxation
+ * parameter, so please use #PreconditionRelaxation# for these.
+ * @author Guido Kanschat, 1999
+ */
+
+template<class MATRIX, class VECTOR>
+class PreconditionUseMatrix
+{
+ public:
+ /**
+ * Type of the preconditioning
+ * function of the matrix.
+ */
+ typedef void ( MATRIX::* function_ptr)(VECTOR&, const VECTOR&) const;
+
+ /**
+ * Constructor.
+ * This constructor stores a
+ * reference to the matrix object
+ * for later use and selects a
+ * preconditioning method, which
+ * must be a member function of
+ * that matrix.
+ */
+ PreconditionUseMatrix(const MATRIX & M, function_ptr method);
+
+ /**
+ * Execute preconditioning.
+ */
+ void operator() (VECTOR&, const VECTOR&) const;
+
+ private:
+ /**
+ * Pointer to the matrix in use.
+ */
+ const MATRIX& matrix;
+ /**
+ * Pointer to the preconditioning
+ * function.
+ */
+ const function_ptr precondition;
+};
+
+/**
+ * Preconditioner for builtin relaxation methods.
+ * Uses methods of #SparseMatrix#.
+ * @author Guido Kanschat, 1999
+ */
+template<class MATRIX, class VECTOR>
+class PreconditionRelaxation
+{
+ public:
+ /**
+ * Type of the preconditioning
+ * function of the matrix.
+ */
+ typedef void ( MATRIX::* function_ptr)(VECTOR&, const VECTOR&,
+ typename MATRIX::entry_type) const;
+
+ /**
+ * Constructor.
+ * This constructor stores a
+ * reference to the matrix object
+ * for later use and selects a
+ * preconditioning method, which
+ * must be a member function of
+ * that matrix.
+ */
+ PreconditionRelaxation(const MATRIX & M, function_ptr method,
+ double omega = 1.);
+
+ /**
+ * Execute preconditioning.
+ */
+ void operator() (VECTOR&, const VECTOR&) const;
+
+ private:
+ /**
+ * Pointer to the matrix in use.
+ */
+ const MATRIX& matrix;
+ /**
+ * Pointer to the preconditioning
+ * function.
+ */
+ const function_ptr precondition;
+ /**
+ * Relaxation parameter.
+ */
+ double omega;
+};
+
+
+template<class VECTOR>
+void
+PreconditionIdentity<VECTOR>::operator() (VECTOR& dst, const VECTOR& src) const
+{
+ dst = src;
+}
+
+template<class MATRIX, class VECTOR>
+PreconditionUseMatrix<MATRIX, VECTOR>::PreconditionUseMatrix(const MATRIX& M,
+ function_ptr method)
+ :
+ matrix(M), precondition(method)
+{}
+
+
+template<class MATRIX, class VECTOR>
+void
+PreconditionUseMatrix<MATRIX, VECTOR>::operator() (VECTOR& dst,
+ const VECTOR& src) const
+{
+ (matrix.*precondition)(dst, src);
+}
+
+template<class MATRIX, class VECTOR>
+PreconditionRelaxation<MATRIX, VECTOR>::PreconditionRelaxation(const MATRIX& M,
+ function_ptr method,
+ double omega)
+ :
+ matrix(M), precondition(method), omega(omega)
+{}
+
+
+template<class MATRIX, class VECTOR>
+void
+PreconditionRelaxation<MATRIX, VECTOR>::operator() (VECTOR& dst,
+ const VECTOR& src) const
+{
+ (matrix.*precondition)(dst, src, omega);
+}
+
+#endif
/**
- * Base class for iterative solvers. This class defines possible
+ * Base class for iterative solvers.
+ *
+ * HAS TO BE UPDATED!
+ *
+ * This class defines possible
* return states of linear solvers and provides interfaces to a memory
* pool and the control object.
*
* object.
*/
Solver (SolverControl &, VectorMemory<Vector> &);
-
- /**
- * Destructor. This one is virtual,
- * overload it in derived classes if
- * necessary.
- */
- virtual ~Solver() {}
-
+
/**
- * Solver procedure.
+ * Solver procedure. This is in
+ * fact only a template for the
+ * solve-function to be
+ * implemented in derived classes
*/
- virtual ReturnState solve (const Matrix &A,
- Vector &x,
- const Vector &b) = 0;
+ template<class Preconditioner>
+ ReturnState solve (const Matrix &A,
+ Vector &x,
+ const Vector &b,
+ const Preconditioner& precondition) const;
/**
* Access to control object.
/**
* Solve primal problem only.
*/
- virtual ReturnState solve (const Matrix &A,
- Vector &x,
- const Vector &b);
+ template<class Preconditioner>
+ ReturnState solve (const Matrix &A,
+ Vector &x,
+ const Vector &b,
+ const Preconditioner& precondition);
protected:
/**
/**
* The iteration loop itself.
*/
- ReturnState iterate();
+ template<class Preconditioner>
+ ReturnState iterate(const Preconditioner& precondition);
};
{}
-template<class Matrix, class Vector> double
+template<class Matrix, class Vector>
+double
SolverBicgstab<Matrix, Vector>::criterion()
{
res = MA->residual(*Vt, *Vx, *Vb);
}
-template < class Matrix, class Vector > SolverControl::State
+template < class Matrix, class Vector >
+SolverControl::State
SolverBicgstab<Matrix, Vector>::start()
{
res = MA->residual(*Vr, *Vx, *Vb);
return state;
}
-template<class Matrix, class Vector> Solver<Matrix,Vector>::ReturnState
-SolverBicgstab<Matrix, Vector>::iterate()
+template<class Matrix, class Vector>
+template<class Preconditioner>
+Solver<Matrix,Vector>::ReturnState
+SolverBicgstab<Matrix, Vector>::iterate(const Preconditioner& precondition)
{
SolverControl::State state = SolverControl::iterate;
alpha = omega = rho = 1.;
beta = rhobar * alpha / (rho * omega);
rho = rhobar;
p.sadd(beta, 1., r, -beta*omega, v);
- MA->precondition(y,p);
+ precondition(y,p);
MA->vmult(v,y);
rhobar = rbar * v;
alpha = rho/rhobar;
s.equ(1., r, -alpha, v);
- MA->precondition(z,s);
+ precondition(z,s);
MA->vmult(t,z);
rhobar = t*s;
omega = rhobar/(t*t);
}
-template<class Matrix, class Vector> Solver<Matrix,Vector>::ReturnState
+template<class Matrix, class Vector>
+template<class Preconditioner>
+Solver<Matrix,Vector>::ReturnState
SolverBicgstab<Matrix, Vector>::solve(const Matrix &A,
Vector &x,
- const Vector &b)
+ const Vector &b,
+ const Preconditioner& precondition)
{
deallog.push("Bicgstab");
Vr = memory.alloc(); Vr->reinit(x);
{
deallog << "Go!" << endl;
if (start() == SolverControl::success) break;
- state = iterate();
+ state = iterate(precondition);
}
while (state == breakdown);
* @author Original implementation by G. Kanschat, R. Becker and F.-T. Suttmeier, reworking and documentation by Wolfgang Bangerth
*/
template<class Matrix, class Vector>
-class SolverPCG : public Solver<Matrix,Vector> {
+class SolverPCG : public Solver<Matrix,Vector>
+{
public:
/**
/**
* Solver method.
*/
- virtual ReturnState solve (const Matrix &A,
- Vector &x,
- const Vector &b);
+ template<class Preconditioner>
+ ReturnState solve (const Matrix &A,
+ Vector &x,
+ const Vector &b,
+ const Preconditioner& precondition);
protected:
/**
template<class Matrix, class Vector>
+template<class Preconditioner>
Solver<Matrix,Vector>::ReturnState
SolverPCG<Matrix,Vector>::solve (const Matrix &A,
Vector &x,
- const Vector &b) {
+ const Vector &b,
+ const Preconditioner& precondition)
+{
SolverControl::State conv=SolverControl::iterate;
// Memory allocation
};
g.scale(-1.);
- A.precondition(h,g);
+ precondition(h,g);
d.equ(-1.,h);
conv = control().check(it,res);
if (conv) break;
- A.precondition(h,g);
+ precondition(h,g);
beta = gh;
gh = g*h;
+++ /dev/null
-/*---------------------------- solver_gauss_seidel.h ---------------------------*/
-/* $Id$ */
-#ifndef __solver_gauss_seidel_H
-#define __solver_gauss_seidel_H
-/*---------------------------- solver_gauss_seidel.h ---------------------------*/
-
-#include <base/trace.h>
-#include <lac/solver_control.h>
-#include <lac/solver.h>
-
-/**
- * Implementation of the Gauss-Seidel method. The stopping criterion
- * is the norm of the defect, i.e. if x is the iteration vector, b the
- * rhs and A the matrix, the $res = \| A x - b \|$.
- */
-template<class Matrix, class Vector>
-class SolverGaussSeidel : public Solver<Matrix, Vector> {
- public:
- /**
- * Constructor. Takes a SolverControl and some
- VectorMemory. VectorMemory is not used.
- */
- SolverGaussSeidel (SolverControl &cn, VectorMemory<Vector> &mem) :
- Solver<Matrix,Vector> (cn,mem)
- {};
-
- /**
- * Solver method. Just specify the vectors A, x and b and let it
- * run. x should contain a reasonable starting value.
- */
- virtual ReturnState solve (const Matrix &A,
- Vector &x,
- const Vector &b);
-
- protected:
- /**
- * Implementation of the computation of
- * the norm of the residual.
- */
- virtual double criterion();
-
- /**
- * Within the iteration loop, the
- * square of the residual vector is
- * stored in this variable. The
- * function #criterion# uses this
- * variable to compute the convergence
- * value, which in this class is the
- * norm of the residual vector and thus
- * the square root of the #res2# value.
- */
- double res2;
-};
-
-/*----------------- Implementation of the Gauss-Seidel Method ------------------------*/
-
-template<class Matrix,class Vector>
-SolverGaussSeidel<Matrix,Vector>::ReturnState
-SolverGaussSeidel<Matrix,Vector>::solve (const Matrix &A,
- Vector &x,
- const Vector &b) {
-
- deallog.push("lac/solver_gauss_seidel.h");
- deallog.push("SolverGaussSeidel::solve");
-
- unsigned int i;
- unsigned int n = b.size();
- double res;
-
- SolverControl::State conv=SolverControl::iterate;
-
- Vector defect(n);
-
- deallog.push("Main loop");
-
- // Main loop
- for(int iter=0; conv==SolverControl::iterate; iter++)
- {
- // Calculate defect, i.e. defect = A x - b
- A.vmult(defect,x);
- defect.add(-1,b);
-
- // Calculate residual
- res2 = defect * defect;
-
- // Apply Gauss-Seidel preconditioner
- for (i=0;i<n;i++)
- {
- defect(i) = defect(i) / A(i,i);
- }
-
- // Correct defect
- x.add(-1,defect);
-
- // Check residual
- res = criterion();
- conv = (control().check(iter, res));
- if (conv != SolverControl::iterate)
- break;
- }
-
- // Output
-
- deallog.pop();
-
- if (conv == SolverControl::failure)
- {
- TRACEMSG("*** exceeded maximum number of iterations");
- deallog.pop();
- deallog.pop();
- return exceeded;
- }
- else
- {
- TRACEMSG("success");
- deallog.pop();
- deallog.pop();
- return success;
- }
-}
-
-template<class Matrix,class Vector>
-inline double
-SolverGaussSeidel<Matrix,Vector>::criterion()
-{
- return sqrt(res2);
-}
-
-/*---------------------------- solver_gauss_seidel.h ---------------------------*/
-/* end of #ifndef __solver_gauss_seidel_H */
-#endif
-/*---------------------------- solver_gauss_seidel.h ---------------------------*/
* computational speed against memory. One of the few other
* possibilities is to use a good preconditioner.
*
- * @author Original implementation by the DEAL authors; adapted, cleaned and documented by Wolfgang Bangerth */
+ * @author Wolfgang Bangerth
+ */
template<class Matrix, class Vector>
-class SolverGMRES : public Solver<Matrix, Vector> {
+class SolverGMRES : public Solver<Matrix, Vector>
+{
public:
/**
* Constructor.
/**
* Solver method.
*/
- virtual ReturnState solve (const Matrix &A,
- Vector &x,
- const Vector &b);
+ template<class Preconditioner>
+ ReturnState solve (const Matrix &A,
+ Vector &x,
+ const Vector &b,
+ const Preconditioner& precondition);
DeclException1 (ExcTooFewTmpVectors,
int,
template<class Matrix, class Vector>
+template<class Preconditioner>
Solver<Matrix,Vector>::ReturnState
SolverGMRES<Matrix,Vector>::solve (const Matrix& A,
Vector & x,
- const Vector& b)
+ const Vector& b,
+ const Preconditioner& precondition)
{
// this code was written by the fathers of
// DEAL. I take absolutely no guarantees
if (left_precondition)
{
A.residual(p,x,b);
- A.precondition(v,p);
+ precondition(v,p);
} else {
A.residual(v,x,b);
};
if (left_precondition)
{
A.vmult(p, *tmp_vectors[inner_iteration]);
- A.precondition(vv,p);
+ precondition(vv,p);
} else {
- A.precondition(p,*tmp_vectors[inner_iteration]);
+ precondition(p,*tmp_vectors[inner_iteration]);
A.vmult(vv,p);
};
p = 0.;
for (unsigned int i=0; i<dim; ++i)
p.add(h(i), *tmp_vectors[i]);
- A.precondition(v,p);
+ precondition(v,p);
x.add(1.,v);
};
/**
* Solve $Ax=b$ for $x$.
*/
- virtual ReturnState solve (const Matrix &A,
- Vector &x,
- const Vector &b);
+ template<class Preconditioner>
+ ReturnState solve (const Matrix &A,
+ Vector &x,
+ const Vector &b,
+ const Preconditioner& precondition);
/**
* Set the damping-coefficient.
template<class Matrix,class Vector>
+template<class Preconditioner>
Solver<Matrix,Vector>::ReturnState
SolverRichardson<Matrix,Vector>::solve (const Matrix &A,
- Vector &x,
- const Vector &b) {
+ Vector &x,
+ const Vector &b,
+ const Preconditioner& precondition)
+{
SolverControl::State conv=SolverControl::iterate;
// Memory allocation
if (conv != SolverControl::iterate)
break;
- A.precondition(d,r);
+ precondition(d,r);
x.add(omega,d);
}