From f1833de3ac64a86a795bed28dc60c68a78a97533 Mon Sep 17 00:00:00 2001 From: guido Date: Fri, 12 Mar 1999 09:28:02 +0000 Subject: [PATCH] minor changes git-svn-id: https://svn.dealii.org/trunk@999 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/lac/include/lac/ivector.h | 7 +- deal.II/lac/include/lac/solver_bicgstab.h | 203 ++++++++++++++-------- deal.II/lac/source/solver_control.cc | 14 +- 3 files changed, 147 insertions(+), 77 deletions(-) diff --git a/deal.II/lac/include/lac/ivector.h b/deal.II/lac/include/lac/ivector.h index ce6ac8aa94..cc2d22ee48 100644 --- a/deal.II/lac/include/lac/ivector.h +++ b/deal.II/lac/include/lac/ivector.h @@ -158,12 +158,7 @@ protected: }; - - - - - - +/*-------------------------Inline functions -------------------------------*/ inline unsigned int iVector::n() const diff --git a/deal.II/lac/include/lac/solver_bicgstab.h b/deal.II/lac/include/lac/solver_bicgstab.h index 5441fd59c4..f8b7724e34 100644 --- a/deal.II/lac/include/lac/solver_bicgstab.h +++ b/deal.II/lac/include/lac/solver_bicgstab.h @@ -12,18 +12,21 @@ * Bicgstab algorithm by van der Vorst. */ template -class SolverBicgstab : public Solver +class SolverBicgstab // + : public Solver { -public: - SolverBicgstab(SolverControl &cn, VectorMemory &mem) : - Solver(cn,mem) {}; + public: + /** + * Constructor. + */ + SolverBicgstab(SolverControl &cn, VectorMemory &mem); - /** - * Solve primal problem only. - */ - virtual ReturnState solve (const Matrix &A, - Vector &x, - const Vector &b); + /** + * Solve primal problem only. + */ + virtual ReturnState solve (const Matrix &A, + Vector &x, + const Vector &b); protected: /** @@ -31,45 +34,106 @@ public: */ virtual double criterion(); - /** - * Auxiliary vectors. - */ - Vector *Vx, *Vr, *Vrbar, *Vp, *Vy, *Vz, *Vs, *Vt, *Vv; - const Vector *Vb; + /** + * Auxiliary vector. + */ + Vector *Vx; + /** + * Auxiliary vector. + */ + Vector *Vr; + /** + * Auxiliary vector. + */ + Vector *Vrbar; + /** + * Auxiliary vector. + */ + Vector *Vp; + /** + * Auxiliary vector. + */ + Vector *Vy; + /** + * Auxiliary vector. + */ + Vector *Vz; + /** + * Auxiliary vector. + */ + Vector *Vs; + /** + * Auxiliary vector. + */ + Vector *Vt; + /** + * Auxiliary vector. + */ + Vector *Vv; + /** + * Right hand side vector. + */ + const Vector *Vb; - /** - * Pointer to the system matrix. - */ - const Matrix *MA; + /** + * Pointer to the system matrix. + */ + const Matrix *MA; - /** - * Auxiliary values. - */ - double alpha, beta, omega, rho, rhobar; + /** + * Auxiliary value. + */ + double alpha; + /** + * Auxiliary value. + */ + double beta; + /** + * Auxiliary value. + */ + double omega; + /** + * Auxiliary value. + */ + double rho; + /** + * Auxiliary value. + */ + double rhobar; - /** - * Current iteration step. - */ - unsigned step; + /** + * Current iteration step. + */ + unsigned int step; - /** - * Residual. - */ - double res; + /** + * Residual. + */ + double res; -private: - /** - * Everything before the iteration loop. - */ - SolverControl::State start(); - - /** - * The iteration loop itself. - */ - ReturnState iterate(); + private: + /** + * Everything before the iteration loop. + */ + SolverControl::State start(); + + /** + * The iteration loop itself. + */ + ReturnState iterate(); }; +/*-------------------------Inline functions -------------------------------*/ + +template +SolverBicgstab::SolverBicgstab(SolverControl &cn, + VectorMemory &mem) + : + Solver(cn,mem) +{} + + template double SolverBicgstab::criterion() { @@ -105,37 +169,38 @@ SolverBicgstab::iterate() Vector& v = *Vv; do - { - rhobar = r*rbar; - beta = rhobar * alpha / (rho * omega); - rho = rhobar; - p.sadd(beta, 1., r, -beta*omega, v); - MA->precondition(y,p); - MA->vmult(v,y); - rhobar = rbar * v; - - if (fabs(rhobar) < 1.e-19) return ReturnState(breakdown); + { + rhobar = r*rbar; + beta = rhobar * alpha / (rho * omega); + rho = rhobar; + p.sadd(beta, 1., r, -beta*omega, v); + MA->precondition(y,p); + MA->vmult(v,y); + rhobar = rbar * v; + + if (fabs(rhobar) < 1.e-19) return ReturnState(breakdown); - alpha = rho/rhobar; - s.equ(1., r, -alpha, v); - MA->precondition(z,s); - MA->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()); - } + alpha = rho/rhobar; + s.equ(1., r, -alpha, v); + MA->precondition(z,s); + MA->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()); + } while (state == SolverControl::iterate); if (state == SolverControl::success) return success; return exceeded; } + template Solver::ReturnState SolverBicgstab::solve(const Matrix &A, - Vector &x, - const Vector &b) + Vector &x, + const Vector &b) { deallog.push("Bicgstab"); Vr = memory.alloc(); Vr->reinit(x); @@ -156,11 +221,11 @@ SolverBicgstab::solve(const Matrix &A, ReturnState state; do - { - deallog << "Go!" << endl; - if (start() == SolverControl::success) break; - state = iterate(); - } + { + deallog << "Go!" << endl; + if (start() == SolverControl::success) break; + state = iterate(); + } while (state == breakdown); memory.free(Vr); diff --git a/deal.II/lac/source/solver_control.cc b/deal.II/lac/source/solver_control.cc index edd8592fe7..a38ebc5ce7 100644 --- a/deal.II/lac/source/solver_control.cc +++ b/deal.II/lac/source/solver_control.cc @@ -29,8 +29,18 @@ SolverControl::check (const unsigned int step, lstep = step; lvalue = check_value; - if (step>=maxsteps) return failure; - if (check_value <= tol) return success; + if (step>=maxsteps) + { + deallog << "Failure step " << step + << " value " << check_value << endl; + return failure; + } + if (check_value <= tol) + { + deallog << "Convergence step " << step + << " value " << check_value << endl; + return success; + } return iterate; }; -- 2.39.5