--- /dev/null
+/*---------------------------- solver.h ---------------------------*/
+/* $Id$ */
+#ifndef __solver_H
+#define __solver_H
+/*---------------------------- solver.h ---------------------------*/
+
+
+
+// forward declaration
+class SolverControl;
+template <class Vector> class VectorMemory;
+
+
+
+/**
+ * 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.
+ *
+ * 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.
+ * \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 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;
+ * };
+ *
+ *
+ * class Vector
+ * {
+ * public:
+ * // scalar product
+ * double operator * (const Vector& v) const;
+ *
+ * // addition of vectors
+ * // $y = y + x$.
+ * void add (const Vector& x);
+ * // $y = y + ax$.
+ * void add (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);
+ *
+ * // $y = ax$.
+ * void equ (double a, const Vector& x);
+ * // $y = ax + bz$.
+ * void equ (double a, const Vector& x,
+ * double b, const Vector& z);
+ * };
+ * \end{verbatim}
+ */
+template <class Matrix, class Vector>
+class Solver
+{
+ public:
+ /**
+ * Declare possible return values of a
+ * solver object.
+ */
+ enum ReturnState {
+ success=0, exceeded, breakdown
+ };
+
+
+
+ /**
+ * Constructor. Assign a control
+ * object which stores the required
+ * precision and an object to provide
+ * memory.
+ *
+ * Of both objects, a reference is
+ * stored, so it is the user's
+ * responsibility to guarantee that the
+ * lifetime of the two arguments is at
+ * least as long as that of the solver
+ * object.
+ */
+ Solver (SolverControl &, VectorMemory<Vector> &);
+
+ /**
+ * Destructor. This one is virtual,
+ * overload it in derived classes if
+ * necessary.
+ */
+ virtual ~Solver() {}
+
+ /**
+ * Solver procedure.
+ */
+ virtual ReturnState solve (const Matrix &A,
+ Vector &x,
+ const Vector &b) = 0;
+
+ /**
+ * Access to control object.
+ */
+ SolverControl& control() const;
+
+ protected:
+ /**
+ * Calculation of convergence
+ * criterion. To allow further
+ * flexibility in solvers, the
+ * convergence criterion can be
+ * implemented by the user. Each
+ * method has it's standard
+ * criterion (usually some kind of
+ * residual) already implemented.
+ */
+ virtual double criterion() = 0;
+
+ /**
+ * Additional entry point for
+ * examination. Allows access to
+ * internal variables of a solver
+ * at a given point dependent of
+ * the actual solver method. Does
+ * nothing by default and should
+ * only be implemented to test
+ * iterative methods.
+ */
+ virtual void expert() {};
+
+ /**
+ * Control structure.
+ */
+ SolverControl &cntrl;
+
+ /**
+ * Memory for auxilliary vectors.
+ */
+ VectorMemory<Vector> &memory;
+};
+
+
+
+
+
+/**
+ * Base class for non-symmetric linear solvers. A second entry point
+ * allows the simultaneous computation of a dual problem.
+ */
+template<class Matrix, class Vector>
+class SolverDual : public Solver<Matrix, Vector>
+{
+ public:
+ /**
+ * Constructor.
+ */
+ SolverDual(SolverControl&, VectorMemory<Vector>&);
+
+ /**
+ * Solve the original problem
+ * $Ax=b$.
+ */
+ virtual ReturnState solve (const Matrix &A,
+ Vector &x,
+ const Vector &b) = 0;
+
+ /**
+ * Solve the two problems
+ * $Ax=b1$ and $A^Tz=b2$ simultanously.
+ */
+ virtual ReturnState solve (const Matrix &A,
+ Vector &x,
+ const Vector &b1,
+ Vector &z,
+ const Vector &b2) = 0;
+};
+
+
+
+
+
+
+/*-------------------------------- Inline functions ------------------------*/
+
+template <class Matrix, class Vector>
+inline
+SolverControl & Solver<Matrix,Vector>::control() const {
+ return cntrl;
+};
+
+
+
+template<class Matrix, class Vector>
+inline
+Solver<Matrix, Vector>::Solver(SolverControl &cn, VectorMemory<Vector> &mem)
+ : cntrl(cn),
+ memory(mem)
+{};
+
+
+
+
+
+/*---------------------------- solver.h ---------------------------*/
+/* end of #ifndef __solver_H */
+#endif
+/*---------------------------- solver.h ---------------------------*/
--- /dev/null
+/*---------------------------- solver_pcg.h ---------------------------*/
+/* $Id$ */
+#ifndef __solver_pcg_H
+#define __solver_pcg_H
+/*---------------------------- solver_pcg.h ---------------------------*/
+
+
+
+#include <lac/solver.h>
+#include <lac/solver_control.h>
+
+
+
+
+/**
+ * Preconditioned cg method.
+ */
+template<class Matrix, class Vector>
+class SolverPCG : public Solver<Matrix,Vector> {
+ public:
+
+ /**
+ * Constructor.
+ */
+ SolverPCG (SolverControl &cn, VectorMemory<Vector> &mem) :
+ Solver<Matrix,Vector>(cn,mem) {};
+
+ /**
+ * Solver method.
+ */
+ 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();
+
+ /**
+ * Temporary vectors, allocated through
+ * the #VectorMemory# object at the start
+ * of the actual solution process and
+ * deallocated at the end.
+ */
+ Vector *Vr, *Vp, *Vz, *VAp;
+
+ /**
+ * 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 ----------------------------*/
+
+
+template<class Matrix, class Vector>
+double SolverPCG<Matrix,Vector>::criterion()
+{
+ return sqrt(res2);
+};
+
+
+
+template<class Matrix, class Vector>
+Solver<Matrix,Vector>::ReturnState
+SolverPCG<Matrix,Vector>::solve (const Matrix &A,
+ Vector &x,
+ const Vector &b) {
+ SolverControl::State conv=SolverControl::iterate;
+
+ // Memory allocation
+ Vr = memory.alloc(); Vector& g = *Vr;
+ Vp = memory.alloc(); Vector& h = *Vp;
+ Vz = memory.alloc(); Vector& d = *Vz;
+ VAp = memory.alloc(); Vector& Ad = *VAp;
+
+ // Implementation taken from the DEAL
+ // library
+ int it=0;
+ double res,gh,alpha,beta;
+
+ res = A.residual(g,x,b);
+ conv = control().check(0,res);
+ if (conv)
+ {
+ memory.free(Vr);
+ memory.free(Vp);
+ memory.free(Vz);
+ memory.free(VAp);
+
+ return success;
+ };
+
+ g.scale(-1.);
+ A.precondition(h,g);
+
+ d.equ(-1.,h);
+
+ gh = g*h;
+
+ while (conv == SolverControl::iterate)
+ {
+ A.vmult(Ad,d);
+
+ alpha = d*Ad;
+ alpha = gh/alpha;
+
+ g.add(alpha,Ad);
+ x.add(alpha,d );
+ res = sqrt(g*g);
+
+ conv = control().check(it,res);
+ if (conv) break;
+
+ A.precondition(h,g);
+
+ beta = gh;
+ gh = g*h;
+ beta = gh/beta;
+
+ d.sadd(beta,-1.,h);
+ it++;
+ };
+
+
+ // Deallocate Memory
+
+ memory.free(Vr);
+ memory.free(Vp);
+ memory.free(Vz);
+ memory.free(VAp);
+
+ // Output
+
+ if (conv == SolverControl::failure)
+ return exceeded;
+ else
+ return success;
+};
+
+
+
+
+/*---------------------------- solver_pcg.h ---------------------------*/
+/* end of #ifndef __solver_pcg_H */
+#endif
+/*---------------------------- solver_pcg.h ---------------------------*/
--- /dev/null
+/*---------------------------- solver_control.h ---------------------------*/
+/* $Id$ */
+#ifndef __solver_control_H
+#define __solver_control_H
+/*---------------------------- solver_control.h ---------------------------*/
+
+
+
+/**
+ * Control class for iterative solvers.
+ *
+ * Used by iterative methods to
+ * determine whether the iteration should be continued. To this respect,
+ * the virtual function #check()# is called in each iteration
+ * with the current iteration
+ * step and the value indicating convergence (usually the residual).
+ *
+ * After the iteration has terminated, the functions #last_value# and
+ * #last_step# can be used to obtain information about the final state
+ * of the iteration.
+ *
+ * #check()# can be replaced in derived classes to allow for more
+ * sophisticated tests.
+ */
+class SolverControl {
+ public:
+ /**
+ * Return states of the check
+ * function, which indicate the state the
+ * solver is in.
+ *
+ * The possible values of State are
+ * <OL>
+ * <LI> #iterate = 0#: continue
+ * the iteration.
+ * <LI> #success#: the goal is reached,
+ * the iterative method can terminate
+ * successfully.
+ * <LI> #failure#!: the iterative
+ * method should stop because
+ * convergence cannot be achieved or at
+ * least was not achieved within the given
+ * maximal number of iterations.
+ * </OL>
+ */
+ enum State {
+ iterate = 0, success, failure
+ };
+
+ /**
+ * Constructor. The parameters
+ * #n# and #tol# are the
+ * maximum number of iteration
+ * steps before failure and the
+ * tolerance to determine success
+ * of the iteration. */
+ SolverControl (const unsigned int n, const double tol);
+
+ /**
+ * Decide about success or failure
+ * of an iteration. This function
+ * gets the current iteration step
+ * to determine, whether the
+ * allowed number of steps has
+ * been exceeded and returns
+ * #failure# in this case. If
+ * #check_value# is below the
+ * prescribed tolerance, it
+ * returns #success#. In all
+ * other cases #iterate# is
+ * returned to suggest
+ * continuation of the iterative
+ * procedure.
+ *
+ * #check()# additionally
+ * preserves #step# and
+ * #check_value#. These
+ * values are accessible by
+ * #last_value()# and
+ * #last_step()#.
+ *
+ * Derived classes may overload this
+ * function, e.g. to log the convergence
+ * indicators (#check_value#) or to do
+ * other computations.
+ */
+ virtual State check (const unsigned int step, const double check_value);
+
+ /**
+ * Return the convergence value of last
+ * iteration step for which #check# was
+ * called by the solver.
+ */
+ double last_value() const;
+
+ /**
+ * Number of last iteration step.
+ */
+ unsigned int last_step() const;
+
+ protected:
+ /**
+ * Maximum number of steps.
+ */
+ const unsigned int maxsteps;
+
+ /**
+ * Prescribed tolerance to be achieved.
+ */
+ const double tol;
+
+ /**
+ * Last value of the convergence criterion.
+ */
+ double lvalue;
+
+ /**
+ * Last step.
+ */
+ unsigned int lstep;
+};
+
+
+
+
+/**
+ * Specialization of #SolverControl# which returns #success# if either
+ * the specified tolerance is achieved or if the
+ * initial residual (or whatever criterion was chosen by the solver
+ * class) is reduced by a given factor. This is useful in cases where
+ * you don't want to solve exactly, but rather want to gain two digits.
+ */
+class ReductionControl : public SolverControl {
+ public:
+ /**
+ * Constructor. Provide the
+ * reduction factor additional to
+ * the arguments of the Control
+ * constructor.
+ */
+ ReductionControl (const unsigned int maxiter,
+ const double tolerance,
+ const double reduce);
+
+ /**
+ * Decide about success or failure
+ * of an iteration. This function
+ * calls the one in the base
+ * class, but sets the tolerance
+ * to #reduction * initial value#
+ * upon the first iteration.
+ */
+ virtual State check (const unsigned int step,
+ const double check_value);
+
+ /**
+ * Return the initial convergence
+ * criterion.
+ */
+ double initial_value() const;
+
+
+ protected:
+ /**
+ * Desired reduction factor.
+ */
+ const double reduce;
+
+ /**
+ * Initial value.
+ */
+ double initial_val;
+
+ /**
+ * Reduced tolerance. Stop iterations
+ * if either this value is achieved
+ * or if the base class indicates
+ * success.
+ */
+ double reduced_tol;
+};
+
+
+
+
+/*---------------------------- solver_control.h ---------------------------*/
+/* end of #ifndef __solver_control_H */
+#endif
+/*---------------------------- solver_control.h ---------------------------*/
--- /dev/null
+/*---------------------------- vector_memory.h ---------------------------*/
+/* $Id$ */
+#ifndef __vector_memory_H
+#define __vector_memory_H
+/*---------------------------- vector_memory.h ---------------------------*/
+
+
+
+
+/**
+ * Memory management for vectors. This class is used by all
+ * iterative methods to allocate space for auxilliary
+ * vectors. This class is used to avoid interaction with the
+ * operating system whenever a vector is needed. Especially, when
+ * an iterative method is invoked as part of an outer iteration,
+ * this would lead to runtime overhead and memory fragmentation.
+ *
+ * Classes derived from this class implement a more or less
+ * sophisticated management of vectors. One of these has to be
+ * applied by the user according to his needs.
+ */
+template<class Vector>
+class VectorMemory {
+ public:
+ /**
+ * Return new vector from the pool.
+ */
+ virtual Vector* alloc() = 0;
+
+ /**
+ * Return a vector into the pool
+ * for later use.
+ */
+ virtual void free(const Vector*) = 0;
+};
+
+
+
+
+
+/**
+ * Simple memory management. This memory class is just made for
+ * tests. It requires the vector to have a constructor with integer
+ * parameter indicating the size. It just allocates and deletes
+ * vectors as needed from the global heap, i.e. performs no
+ * specially adapted actions to the purpose of this class.
+ */
+template<class Vector>
+class PrimitiveVectorMemory : public VectorMemory<Vector> {
+ const unsigned int size;
+ public:
+ /**
+ * Constructor.
+ * #sz# is the length of vectors to
+ * be generated by alloc.
+ */
+ PrimitiveVectorMemory (const unsigned int sz) : size(sz) {};
+
+ /**
+ * Allocate a vector of the given size
+ * from the global heap.
+ */
+ virtual Vector* alloc() {
+ return new Vector(size);
+ };
+
+ /**
+ * Return a vector to the global heap.
+ */
+ virtual void free (const Vector* v) {
+ delete v;
+ };
+};
+
+
+
+/*---------------------------- vector_memory.h ---------------------------*/
+/* end of #ifndef __vector_memory_H */
+#endif
+/*---------------------------- vector_memory.h ---------------------------*/
--- /dev/null
+// $Id$
+
+
+#include <lac/solver_control.h>
+
+
+
+/*----------------------- SolverControl ---------------------------------*/
+
+
+SolverControl::SolverControl (const unsigned int maxiter,
+ const double tolerance) :
+ maxsteps(maxiter),
+ tol(tolerance),
+ lvalue(1.e300),
+ lstep(0)
+{};
+
+
+
+SolverControl::State
+SolverControl::check (const unsigned int step,
+ const double check_value) {
+ lstep = step;
+ lvalue = check_value;
+ if (step>=maxsteps) return failure;
+ if (check_value <= tol) return success;
+ return iterate;
+};
+
+
+
+double
+SolverControl::last_value() const
+{
+ return lvalue;
+};
+
+
+
+unsigned int
+SolverControl::last_step() const
+{
+ return lstep;
+};
+
+
+
+
+/*----------------------- ReductionControl ---------------------------------*/
+
+
+
+ReductionControl::ReductionControl(const unsigned int n,
+ const double tol,
+ const double red) :
+ SolverControl (n, tol),
+ reduce(red)
+{};
+
+
+
+
+double
+ReductionControl::initial_value() const {
+ return initial_val;
+};
+
+
+
+SolverControl::State
+ReductionControl::check (const unsigned int step,
+ const double check_value)
+{
+ if (step==0)
+ {
+ initial_val = check_value;
+ reduced_tol = check_value * reduce;
+ };
+
+ if (check_value < reduced_tol)
+ return success;
+ else
+ return SolverControl::check(step, check_value);
+};
+