*
* @author Guido Kanschat, 2000
*/
-template <class MATRIX = SparseMatrix<double>,
- class VECTOR = Vector<double> >
-class EigenPower : private Solver<MATRIX,VECTOR>
+template <class VECTOR = Vector<double> >
+class EigenPower : private Solver<VECTOR>
{
public:
/**
* corresponding eigenvector,
* normalized with respect to the l2-norm.
*/
- typename Solver<MATRIX,VECTOR>::ReturnState
+ template <class MATRIX>
+ typename Solver<VECTOR>::ReturnState
solve (double &value,
const MATRIX &A,
VECTOR &x);
*
* This class implements an adaptive version of the inverse iteration by Wieland.
*
+ * There are two choices for the stopping criterion: by default, the
+ * norm of the residual $A x - l x$ is computed. Since this might not
+ * converge to zero for non-symmetric matrices with non-trivial Jordan
+ * blocks, it can be replaced by checking the difference of successive
+ * eigenvalues. Use @p{AdditionalData::use_residual} for switching
+ * these options.
+ *
* @author Guido Kanschat, 2000
*/
-template <class MATRIX = SparseMatrix<double>,
- class VECTOR = Vector<double> >
-class EigenInverse : private Solver<MATRIX, VECTOR>
+template <class VECTOR = Vector<double> >
+class EigenInverse : private Solver<VECTOR>
{
public:
/**
* Standardized data struct to
* pipe additional data to the
- * solver. This solver does not
- * need additional data yet.
+ * solver.
*/
struct AdditionalData
- {};
+ {
+ /**
+ * Flag for the stopping criterion.
+ */
+ bool use_residual;
+ /**
+ * Constructor.
+ */
+ AdditionalData (bool use_residual = true):
+ use_residual(use_residual)
+ {}
+
+ };
/**
* Constructor.
* normalized with respect to the
* l2-norm.
*/
- typename Solver<MATRIX,VECTOR>::ReturnState
+ template <class MATRIX>
+ typename Solver<VECTOR>::ReturnState
solve (double &value,
const MATRIX &A,
VECTOR &x);
//----------------------------------------------------------------------
-template <class MATRIX, class VECTOR>
-EigenPower<MATRIX, VECTOR>::EigenPower (SolverControl &cn,
- VectorMemory<VECTOR> &mem,
- const AdditionalData &data):
- Solver<MATRIX, VECTOR>(cn, mem),
+template <class VECTOR>
+EigenPower<VECTOR>::EigenPower (SolverControl &cn,
+ VectorMemory<VECTOR> &mem,
+ const AdditionalData &data):
+ Solver<VECTOR>(cn, mem),
additional_data(data)
{}
-template <class MATRIX, class VECTOR>
-EigenPower<MATRIX, VECTOR>::~EigenPower ()
+template <class VECTOR>
+EigenPower<VECTOR>::~EigenPower ()
{}
-template <class MATRIX, class VECTOR>
-typename Solver<MATRIX,VECTOR>::ReturnState
-EigenPower<MATRIX, VECTOR>::solve (double &value,
- const MATRIX &A,
- VECTOR &x)
+template <class VECTOR>
+template <class MATRIX>
+typename Solver<VECTOR>::ReturnState
+EigenPower<VECTOR>::solve (double &value,
+ const MATRIX &A,
+ VECTOR &x)
{
SolverControl::State conv=SolverControl::iterate;
//----------------------------------------------------------------------//
-template <class MATRIX, class VECTOR>
-EigenInverse<MATRIX, VECTOR>::EigenInverse (SolverControl &cn,
+template <class VECTOR>
+EigenInverse<VECTOR>::EigenInverse (SolverControl &cn,
VectorMemory<VECTOR> &mem,
const AdditionalData &data):
- Solver<MATRIX, VECTOR>(cn, mem),
+ Solver<VECTOR>(cn, mem),
additional_data(data)
{}
-template <class MATRIX, class VECTOR>
-EigenInverse<MATRIX, VECTOR>::~EigenInverse ()
+template <class VECTOR>
+EigenInverse<VECTOR>::~EigenInverse ()
{}
-template <class MATRIX, class VECTOR>
-typename Solver<MATRIX,VECTOR>::ReturnState
-EigenInverse<MATRIX, VECTOR>::solve (double &value,
- const MATRIX &A,
- VECTOR &x)
+template <class VECTOR>
+template <class MATRIX>
+typename Solver<VECTOR>::ReturnState
+EigenInverse<VECTOR>::solve (double &value,
+ const MATRIX &A,
+ VECTOR &x)
{
deallog.push("Wieland");
// Define solver
ReductionControl inner_control (100, 1.e-16, 1.e-8, false, false);
PreconditionIdentity prec;
- SolverQMRS<ShiftedMatrix <MATRIX>, VECTOR>
+ SolverCG<VECTOR>
solver(inner_control, memory);
// Next step for recomputing the shift
// Update normalized eigenvector
x.equ (1./length, y);
// Compute residual
- y.equ (value, x);
- double res = A.residual (r,x,y);
-
- // Check the residual
- conv = control().check (iter, res);
+ if (additional_data.use_residual)
+ {
+ y.equ (value, x);
+ double res = A.residual (r,x,y);
+ // Check the residual
+ conv = control().check (iter, res);
+ } else {
+ conv = control().check (iter, fabs(1.-old_value/value));
+ }
old_value = value;
}
*
* @author Wolfgang Bangerth, Guido Kanschat, Ralf Hartmann, 1997, 1998, 1999
*/
-template <class Matrix = SparseMatrix<double>,
- class Vector = Vector<double> >
+template <class Vector = Vector<double> >
class Solver
{
public:
/*-------------------------------- Inline functions ------------------------*/
-template <class Matrix, class Vector>
+template <class Vector>
inline
-SolverControl & Solver<Matrix,Vector>::control() const
+SolverControl &
+Solver<Vector>::control() const
{
return cntrl;
};
-template<class Matrix, class Vector>
+template<class Vector>
inline
-Solver<Matrix, Vector>::Solver(SolverControl &cn, VectorMemory<Vector> &mem)
+Solver<Vector>::Solver(SolverControl &cn, VectorMemory<Vector> &mem)
: cntrl(cn),
memory(mem)
{};
* has a default argument, so you may call it without the additional
* parameter.
*/
-template <class Matrix = SparseMatrix<double>, class Vector = Vector<double> >
-class SolverBicgstab : public Solver<Matrix,Vector>
+template <class VECTOR = Vector<double> >
+class SolverBicgstab : private Solver<VECTOR>
{
public:
/**
* Constructor.
*/
SolverBicgstab (SolverControl &cn,
- VectorMemory<Vector> &mem,
+ VectorMemory<VECTOR> &mem,
const AdditionalData &data=AdditionalData());
/**
/**
* Solve primal problem only.
*/
- template<class Preconditioner>
- typename Solver<Matrix,Vector>::ReturnState
- solve (const Matrix &A,
- Vector &x,
- const Vector &b,
+ template<class MATRIX, class Preconditioner>
+ typename Solver<VECTOR>::ReturnState
+ solve (const MATRIX &A,
+ VECTOR &x,
+ const VECTOR &b,
const Preconditioner& precondition);
protected:
/**
* Computation of the stopping criterion.
*/
- virtual double criterion();
+ template <class MATRIX>
+ double criterion (const MATRIX& A, const VECTOR& x, const VECTOR& b);
/**
* Interface for derived class.
* convergence history.
*/
virtual void print_vectors(const unsigned int step,
- const Vector& x,
- const Vector& r,
- const Vector& d) const;
+ const VECTOR& x,
+ const VECTOR& r,
+ const VECTOR& d) const;
/**
* Auxiliary vector.
*/
- Vector *Vx;
+ VECTOR *Vx;
/**
* Auxiliary vector.
*/
- Vector *Vr;
+ VECTOR *Vr;
/**
* Auxiliary vector.
*/
- Vector *Vrbar;
+ VECTOR *Vrbar;
/**
* Auxiliary vector.
*/
- Vector *Vp;
+ VECTOR *Vp;
/**
* Auxiliary vector.
*/
- Vector *Vy;
+ VECTOR *Vy;
/**
* Auxiliary vector.
*/
- Vector *Vz;
+ VECTOR *Vz;
/**
* Auxiliary vector.
*/
- Vector *Vs;
+ VECTOR *Vs;
/**
* Auxiliary vector.
*/
- Vector *Vt;
+ VECTOR *Vt;
/**
* Auxiliary vector.
*/
- Vector *Vv;
+ VECTOR *Vv;
/**
* Right hand side vector.
*/
- const Vector *Vb;
-
- /**
- * Pointer to the system matrix.
- */
- const Matrix *MA;
+ const VECTOR *Vb;
/**
* Auxiliary value.
/**
* Everything before the iteration loop.
*/
- SolverControl::State start();
+ template <class MATRIX>
+ SolverControl::State start(const MATRIX& A);
/**
* The iteration loop itself.
*/
- template<class Preconditioner>
- typename Solver<Matrix,Vector>::ReturnState
- iterate(const Preconditioner& precondition);
+ template<class MATRIX, class PRECONDITIONER>
+ typename Solver<VECTOR>::ReturnState
+ iterate(const MATRIX& A, const PRECONDITIONER& precondition);
};
/*-------------------------Inline functions -------------------------------*/
-template<class Matrix, class Vector>
-SolverBicgstab<Matrix, Vector>::SolverBicgstab (SolverControl &cn,
- VectorMemory<Vector> &mem,
- const AdditionalData &)
+template<class VECTOR>
+SolverBicgstab<VECTOR>::SolverBicgstab (SolverControl &cn,
+ VectorMemory<VECTOR> &mem,
+ const AdditionalData &)
:
- Solver<Matrix,Vector>(cn,mem)
+ Solver<VECTOR>(cn,mem)
{}
-template<class Matrix, class Vector>
-SolverBicgstab<Matrix, Vector>::~SolverBicgstab ()
+template<class VECTOR>
+SolverBicgstab<VECTOR>::~SolverBicgstab ()
{}
-template<class Matrix, class Vector>
+template <class VECTOR>
+template <class MATRIX>
double
-SolverBicgstab<Matrix, Vector>::criterion()
+SolverBicgstab<VECTOR>::criterion (const MATRIX& A, const VECTOR& x, const VECTOR& b)
{
- res = MA->residual(*Vt, *Vx, *Vb);
+ res = A.residual(*Vt, x, b);
return res;
}
-template < class Matrix, class Vector >
+template <class VECTOR >
+template <class MATRIX>
SolverControl::State
-SolverBicgstab<Matrix, Vector>::start()
+SolverBicgstab<VECTOR>::start(const MATRIX& A)
{
- res = MA->residual(*Vr, *Vx, *Vb);
+ res = A.residual(*Vr, *Vx, *Vb);
Vp->reinit(*Vx);
Vv->reinit(*Vx);
Vrbar->equ(1.,*Vr);
-template<class Matrix, class Vector>
+template<class VECTOR>
void
-SolverBicgstab<Matrix,Vector>::print_vectors(const unsigned int,
- const Vector&,
- const Vector&,
- const Vector&) const
+SolverBicgstab<VECTOR>::print_vectors(const unsigned int,
+ const VECTOR&,
+ const VECTOR&,
+ const VECTOR&) const
{}
-template<class Matrix, class Vector>
-template<class Preconditioner>
-typename Solver<Matrix,Vector>::ReturnState
-SolverBicgstab<Matrix, Vector>::iterate(const Preconditioner& precondition)
+template<class VECTOR>
+template<class MATRIX, class PRECONDITIONER>
+typename Solver<VECTOR>::ReturnState
+SolverBicgstab<VECTOR>::iterate(const MATRIX& A,
+ const PRECONDITIONER& precondition)
{
SolverControl::State state = SolverControl::iterate;
alpha = omega = rho = 1.;
- Vector& r = *Vr;
- Vector& rbar = *Vrbar;
- Vector& p = *Vp;
- Vector& y = *Vy;
- Vector& z = *Vz;
- Vector& s = *Vs;
- Vector& t = *Vt;
- Vector& v = *Vv;
+ VECTOR& r = *Vr;
+ VECTOR& rbar = *Vrbar;
+ VECTOR& p = *Vp;
+ VECTOR& y = *Vy;
+ VECTOR& z = *Vz;
+ VECTOR& s = *Vs;
+ VECTOR& t = *Vt;
+ VECTOR& v = *Vv;
do
{
rho = rhobar;
p.sadd(beta, 1., r, -beta*omega, v);
precondition(y,p);
- MA->vmult(v,y);
+ A.vmult(v,y);
rhobar = rbar * v;
alpha = rho/rhobar;
s.equ(1., r, -alpha, v);
precondition(z,s);
- MA->vmult(t,z);
+ A.vmult(t,z);
rhobar = t*s;
omega = rhobar/(t*t);
Vx->add(alpha, y, omega, z);
r.equ(1., s, -omega, t);
- state = control().check(++step, criterion());
+ state = control().check(++step, criterion(A, *Vx, *Vb));
print_vectors(step, *Vx, r, y);
}
while (state == SolverControl::iterate);
}
-template<class Matrix, class Vector>
-template<class Preconditioner>
-typename Solver<Matrix,Vector>::ReturnState
-SolverBicgstab<Matrix, Vector>::solve(const Matrix &A,
- Vector &x,
- const Vector &b,
- const Preconditioner& precondition)
+template<class VECTOR>
+template<class MATRIX, class PRECONDITIONER>
+typename Solver<VECTOR>::ReturnState
+SolverBicgstab<VECTOR>::solve(const MATRIX &A,
+ VECTOR &x,
+ const VECTOR &b,
+ const PRECONDITIONER& precondition)
{
deallog.push("Bicgstab");
Vr = memory.alloc(); Vr->reinit(x);
Vt = memory.alloc(); Vt->reinit(x);
Vv = memory.alloc();
- MA = &A;
Vx = &x;
Vb = &b;
{
if (step)
deallog << "Restart step " << step << endl;
- if (start() == SolverControl::success) break;
- state = iterate(precondition);
+ if (start(A) == SolverControl::success) break;
+ state = iterate(A, 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 = SparseMatrix<double>, class Vector = Vector<double> >
-class SolverCG : public Solver<Matrix,Vector>
+template <class VECTOR = Vector<double> >
+class SolverCG : private Solver<VECTOR>
{
public:
/**
* Constructor.
*/
SolverCG (SolverControl &cn,
- VectorMemory<Vector> &mem,
+ VectorMemory<VECTOR> &mem,
const AdditionalData &data=AdditionalData());
/**
/**
* Solver method.
*/
- template<class Preconditioner>
- typename Solver<Matrix,Vector>::ReturnState
- solve (const Matrix &A,
- Vector &x,
- const Vector &b,
+ template<class MATRIX, class Preconditioner>
+ typename Solver<VECTOR>::ReturnState
+ solve (const MATRIX &A,
+ VECTOR &x,
+ const VECTOR &b,
const Preconditioner& precondition);
protected:
* convergence history.
*/
virtual void print_vectors(const unsigned int step,
- const Vector& x,
- const Vector& r,
- const Vector& d) const;
+ const VECTOR& x,
+ const VECTOR& r,
+ const VECTOR& d) const;
/**
* Temporary vectors, allocated through
* of the actual solution process and
* deallocated at the end.
*/
- Vector *Vr;
- Vector *Vp;
- Vector *Vz;
- Vector *VAp;
+ VECTOR *Vr;
+ VECTOR *Vp;
+ VECTOR *Vz;
+ VECTOR *VAp;
/**
* Within the iteration loop, the
/*------------------------- Implementation ----------------------------*/
-template<class Matrix, class Vector>
-SolverCG<Matrix,Vector>::SolverCG(SolverControl &cn,
- VectorMemory<Vector> &mem,
- const AdditionalData &)
+template<class VECTOR>
+SolverCG<VECTOR>::SolverCG(SolverControl &cn,
+ VectorMemory<VECTOR> &mem,
+ const AdditionalData &)
:
- Solver<Matrix,Vector>(cn,mem)
+ Solver<VECTOR>(cn,mem)
{}
-template<class Matrix, class Vector>
-SolverCG<Matrix,Vector>::~SolverCG ()
+template<class VECTOR>
+SolverCG<VECTOR>::~SolverCG ()
{}
-template<class Matrix, class Vector>
+template<class VECTOR>
double
-SolverCG<Matrix,Vector>::criterion()
+SolverCG<VECTOR>::criterion()
{
return sqrt(res2);
}
-template<class Matrix, class Vector>
+template<class VECTOR>
void
-SolverCG<Matrix,Vector>::print_vectors(const unsigned int,
- const Vector&,
- const Vector&,
- const Vector&) const
+SolverCG<VECTOR>::print_vectors(const unsigned int,
+ const VECTOR&,
+ const VECTOR&,
+ const VECTOR&) const
{}
-template<class Matrix, class Vector>
-template<class Preconditioner>
-typename Solver<Matrix,Vector>::ReturnState
-SolverCG<Matrix,Vector>::solve (const Matrix &A,
- Vector &x,
- const Vector &b,
- const Preconditioner& precondition)
+template<class VECTOR>
+template<class MATRIX, class Preconditioner>
+typename Solver<VECTOR>::ReturnState
+SolverCG<VECTOR>::solve (const MATRIX &A,
+ VECTOR &x,
+ const VECTOR &b,
+ const Preconditioner& precondition)
{
SolverControl::State conv=SolverControl::iterate;
Vz = memory.alloc();
VAp = memory.alloc();
// define some aliases for simpler access
- Vector& g = *Vr;
- Vector& h = *Vp;
- Vector& d = *Vz;
- Vector& Ad = *VAp;
+ VECTOR& g = *Vr;
+ VECTOR& h = *Vp;
+ VECTOR& d = *Vz;
+ VECTOR& Ad = *VAp;
// resize the vectors, but do not set
// the values since they'd be overwritten
// soon anyway.
*
* @author Wolfgang Bangerth
*/
-template <class MATRIX = SparseMatrix<double>,
- class VECTOR = Vector<double> >
-class SolverGMRES : public Solver<MATRIX, VECTOR>
+template <class VECTOR = Vector<double> >
+class SolverGMRES : private Solver<VECTOR>
{
public:
/**
/**
* Solver method.
*/
- template<class Preconditioner>
- typename Solver<MATRIX,VECTOR>::ReturnState solve (const MATRIX &A,
+ template<class MATRIX, class PRECONDITIONER>
+ typename Solver<VECTOR>::ReturnState solve (const MATRIX &A,
VECTOR &x,
const VECTOR &b,
- const Preconditioner& precondition);
+ const PRECONDITIONER& precondition);
DeclException1 (ExcTooFewTmpVectors,
int,
/* --------------------- Inline and template functions ------------------- */
-template <class MATRIX, class VECTOR>
-SolverGMRES<MATRIX,VECTOR>::SolverGMRES (SolverControl &cn,
+template <class VECTOR>
+SolverGMRES<VECTOR>::SolverGMRES (SolverControl &cn,
VectorMemory<VECTOR> &mem,
const AdditionalData &data) :
- Solver<MATRIX,VECTOR> (cn,mem),
+ Solver<VECTOR> (cn,mem),
additional_data(data)
{
Assert (additional_data.max_n_tmp_vectors >= 10,
};
-template <class MATRIX, class VECTOR>
+template <class VECTOR>
inline
void
-SolverGMRES<MATRIX,VECTOR>::givens_rotation (Vector<double> &h,
- Vector<double> &b,
- Vector<double> &ci,
- Vector<double> &si,
- int col) const
+SolverGMRES<VECTOR>::givens_rotation (Vector<double> &h,
+ Vector<double> &b,
+ Vector<double> &ci,
+ Vector<double> &si,
+ int col) const
{
for (int i=0 ; i<col ; i++)
{
}
-template<class MATRIX, class VECTOR>
-template<class Preconditioner>
-typename Solver<MATRIX,VECTOR>::ReturnState
-SolverGMRES<MATRIX,VECTOR>::solve (const MATRIX& A,
+template<class VECTOR>
+template<class MATRIX, class PRECONDITIONER>
+typename Solver<VECTOR>::ReturnState
+SolverGMRES<VECTOR>::solve (const MATRIX& A,
VECTOR & x,
const VECTOR& b,
- const Preconditioner& precondition)
+ const PRECONDITIONER& precondition)
{
// this code was written by the fathers of
// DEAL. I take absolutely no guarantees
};
-template<class MATRIX, class VECTOR>
+template<class VECTOR>
double
-SolverGMRES<MATRIX,VECTOR>::criterion ()
+SolverGMRES<VECTOR>::criterion ()
{
// dummy implementation. this function is
// not needed for the present implementation
*
* @author Thomas Richter, 2000
*/
-template <class Matrix = SparseMatrix<double>, class Vector = Vector<double> >
-class SolverMinRes : public Solver<Matrix,Vector>
+template <class VECTOR = Vector<double> >
+class SolverMinRes : private Solver<VECTOR>
{
public:
/**
* Constructor.
*/
SolverMinRes (SolverControl &cn,
- VectorMemory<Vector> &mem,
+ VectorMemory<VECTOR> &mem,
const AdditionalData &data=AdditionalData());
/**
/**
* Solver method.
*/
- template<class Preconditioner>
- typename Solver<Matrix,Vector>::ReturnState
- solve (const Matrix &A,
- Vector &x,
- const Vector &b,
+ template<class MATRIX, class Preconditioner>
+ typename Solver<VECTOR>::ReturnState
+ solve (const MATRIX &A,
+ VECTOR &x,
+ const VECTOR &b,
const Preconditioner& precondition);
/**
* convergence history.
*/
virtual void print_vectors(const unsigned int step,
- const Vector& x,
- const Vector& r,
- const Vector& d) const;
+ const VECTOR& x,
+ const VECTOR& r,
+ const VECTOR& d) const;
/**
* Temporary vectors, allocated through
* of the actual solution process and
* deallocated at the end.
*/
- Vector *Vu0, *Vu1, *Vu2;
- Vector *Vm0, *Vm1, *Vm2;
- Vector *Vv;
+ VECTOR *Vu0, *Vu1, *Vu2;
+ VECTOR *Vm0, *Vm1, *Vm2;
+ VECTOR *Vv;
/**
* Within the iteration loop, the
/*------------------------- Implementation ----------------------------*/
-template<class Matrix, class Vector>
-SolverMinRes<Matrix,Vector>::SolverMinRes (SolverControl &cn,
- VectorMemory<Vector> &mem,
- const AdditionalData &)
+template<class VECTOR>
+SolverMinRes<VECTOR>::SolverMinRes (SolverControl &cn,
+ VectorMemory<VECTOR> &mem,
+ const AdditionalData &)
:
- Solver<Matrix,Vector>(cn,mem)
+ Solver<VECTOR>(cn,mem)
{}
-template<class Matrix, class Vector>
-SolverMinRes<Matrix,Vector>::~SolverMinRes ()
+template<class VECTOR>
+SolverMinRes<VECTOR>::~SolverMinRes ()
{}
-template<class Matrix, class Vector>
+template<class VECTOR>
double
-SolverMinRes<Matrix,Vector>::criterion()
+SolverMinRes<VECTOR>::criterion()
{
return res2;
};
-template<class Matrix, class Vector>
+template<class VECTOR>
void
-SolverMinRes<Matrix,Vector>::print_vectors(const unsigned int,
- const Vector&,
- const Vector&,
- const Vector&) const
+SolverMinRes<VECTOR>::print_vectors(const unsigned int,
+ const VECTOR&,
+ const VECTOR&,
+ const VECTOR&) const
{}
-template<class Matrix, class Vector>
-template<class Preconditioner>
-typename Solver<Matrix,Vector>::ReturnState
-SolverMinRes<Matrix,Vector>::solve (const Matrix &A,
- Vector &x,
- const Vector &b,
- const Preconditioner& precondition)
+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)
{
SolverControl::State conv=SolverControl::iterate;
Vm1 = memory.alloc();
Vm2 = memory.alloc();
// define some aliases for simpler access
- typedef Vector vecref;
+ typedef VECTOR vecref;
vecref u[3] = {*Vu0, *Vu1, *Vu2};
vecref m[3] = {*Vm0, *Vm1, *Vm2};
vecref v = *Vv;
*
* @author Guido Kanschat, 1999
*/
-template <class Matrix = SparseMatrix<double>, class Vector = Vector<double> >
-class SolverQMRS : public Solver<Matrix,Vector>
+template <class VECTOR = Vector<double> >
+class SolverQMRS : private Solver<VECTOR>
{
public:
/**
* Constructor.
*/
SolverQMRS (SolverControl &cn,
- VectorMemory<Vector> &mem,
+ VectorMemory<VECTOR> &mem,
const AdditionalData &data=AdditionalData());
/**
* Solver method.
*/
- template<class Preconditioner>
- typename Solver<Matrix,Vector>::ReturnState
- solve (const Matrix &A,
- Vector &x,
- const Vector &b,
+ template<class MATRIX, class Preconditioner>
+ typename Solver<VECTOR>::ReturnState
+ solve (const MATRIX &A,
+ VECTOR &x,
+ const VECTOR &b,
const Preconditioner& precondition);
/**
* convergence history.
*/
virtual void print_vectors(const unsigned int step,
- const Vector& x,
- const Vector& r,
- const Vector& d) const;
+ const VECTOR& x,
+ const VECTOR& r,
+ const VECTOR& d) const;
protected:
/**
* Implementation of the computation of
* of the actual solution process and
* deallocated at the end.
*/
- Vector *Vv;
- Vector *Vp;
- Vector *Vq;
- Vector *Vt;
- Vector *Vd;
+ VECTOR *Vv;
+ VECTOR *Vp;
+ VECTOR *Vq;
+ VECTOR *Vt;
+ VECTOR *Vd;
/**
* Iteration vector.
*/
- Vector *Vx;
+ VECTOR *Vx;
/**
* RHS vector.
*/
- const Vector *Vb;
+ const VECTOR *Vb;
- /**
- * Pointer to the matrix to be inverted.
- */
- const Matrix* MA;
/**
* Within the iteration loop, the
* square of the residual vector is
/**
* The iteration loop itself.
*/
- template<class Preconditioner>
- typename Solver<Matrix,Vector>::ReturnState
- iterate(const Preconditioner& precondition);
+ template<class MATRIX, class Preconditioner>
+ typename Solver<VECTOR>::ReturnState
+ iterate(const MATRIX& A, const Preconditioner& precondition);
/**
* The current iteration step.
*/
/*------------------------- Implementation ----------------------------*/
-template<class Matrix, class Vector>
-SolverQMRS<Matrix,Vector>::SolverQMRS(SolverControl &cn,
- VectorMemory<Vector> &mem,
- const AdditionalData &data) :
- Solver<Matrix,Vector>(cn,mem),
+template<class VECTOR>
+SolverQMRS<VECTOR>::SolverQMRS(SolverControl &cn,
+ VectorMemory<VECTOR> &mem,
+ const AdditionalData &data) :
+ Solver<VECTOR>(cn,mem),
additional_data(data)
{};
-template<class Matrix, class Vector>
+template<class VECTOR>
double
-SolverQMRS<Matrix,Vector>::criterion()
+SolverQMRS<VECTOR>::criterion()
{
return sqrt(res2);
}
-template<class Matrix, class Vector>
+template<class VECTOR>
void
-SolverQMRS<Matrix,Vector>::print_vectors(const unsigned int,
- const Vector&,
- const Vector&,
- const Vector&) const
+SolverQMRS<VECTOR>::print_vectors(const unsigned int,
+ const VECTOR&,
+ const VECTOR&,
+ const VECTOR&) const
{}
-template<class Matrix, class Vector>
-template<class Preconditioner>
-typename Solver<Matrix,Vector>::ReturnState
-SolverQMRS<Matrix,Vector>::solve (const Matrix &A,
- Vector &x,
- const Vector &b,
- const Preconditioner& precondition)
+template<class VECTOR>
+template<class MATRIX, class Preconditioner>
+typename Solver<VECTOR>::ReturnState
+SolverQMRS<VECTOR>::solve (const MATRIX &A,
+ VECTOR &x,
+ const VECTOR &b,
+ const Preconditioner& precondition)
{
deallog.push("QMRS");
Vt = memory.alloc();
Vd = memory.alloc();
- MA = &A;
Vx = &x;
Vb = &b;
// resize the vectors, but do not set
{
if (step)
deallog << "Restart step " << step << endl;
- state = iterate(precondition);
+ state = iterate(A, precondition);
}
while (state == breakdown);
return state;
};
-template<class Matrix, class Vector>
-template<class Preconditioner>
-typename Solver<Matrix,Vector>::ReturnState
-SolverQMRS<Matrix,Vector>::iterate(const Preconditioner& precondition)
+template<class VECTOR>
+template<class MATRIX, class Preconditioner>
+typename Solver<VECTOR>::ReturnState
+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.
SolverControl::State state = SolverControl::iterate;
// define some aliases for simpler access
- Vector& v = *Vv;
- Vector& p = *Vp;
- Vector& q = *Vq;
- Vector& t = *Vt;
- Vector& d = *Vd;
- Vector& x = *Vx;
- const Vector& b = *Vb;
-
- const Matrix& A = *MA;
+ VECTOR& v = *Vv;
+ VECTOR& p = *Vp;
+ VECTOR& q = *Vq;
+ VECTOR& t = *Vt;
+ VECTOR& d = *Vd;
+ VECTOR& x = *Vx;
+ const VECTOR& b = *Vb;
int it=0;
*
* @author Ralf Hartmann
*/
-template <class Matrix = SparseMatrix<double>, class Vector = Vector<double> >
-class SolverRichardson : public Solver<Matrix, Vector>
+template <class VECTOR = Vector<double> >
+class SolverRichardson : private Solver<VECTOR>
{
public:
/**
* Constructor.
*/
SolverRichardson (SolverControl &cn,
- VectorMemory<Vector> &mem,
+ VectorMemory<VECTOR> &mem,
const AdditionalData &data=AdditionalData());
/**
* Solve $Ax=b$ for $x$.
*/
- template<class Preconditioner>
- typename Solver<Matrix,Vector>::ReturnState solve (const Matrix &A,
- Vector &x,
- const Vector &b,
+ template<class MATRIX, class Preconditioner>
+ typename Solver<VECTOR>::ReturnState solve (const MATRIX &A,
+ VECTOR &x,
+ const VECTOR &b,
const Preconditioner& precondition);
/**
* convergence history.
*/
virtual void print_vectors(const unsigned int step,
- const Vector& x,
- const Vector& r,
- const Vector& d) const;
+ const VECTOR& x,
+ const VECTOR& r,
+ const VECTOR& d) const;
protected:
/**
* of the actual solution process and
* deallocated at the end.
*/
- Vector *Vr;
- Vector *Vd;
+ VECTOR *Vr;
+ VECTOR *Vd;
/**
* Damping parameter.
/*----------------- Implementation of the Richardson Method ------------------*/
-template<class Matrix, class Vector>
-SolverRichardson<Matrix,Vector>::SolverRichardson(SolverControl &cn,
- VectorMemory<Vector> &mem,
- const AdditionalData &data):
- Solver<Matrix,Vector> (cn,mem),
+template<class VECTOR>
+SolverRichardson<VECTOR>::SolverRichardson(SolverControl &cn,
+ VectorMemory<VECTOR> &mem,
+ const AdditionalData &data):
+ Solver<VECTOR> (cn,mem),
additional_data(data) {};
-template<class Matrix,class Vector>
-template<class Preconditioner>
-typename Solver<Matrix,Vector>::ReturnState
-SolverRichardson<Matrix,Vector>::solve (const Matrix &A,
- Vector &x,
- const Vector &b,
- const Preconditioner& precondition)
+template<class VECTOR>
+template<class MATRIX, class Preconditioner>
+typename Solver<VECTOR>::ReturnState
+SolverRichardson<VECTOR>::solve (const MATRIX &A,
+ VECTOR &x,
+ const VECTOR &b,
+ const Preconditioner& precondition)
{
SolverControl::State conv=SolverControl::iterate;
// Memory allocation
- Vr = memory.alloc(); Vector& r = *Vr; r.reinit(x);
- Vd = memory.alloc(); Vector& d = *Vd; d.reinit(x);
+ Vr = memory.alloc(); VECTOR& r = *Vr; r.reinit(x);
+ Vd = memory.alloc(); VECTOR& d = *Vd; d.reinit(x);
deallog.push("Richardson");
}
-template<class Matrix, class Vector>
+template<class VECTOR>
void
-SolverRichardson<Matrix,Vector>::print_vectors(const unsigned int,
- const Vector&,
- const Vector&,
- const Vector&) const
+SolverRichardson<VECTOR>::print_vectors(const unsigned int,
+ const VECTOR&,
+ const VECTOR&,
+ const VECTOR&) const
{}
-template<class Matrix,class Vector>
+template<class VECTOR>
inline double
-SolverRichardson<Matrix,Vector>::criterion()
+SolverRichardson<VECTOR>::criterion()
{
return res;
}
-template<class Matrix,class Vector>
+template<class VECTOR>
inline void
-SolverRichardson<Matrix,Vector>::set_omega(double om)
+SolverRichardson<VECTOR>::set_omega(double om)
{
additional_data.omega=om;
}
*
* @author Ralf Hartmann, 1999
*/
-template <class Matrix = SparseMatrix<double>,
- class Vector = Vector<double> >
+template <class Vector = Vector<double> >
class SolverSelector
{
public:
* was specified in the constructor.
*
*/
- template<class Preconditioner>
- typename Solver<Matrix,Vector>::ReturnState solve(const Matrix &A,
+ template<class Matrix, class Preconditioner>
+ typename Solver<Vector>::ReturnState solve(const Matrix &A,
Vector &x,
const Vector &b,
const Preconditioner &precond) const;
* Set the additional data. For more
* info see the #Solver# class.
*/
- void set_data(const typename SolverRichardson<Matrix,Vector>
+ void set_data(const typename SolverRichardson<Vector>
::AdditionalData &data);
/**
* Set the additional data. For more
* info see the #Solver# class.
*/
- void set_data(const typename SolverCG<Matrix,Vector>
+ void set_data(const typename SolverCG<Vector>
::AdditionalData &data);
/**
* Set the additional data. For more
* info see the #Solver# class.
*/
- void set_data(const typename SolverBicgstab<Matrix,Vector>
+ void set_data(const typename SolverBicgstab<Vector>
::AdditionalData &data);
/**
* Set the additional data. For more
* info see the #Solver# class.
*/
- void set_data(const typename SolverGMRES<Matrix,Vector>
+ void set_data(const typename SolverGMRES<Vector>
::AdditionalData &data);
/**
/**
* Stores the additional data.
*/
- typename SolverRichardson<Matrix,Vector>::AdditionalData richardson_data;
+ typename SolverRichardson<Vector>::AdditionalData richardson_data;
/**
* Stores the additional data.
*/
- typename SolverCG<Matrix,Vector>::AdditionalData cg_data;
+ typename SolverCG<Vector>::AdditionalData cg_data;
/**
* Stores the additional data.
*/
- typename SolverBicgstab<Matrix,Vector>::AdditionalData bicgstab_data;
+ typename SolverBicgstab<Vector>::AdditionalData bicgstab_data;
/**
* Stores the additional data.
*/
- typename SolverGMRES<Matrix,Vector>::AdditionalData gmres_data;
+ typename SolverGMRES<Vector>::AdditionalData gmres_data;
};
/* --------------------- Inline and template functions ------------------- */
-template <class Matrix, class Vector>
-SolverSelector<Matrix, Vector>::SolverSelector(string solver_name,
- SolverControl &control,
- VectorMemory<Vector> &vector_memory) :
+template <class Vector>
+SolverSelector<Vector>::SolverSelector(string solver_name,
+ SolverControl &control,
+ VectorMemory<Vector> &vector_memory) :
solver_name(solver_name),
control(&control),
vector_memory(&vector_memory)
{};
-template <class Matrix, class Vector>
-SolverSelector<Matrix, Vector>::~SolverSelector()
+template <class Vector>
+SolverSelector<Vector>::~SolverSelector()
{};
-template <class Matrix, class Vector>
-template<class Preconditioner>
-typename Solver<Matrix,Vector>::ReturnState
-SolverSelector<Matrix, Vector>::solve(const Matrix &A,
- Vector &x,
- const Vector &b,
- const Preconditioner &precond) const
+template <class Vector>
+template<class Matrix, class Preconditioner>
+typename Solver<Vector>::ReturnState
+SolverSelector<Vector>::solve(const Matrix &A,
+ Vector &x,
+ const Vector &b,
+ const Preconditioner &precond) const
{
if (solver_name=="richardson")
{
- SolverRichardson<Matrix,Vector> solver(*control,*vector_memory,
+ SolverRichardson<Vector> solver(*control,*vector_memory,
richardson_data);
return solver.solve(A,x,b,precond);
}
else if (solver_name=="cg")
{
- SolverCG<Matrix,Vector> solver(*control,*vector_memory,
+ SolverCG<Vector> solver(*control,*vector_memory,
cg_data);
return solver.solve(A,x,b,precond);
}
else if (solver_name=="bicgstab")
{
- SolverBicgstab<Matrix,Vector> solver(*control,*vector_memory,
+ SolverBicgstab<Vector> solver(*control,*vector_memory,
bicgstab_data);
return solver.solve(A,x,b,precond);
}
else if (solver_name=="gmres")
{
- SolverGMRES<Matrix,Vector> solver(*control,*vector_memory,
+ SolverGMRES<Vector> solver(*control,*vector_memory,
gmres_data);
return solver.solve(A,x,b,precond);
}
else
Assert(false,ExcSolverDoesNotExist(solver_name));
- return Solver<Matrix,Vector>::breakdown;
+ return Solver<Vector>::breakdown;
};
-template <class Matrix, class Vector>
-string SolverSelector<Matrix, Vector>::get_solver_names()
+template <class Vector>
+string SolverSelector<Vector>::get_solver_names()
{
return "richardson|cg|bicgstab|gmres";
};
-template <class Matrix, class Vector>
-void SolverSelector<Matrix, Vector>::set_data(
- const typename SolverGMRES<Matrix,Vector>::AdditionalData &data)
+template <class Vector>
+void SolverSelector<Vector>::set_data(
+ const typename SolverGMRES<Vector>::AdditionalData &data)
{
gmres_data=data;
}
-template <class Matrix, class Vector>
-void SolverSelector<Matrix, Vector>::set_data(
- const typename SolverRichardson<Matrix,Vector>::AdditionalData &data)
+template <class Vector>
+void SolverSelector<Vector>::set_data(
+ const typename SolverRichardson<Vector>::AdditionalData &data)
{
richardson_data=data;
}
-template <class Matrix, class Vector>
-void SolverSelector<Matrix, Vector>::set_data(
- const typename SolverCG<Matrix,Vector>::AdditionalData &data)
+template <class Vector>
+void SolverSelector<Vector>::set_data(
+ const typename SolverCG<Vector>::AdditionalData &data)
{
cg_data=data;
}
-template <class Matrix, class Vector>
-void SolverSelector<Matrix, Vector>::set_data(
- const typename SolverBicgstab<Matrix,Vector>::AdditionalData &data)
+template <class Vector>
+void SolverSelector<Vector>::set_data(
+ const typename SolverBicgstab<Vector>::AdditionalData &data)
{
bicgstab_data=data;
};
PreconditionMG<FDMG, Vector<TYPE> >
precondition(multigrid, smoother, smoother, coarsegrid);
-// SolverRichardson<SparseMatrix<double> , Vector<TYPE> > solver(control, mem);
- SolverCG<SparseMatrix<double> , Vector<TYPE> > solver(control, mem);
+// SolverRichardson<Vector<TYPE> > solver(control, mem);
+ SolverCG<Vector<TYPE> > solver(control, mem);
Vector<TYPE> u(dim);
Vector<TYPE> f(dim);
projects / dealii.git / commitdiff