From f37ee0a6a5140a57cc6195eac2372ed149c030ad Mon Sep 17 00:00:00 2001 From: Guido Kanschat Date: Tue, 25 Apr 2000 19:08:00 +0000 Subject: [PATCH] Solver classes without matrix template git-svn-id: https://svn.dealii.org/trunk@2770 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/lac/include/lac/eigen.h | 102 ++++++++------ deal.II/lac/include/lac/solver.h | 12 +- deal.II/lac/include/lac/solver_bicgstab.h | 147 ++++++++++---------- deal.II/lac/include/lac/solver_cg.h | 80 +++++------ deal.II/lac/include/lac/solver_gmres.h | 43 +++--- deal.II/lac/include/lac/solver_minres.h | 72 +++++----- deal.II/lac/include/lac/solver_qmrs.h | 112 +++++++-------- deal.II/lac/include/lac/solver_richardson.h | 70 +++++----- deal.II/lac/include/lac/solver_selector.h | 87 ++++++------ tests/lac/mg.cc | 4 +- 10 files changed, 372 insertions(+), 357 deletions(-) diff --git a/deal.II/lac/include/lac/eigen.h b/deal.II/lac/include/lac/eigen.h index e12dcdea72..e9489ccae5 100644 --- a/deal.II/lac/include/lac/eigen.h +++ b/deal.II/lac/include/lac/eigen.h @@ -89,9 +89,8 @@ class ShiftedMatrix * * @author Guido Kanschat, 2000 */ -template , - class VECTOR = Vector > -class EigenPower : private Solver +template > +class EigenPower : private Solver { public: /** @@ -139,7 +138,8 @@ class EigenPower : private Solver * corresponding eigenvector, * normalized with respect to the l2-norm. */ - typename Solver::ReturnState + template + typename Solver::ReturnState solve (double &value, const MATRIX &A, VECTOR &x); @@ -158,21 +158,38 @@ class EigenPower : private Solver * * 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 VECTOR = Vector > -class EigenInverse : private Solver +template > +class EigenInverse : private Solver { 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. @@ -200,7 +217,8 @@ class EigenInverse : private Solver * normalized with respect to the * l2-norm. */ - typename Solver::ReturnState + template + typename Solver::ReturnState solve (double &value, const MATRIX &A, VECTOR &x); @@ -267,25 +285,26 @@ ShiftedMatrix::residual (VECTOR& dst, //---------------------------------------------------------------------- -template -EigenPower::EigenPower (SolverControl &cn, - VectorMemory &mem, - const AdditionalData &data): - Solver(cn, mem), +template +EigenPower::EigenPower (SolverControl &cn, + VectorMemory &mem, + const AdditionalData &data): + Solver(cn, mem), additional_data(data) {} -template -EigenPower::~EigenPower () +template +EigenPower::~EigenPower () {} -template -typename Solver::ReturnState -EigenPower::solve (double &value, - const MATRIX &A, - VECTOR &x) +template +template +typename Solver::ReturnState +EigenPower::solve (double &value, + const MATRIX &A, + VECTOR &x) { SolverControl::State conv=SolverControl::iterate; @@ -350,25 +369,26 @@ EigenPower::solve (double &value, //----------------------------------------------------------------------// -template -EigenInverse::EigenInverse (SolverControl &cn, +template +EigenInverse::EigenInverse (SolverControl &cn, VectorMemory &mem, const AdditionalData &data): - Solver(cn, mem), + Solver(cn, mem), additional_data(data) {} -template -EigenInverse::~EigenInverse () +template +EigenInverse::~EigenInverse () {} -template -typename Solver::ReturnState -EigenInverse::solve (double &value, - const MATRIX &A, - VECTOR &x) +template +template +typename Solver::ReturnState +EigenInverse::solve (double &value, + const MATRIX &A, + VECTOR &x) { deallog.push("Wieland"); @@ -380,7 +400,7 @@ EigenInverse::solve (double &value, // Define solver ReductionControl inner_control (100, 1.e-16, 1.e-8, false, false); PreconditionIdentity prec; - SolverQMRS, VECTOR> + SolverCG solver(inner_control, memory); // Next step for recomputing the shift @@ -433,11 +453,15 @@ EigenInverse::solve (double &value, // 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; } diff --git a/deal.II/lac/include/lac/solver.h b/deal.II/lac/include/lac/solver.h index 4c9fc9964c..88156952e9 100644 --- a/deal.II/lac/include/lac/solver.h +++ b/deal.II/lac/include/lac/solver.h @@ -129,8 +129,7 @@ * * @author Wolfgang Bangerth, Guido Kanschat, Ralf Hartmann, 1997, 1998, 1999 */ -template , - class Vector = Vector > +template > class Solver { public: @@ -179,17 +178,18 @@ class Solver /*-------------------------------- Inline functions ------------------------*/ -template +template inline -SolverControl & Solver::control() const +SolverControl & +Solver::control() const { return cntrl; }; -template +template inline -Solver::Solver(SolverControl &cn, VectorMemory &mem) +Solver::Solver(SolverControl &cn, VectorMemory &mem) : cntrl(cn), memory(mem) {}; diff --git a/deal.II/lac/include/lac/solver_bicgstab.h b/deal.II/lac/include/lac/solver_bicgstab.h index 3e6020b3a1..c1838ed788 100644 --- a/deal.II/lac/include/lac/solver_bicgstab.h +++ b/deal.II/lac/include/lac/solver_bicgstab.h @@ -36,8 +36,8 @@ * has a default argument, so you may call it without the additional * parameter. */ -template , class Vector = Vector > -class SolverBicgstab : public Solver +template > +class SolverBicgstab : private Solver { public: /** @@ -52,7 +52,7 @@ class SolverBicgstab : public Solver * Constructor. */ SolverBicgstab (SolverControl &cn, - VectorMemory &mem, + VectorMemory &mem, const AdditionalData &data=AdditionalData()); /** @@ -63,18 +63,19 @@ class SolverBicgstab : public Solver /** * Solve primal problem only. */ - template - typename Solver::ReturnState - solve (const Matrix &A, - Vector &x, - const Vector &b, + template + typename Solver::ReturnState + solve (const MATRIX &A, + VECTOR &x, + const VECTOR &b, const Preconditioner& precondition); protected: /** * Computation of the stopping criterion. */ - virtual double criterion(); + template + double criterion (const MATRIX& A, const VECTOR& x, const VECTOR& b); /** * Interface for derived class. @@ -86,55 +87,50 @@ class SolverBicgstab : public Solver * 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. @@ -171,50 +167,53 @@ class SolverBicgstab : public Solver /** * Everything before the iteration loop. */ - SolverControl::State start(); + template + SolverControl::State start(const MATRIX& A); /** * The iteration loop itself. */ - template - typename Solver::ReturnState - iterate(const Preconditioner& precondition); + template + typename Solver::ReturnState + iterate(const MATRIX& A, const PRECONDITIONER& precondition); }; /*-------------------------Inline functions -------------------------------*/ -template -SolverBicgstab::SolverBicgstab (SolverControl &cn, - VectorMemory &mem, - const AdditionalData &) +template +SolverBicgstab::SolverBicgstab (SolverControl &cn, + VectorMemory &mem, + const AdditionalData &) : - Solver(cn,mem) + Solver(cn,mem) {} -template -SolverBicgstab::~SolverBicgstab () +template +SolverBicgstab::~SolverBicgstab () {} -template +template +template double -SolverBicgstab::criterion() +SolverBicgstab::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 +template SolverControl::State -SolverBicgstab::start() +SolverBicgstab::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); @@ -224,32 +223,33 @@ SolverBicgstab::start() -template +template void -SolverBicgstab::print_vectors(const unsigned int, - const Vector&, - const Vector&, - const Vector&) const +SolverBicgstab::print_vectors(const unsigned int, + const VECTOR&, + const VECTOR&, + const VECTOR&) const {} -template -template -typename Solver::ReturnState -SolverBicgstab::iterate(const Preconditioner& precondition) +template +template +typename Solver::ReturnState +SolverBicgstab::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 { @@ -258,7 +258,7 @@ SolverBicgstab::iterate(const Preconditioner& precondition) 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; @@ -269,13 +269,13 @@ SolverBicgstab::iterate(const Preconditioner& precondition) 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); @@ -284,13 +284,13 @@ SolverBicgstab::iterate(const Preconditioner& precondition) } -template -template -typename Solver::ReturnState -SolverBicgstab::solve(const Matrix &A, - Vector &x, - const Vector &b, - const Preconditioner& precondition) +template +template +typename Solver::ReturnState +SolverBicgstab::solve(const MATRIX &A, + VECTOR &x, + const VECTOR &b, + const PRECONDITIONER& precondition) { deallog.push("Bicgstab"); Vr = memory.alloc(); Vr->reinit(x); @@ -302,7 +302,6 @@ SolverBicgstab::solve(const Matrix &A, Vt = memory.alloc(); Vt->reinit(x); Vv = memory.alloc(); - MA = &A; Vx = &x; Vb = &b; @@ -314,8 +313,8 @@ SolverBicgstab::solve(const Matrix &A, { 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); diff --git a/deal.II/lac/include/lac/solver_cg.h b/deal.II/lac/include/lac/solver_cg.h index 5a30403b45..81c49d474b 100644 --- a/deal.II/lac/include/lac/solver_cg.h +++ b/deal.II/lac/include/lac/solver_cg.h @@ -43,8 +43,8 @@ * * @author Original implementation by G. Kanschat, R. Becker and F.-T. Suttmeier, reworking and documentation by Wolfgang Bangerth */ -template , class Vector = Vector > -class SolverCG : public Solver +template > +class SolverCG : private Solver { public: /** @@ -59,7 +59,7 @@ class SolverCG : public Solver * Constructor. */ SolverCG (SolverControl &cn, - VectorMemory &mem, + VectorMemory &mem, const AdditionalData &data=AdditionalData()); /** @@ -70,11 +70,11 @@ class SolverCG : public Solver /** * Solver method. */ - template - typename Solver::ReturnState - solve (const Matrix &A, - Vector &x, - const Vector &b, + template + typename Solver::ReturnState + solve (const MATRIX &A, + VECTOR &x, + const VECTOR &b, const Preconditioner& precondition); protected: @@ -96,9 +96,9 @@ class SolverCG : public Solver * 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 @@ -106,10 +106,10 @@ class SolverCG : public Solver * 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 @@ -128,48 +128,48 @@ class SolverCG : public Solver /*------------------------- Implementation ----------------------------*/ -template -SolverCG::SolverCG(SolverControl &cn, - VectorMemory &mem, - const AdditionalData &) +template +SolverCG::SolverCG(SolverControl &cn, + VectorMemory &mem, + const AdditionalData &) : - Solver(cn,mem) + Solver(cn,mem) {} -template -SolverCG::~SolverCG () +template +SolverCG::~SolverCG () {} -template +template double -SolverCG::criterion() +SolverCG::criterion() { return sqrt(res2); } -template +template void -SolverCG::print_vectors(const unsigned int, - const Vector&, - const Vector&, - const Vector&) const +SolverCG::print_vectors(const unsigned int, + const VECTOR&, + const VECTOR&, + const VECTOR&) const {} -template -template -typename Solver::ReturnState -SolverCG::solve (const Matrix &A, - Vector &x, - const Vector &b, - const Preconditioner& precondition) +template +template +typename Solver::ReturnState +SolverCG::solve (const MATRIX &A, + VECTOR &x, + const VECTOR &b, + const Preconditioner& precondition) { SolverControl::State conv=SolverControl::iterate; @@ -181,10 +181,10 @@ SolverCG::solve (const Matrix &A, 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. diff --git a/deal.II/lac/include/lac/solver_gmres.h b/deal.II/lac/include/lac/solver_gmres.h index d20e0bdb5f..2dbf439fe8 100644 --- a/deal.II/lac/include/lac/solver_gmres.h +++ b/deal.II/lac/include/lac/solver_gmres.h @@ -71,9 +71,8 @@ * * @author Wolfgang Bangerth */ -template , - class VECTOR = Vector > -class SolverGMRES : public Solver +template > +class SolverGMRES : private Solver { public: /** @@ -109,11 +108,11 @@ class SolverGMRES : public Solver /** * Solver method. */ - template - typename Solver::ReturnState solve (const MATRIX &A, + template + typename Solver::ReturnState solve (const MATRIX &A, VECTOR &x, const VECTOR &b, - const Preconditioner& precondition); + const PRECONDITIONER& precondition); DeclException1 (ExcTooFewTmpVectors, int, @@ -149,11 +148,11 @@ class SolverGMRES : public Solver /* --------------------- Inline and template functions ------------------- */ -template -SolverGMRES::SolverGMRES (SolverControl &cn, +template +SolverGMRES::SolverGMRES (SolverControl &cn, VectorMemory &mem, const AdditionalData &data) : - Solver (cn,mem), + Solver (cn,mem), additional_data(data) { Assert (additional_data.max_n_tmp_vectors >= 10, @@ -161,14 +160,14 @@ SolverGMRES::SolverGMRES (SolverControl &cn, }; -template +template inline void -SolverGMRES::givens_rotation (Vector &h, - Vector &b, - Vector &ci, - Vector &si, - int col) const +SolverGMRES::givens_rotation (Vector &h, + Vector &b, + Vector &ci, + Vector &si, + int col) const { for (int i=0 ; i::givens_rotation (Vector &h, } -template -template -typename Solver::ReturnState -SolverGMRES::solve (const MATRIX& A, +template +template +typename Solver::ReturnState +SolverGMRES::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 @@ -400,9 +399,9 @@ unsigned int dim = 0; }; -template +template double -SolverGMRES::criterion () +SolverGMRES::criterion () { // dummy implementation. this function is // not needed for the present implementation diff --git a/deal.II/lac/include/lac/solver_minres.h b/deal.II/lac/include/lac/solver_minres.h index 0ef4b63d1b..6bf5e26a35 100644 --- a/deal.II/lac/include/lac/solver_minres.h +++ b/deal.II/lac/include/lac/solver_minres.h @@ -47,8 +47,8 @@ * * @author Thomas Richter, 2000 */ -template , class Vector = Vector > -class SolverMinRes : public Solver +template > +class SolverMinRes : private Solver { public: /** @@ -63,7 +63,7 @@ class SolverMinRes : public Solver * Constructor. */ SolverMinRes (SolverControl &cn, - VectorMemory &mem, + VectorMemory &mem, const AdditionalData &data=AdditionalData()); /** @@ -74,11 +74,11 @@ class SolverMinRes : public Solver /** * Solver method. */ - template - typename Solver::ReturnState - solve (const Matrix &A, - Vector &x, - const Vector &b, + template + typename Solver::ReturnState + solve (const MATRIX &A, + VECTOR &x, + const VECTOR &b, const Preconditioner& precondition); /** @@ -103,9 +103,9 @@ class SolverMinRes : public Solver * 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 @@ -113,9 +113,9 @@ class SolverMinRes : public Solver * 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 @@ -134,46 +134,46 @@ class SolverMinRes : public Solver /*------------------------- Implementation ----------------------------*/ -template -SolverMinRes::SolverMinRes (SolverControl &cn, - VectorMemory &mem, - const AdditionalData &) +template +SolverMinRes::SolverMinRes (SolverControl &cn, + VectorMemory &mem, + const AdditionalData &) : - Solver(cn,mem) + Solver(cn,mem) {} -template -SolverMinRes::~SolverMinRes () +template +SolverMinRes::~SolverMinRes () {} -template +template double -SolverMinRes::criterion() +SolverMinRes::criterion() { return res2; }; -template +template void -SolverMinRes::print_vectors(const unsigned int, - const Vector&, - const Vector&, - const Vector&) const +SolverMinRes::print_vectors(const unsigned int, + const VECTOR&, + const VECTOR&, + const VECTOR&) const {} -template -template -typename Solver::ReturnState -SolverMinRes::solve (const Matrix &A, - Vector &x, - const Vector &b, - const Preconditioner& precondition) +template +template +typename Solver::ReturnState +SolverMinRes::solve (const MATRIX &A, + VECTOR &x, + const VECTOR &b, + const Preconditioner& precondition) { SolverControl::State conv=SolverControl::iterate; @@ -192,7 +192,7 @@ SolverMinRes::solve (const Matrix &A, 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; diff --git a/deal.II/lac/include/lac/solver_qmrs.h b/deal.II/lac/include/lac/solver_qmrs.h index e41631fdd7..c6f181d570 100644 --- a/deal.II/lac/include/lac/solver_qmrs.h +++ b/deal.II/lac/include/lac/solver_qmrs.h @@ -48,8 +48,8 @@ * * @author Guido Kanschat, 1999 */ -template , class Vector = Vector > -class SolverQMRS : public Solver +template > +class SolverQMRS : private Solver { public: /** @@ -96,17 +96,17 @@ class SolverQMRS : public Solver * Constructor. */ SolverQMRS (SolverControl &cn, - VectorMemory &mem, + VectorMemory &mem, const AdditionalData &data=AdditionalData()); /** * Solver method. */ - template - typename Solver::ReturnState - solve (const Matrix &A, - Vector &x, - const Vector &b, + template + typename Solver::ReturnState + solve (const MATRIX &A, + VECTOR &x, + const VECTOR &b, const Preconditioner& precondition); /** @@ -119,9 +119,9 @@ class SolverQMRS : public Solver * 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 @@ -135,24 +135,20 @@ class SolverQMRS : public Solver * 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 @@ -172,9 +168,9 @@ class SolverQMRS : public Solver /** * The iteration loop itself. */ - template - typename Solver::ReturnState - iterate(const Preconditioner& precondition); + template + typename Solver::ReturnState + iterate(const MATRIX& A, const Preconditioner& precondition); /** * The current iteration step. */ @@ -185,41 +181,41 @@ class SolverQMRS : public Solver /*------------------------- Implementation ----------------------------*/ -template -SolverQMRS::SolverQMRS(SolverControl &cn, - VectorMemory &mem, - const AdditionalData &data) : - Solver(cn,mem), +template +SolverQMRS::SolverQMRS(SolverControl &cn, + VectorMemory &mem, + const AdditionalData &data) : + Solver(cn,mem), additional_data(data) {}; -template +template double -SolverQMRS::criterion() +SolverQMRS::criterion() { return sqrt(res2); } -template +template void -SolverQMRS::print_vectors(const unsigned int, - const Vector&, - const Vector&, - const Vector&) const +SolverQMRS::print_vectors(const unsigned int, + const VECTOR&, + const VECTOR&, + const VECTOR&) const {} -template -template -typename Solver::ReturnState -SolverQMRS::solve (const Matrix &A, - Vector &x, - const Vector &b, - const Preconditioner& precondition) +template +template +typename Solver::ReturnState +SolverQMRS::solve (const MATRIX &A, + VECTOR &x, + const VECTOR &b, + const Preconditioner& precondition) { deallog.push("QMRS"); @@ -230,7 +226,6 @@ SolverQMRS::solve (const Matrix &A, Vt = memory.alloc(); Vd = memory.alloc(); - MA = &A; Vx = &x; Vb = &b; // resize the vectors, but do not set @@ -249,7 +244,7 @@ SolverQMRS::solve (const Matrix &A, { if (step) deallog << "Restart step " << step << endl; - state = iterate(precondition); + state = iterate(A, precondition); } while (state == breakdown); @@ -268,10 +263,11 @@ SolverQMRS::solve (const Matrix &A, return state; }; -template -template -typename Solver::ReturnState -SolverQMRS::iterate(const Preconditioner& precondition) +template +template +typename Solver::ReturnState +SolverQMRS::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. @@ -282,15 +278,13 @@ SolverQMRS::iterate(const Preconditioner& precondition) 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; diff --git a/deal.II/lac/include/lac/solver_richardson.h b/deal.II/lac/include/lac/solver_richardson.h index e9ae82f514..eb5a6332e9 100644 --- a/deal.II/lac/include/lac/solver_richardson.h +++ b/deal.II/lac/include/lac/solver_richardson.h @@ -37,8 +37,8 @@ * * @author Ralf Hartmann */ -template , class Vector = Vector > -class SolverRichardson : public Solver +template > +class SolverRichardson : private Solver { public: /** @@ -67,16 +67,16 @@ class SolverRichardson : public Solver * Constructor. */ SolverRichardson (SolverControl &cn, - VectorMemory &mem, + VectorMemory &mem, const AdditionalData &data=AdditionalData()); /** * Solve $Ax=b$ for $x$. */ - template - typename Solver::ReturnState solve (const Matrix &A, - Vector &x, - const Vector &b, + template + typename Solver::ReturnState solve (const MATRIX &A, + VECTOR &x, + const VECTOR &b, const Preconditioner& precondition); /** @@ -95,9 +95,9 @@ class SolverRichardson : public Solver * 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: /** @@ -112,8 +112,8 @@ class SolverRichardson : public Solver * of the actual solution process and * deallocated at the end. */ - Vector *Vr; - Vector *Vd; + VECTOR *Vr; + VECTOR *Vd; /** * Damping parameter. @@ -137,27 +137,27 @@ class SolverRichardson : public Solver /*----------------- Implementation of the Richardson Method ------------------*/ -template -SolverRichardson::SolverRichardson(SolverControl &cn, - VectorMemory &mem, - const AdditionalData &data): - Solver (cn,mem), +template +SolverRichardson::SolverRichardson(SolverControl &cn, + VectorMemory &mem, + const AdditionalData &data): + Solver (cn,mem), additional_data(data) {}; -template -template -typename Solver::ReturnState -SolverRichardson::solve (const Matrix &A, - Vector &x, - const Vector &b, - const Preconditioner& precondition) +template +template +typename Solver::ReturnState +SolverRichardson::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"); @@ -188,27 +188,27 @@ SolverRichardson::solve (const Matrix &A, } -template +template void -SolverRichardson::print_vectors(const unsigned int, - const Vector&, - const Vector&, - const Vector&) const +SolverRichardson::print_vectors(const unsigned int, + const VECTOR&, + const VECTOR&, + const VECTOR&) const {} -template +template inline double -SolverRichardson::criterion() +SolverRichardson::criterion() { return res; } -template +template inline void -SolverRichardson::set_omega(double om) +SolverRichardson::set_omega(double om) { additional_data.omega=om; } diff --git a/deal.II/lac/include/lac/solver_selector.h b/deal.II/lac/include/lac/solver_selector.h index aa47cfbfc6..1122e2a704 100644 --- a/deal.II/lac/include/lac/solver_selector.h +++ b/deal.II/lac/include/lac/solver_selector.h @@ -85,8 +85,7 @@ * * @author Ralf Hartmann, 1999 */ -template , - class Vector = Vector > +template > class SolverSelector { public: @@ -112,8 +111,8 @@ class SolverSelector * was specified in the constructor. * */ - template - typename Solver::ReturnState solve(const Matrix &A, + template + typename Solver::ReturnState solve(const Matrix &A, Vector &x, const Vector &b, const Preconditioner &precond) const; @@ -122,28 +121,28 @@ class SolverSelector * Set the additional data. For more * info see the #Solver# class. */ - void set_data(const typename SolverRichardson + void set_data(const typename SolverRichardson ::AdditionalData &data); /** * Set the additional data. For more * info see the #Solver# class. */ - void set_data(const typename SolverCG + void set_data(const typename SolverCG ::AdditionalData &data); /** * Set the additional data. For more * info see the #Solver# class. */ - void set_data(const typename SolverBicgstab + void set_data(const typename SolverBicgstab ::AdditionalData &data); /** * Set the additional data. For more * info see the #Solver# class. */ - void set_data(const typename SolverGMRES + void set_data(const typename SolverGMRES ::AdditionalData &data); /** @@ -183,116 +182,116 @@ class SolverSelector /** * Stores the additional data. */ - typename SolverRichardson::AdditionalData richardson_data; + typename SolverRichardson::AdditionalData richardson_data; /** * Stores the additional data. */ - typename SolverCG::AdditionalData cg_data; + typename SolverCG::AdditionalData cg_data; /** * Stores the additional data. */ - typename SolverBicgstab::AdditionalData bicgstab_data; + typename SolverBicgstab::AdditionalData bicgstab_data; /** * Stores the additional data. */ - typename SolverGMRES::AdditionalData gmres_data; + typename SolverGMRES::AdditionalData gmres_data; }; /* --------------------- Inline and template functions ------------------- */ -template -SolverSelector::SolverSelector(string solver_name, - SolverControl &control, - VectorMemory &vector_memory) : +template +SolverSelector::SolverSelector(string solver_name, + SolverControl &control, + VectorMemory &vector_memory) : solver_name(solver_name), control(&control), vector_memory(&vector_memory) {}; -template -SolverSelector::~SolverSelector() +template +SolverSelector::~SolverSelector() {}; -template -template -typename Solver::ReturnState -SolverSelector::solve(const Matrix &A, - Vector &x, - const Vector &b, - const Preconditioner &precond) const +template +template +typename Solver::ReturnState +SolverSelector::solve(const Matrix &A, + Vector &x, + const Vector &b, + const Preconditioner &precond) const { if (solver_name=="richardson") { - SolverRichardson solver(*control,*vector_memory, + SolverRichardson solver(*control,*vector_memory, richardson_data); return solver.solve(A,x,b,precond); } else if (solver_name=="cg") { - SolverCG solver(*control,*vector_memory, + SolverCG solver(*control,*vector_memory, cg_data); return solver.solve(A,x,b,precond); } else if (solver_name=="bicgstab") { - SolverBicgstab solver(*control,*vector_memory, + SolverBicgstab solver(*control,*vector_memory, bicgstab_data); return solver.solve(A,x,b,precond); } else if (solver_name=="gmres") { - SolverGMRES solver(*control,*vector_memory, + SolverGMRES solver(*control,*vector_memory, gmres_data); return solver.solve(A,x,b,precond); } else Assert(false,ExcSolverDoesNotExist(solver_name)); - return Solver::breakdown; + return Solver::breakdown; }; -template -string SolverSelector::get_solver_names() +template +string SolverSelector::get_solver_names() { return "richardson|cg|bicgstab|gmres"; }; -template -void SolverSelector::set_data( - const typename SolverGMRES::AdditionalData &data) +template +void SolverSelector::set_data( + const typename SolverGMRES::AdditionalData &data) { gmres_data=data; } -template -void SolverSelector::set_data( - const typename SolverRichardson::AdditionalData &data) +template +void SolverSelector::set_data( + const typename SolverRichardson::AdditionalData &data) { richardson_data=data; } -template -void SolverSelector::set_data( - const typename SolverCG::AdditionalData &data) +template +void SolverSelector::set_data( + const typename SolverCG::AdditionalData &data) { cg_data=data; } -template -void SolverSelector::set_data( - const typename SolverBicgstab::AdditionalData &data) +template +void SolverSelector::set_data( + const typename SolverBicgstab::AdditionalData &data) { bicgstab_data=data; }; diff --git a/tests/lac/mg.cc b/tests/lac/mg.cc index 9c8e3541a7..19a864ceb9 100644 --- a/tests/lac/mg.cc +++ b/tests/lac/mg.cc @@ -155,8 +155,8 @@ for (unsigned int level = 0; level <= maxlevel; ++level) PreconditionMG > precondition(multigrid, smoother, smoother, coarsegrid); -// SolverRichardson , Vector > solver(control, mem); - SolverCG , Vector > solver(control, mem); +// SolverRichardson > solver(control, mem); + SolverCG > solver(control, mem); Vector u(dim); Vector f(dim); -- 2.39.5