From: guido Date: Tue, 8 Sep 1998 15:16:22 +0000 (+0000) Subject: new Bicgstab X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=5e35e9dcae74870c9034c52c993029738c88387d;p=dealii-svn.git new Bicgstab Iterative methods perform reinit on auxiliary vectors PrimitiveVectorMemory does not take a size git-svn-id: https://svn.dealii.org/trunk@590 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/lac/include/lac/solver_bicgstab.h b/deal.II/lac/include/lac/solver_bicgstab.h new file mode 100644 index 0000000000..5441fd59c4 --- /dev/null +++ b/deal.II/lac/include/lac/solver_bicgstab.h @@ -0,0 +1,183 @@ +/*---------------------------- solver_bicgstab.h ---------------------------*/ +/* $Id$ */ +#ifndef __solver_bicgstab_H +#define __solver_bicgstab_H +/*---------------------------- solver_bicgstab.h ---------------------------*/ + +#include +#include +#include + +/** + * Bicgstab algorithm by van der Vorst. + */ +template +class SolverBicgstab : public Solver +{ +public: + SolverBicgstab(SolverControl &cn, VectorMemory &mem) : + Solver(cn,mem) {}; + + /** + * Solve primal problem only. + */ + virtual ReturnState solve (const Matrix &A, + Vector &x, + const Vector &b); + + protected: + /** + * Computation of the stopping criterion. + */ + virtual double criterion(); + + /** + * Auxiliary vectors. + */ + Vector *Vx, *Vr, *Vrbar, *Vp, *Vy, *Vz, *Vs, *Vt, *Vv; + const Vector *Vb; + + /** + * Pointer to the system matrix. + */ + const Matrix *MA; + + /** + * Auxiliary values. + */ + double alpha, beta, omega, rho, rhobar; + + /** + * Current iteration step. + */ + unsigned step; + + /** + * Residual. + */ + double res; + +private: + /** + * Everything before the iteration loop. + */ + SolverControl::State start(); + + /** + * The iteration loop itself. + */ + ReturnState iterate(); + +}; + +template double +SolverBicgstab::criterion() +{ + res = MA->residual(*Vt, *Vx, *Vb); + return res; +} + + +template < class Matrix, class Vector > SolverControl::State +SolverBicgstab::start() +{ + res = MA->residual(*Vr, *Vx, *Vb); + Vp->reinit(*Vx); + Vv->reinit(*Vx); + Vrbar->equ(1.,*Vr); + SolverControl::State state = control().check(step, res); + return state; +} + +template Solver::ReturnState +SolverBicgstab::iterate() +{ + 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; + + 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); + + 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) +{ + deallog.push("Bicgstab"); + Vr = memory.alloc(); Vr->reinit(x); + Vrbar = memory.alloc(); Vrbar->reinit(x); + Vp = memory.alloc(); + Vy = memory.alloc(); Vy->reinit(x); + Vz = memory.alloc(); Vz->reinit(x); + Vs = memory.alloc(); Vs->reinit(x); + Vt = memory.alloc(); Vt->reinit(x); + Vv = memory.alloc(); + + MA = &A; + Vx = &x; + Vb = &b; + + step = 0; + + ReturnState state; + + do + { + deallog << "Go!" << endl; + if (start() == SolverControl::success) break; + state = iterate(); + } + while (state == breakdown); + + memory.free(Vr); + memory.free(Vrbar); + memory.free(Vp); + memory.free(Vy); + memory.free(Vz); + memory.free(Vs); + memory.free(Vt); + memory.free(Vv); + + deallog.pop(); + return state; +} + + +/*---------------------------- solver_bicgstab.h ---------------------------*/ +/* end of #ifndef __solver_bicgstab_H */ +#endif +/*---------------------------- solver_bicgstab.h ---------------------------*/ diff --git a/deal.II/lac/include/lac/solver_cg.h b/deal.II/lac/include/lac/solver_cg.h index 548b5489f4..e7e8cf6b60 100644 --- a/deal.II/lac/include/lac/solver_cg.h +++ b/deal.II/lac/include/lac/solver_cg.h @@ -83,10 +83,10 @@ SolverPCG::solve (const Matrix &A, 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; + Vr = memory.alloc(); Vector& g = *Vr; g.reinit(x); + Vp = memory.alloc(); Vector& h = *Vp; h.reinit(x); + Vz = memory.alloc(); Vector& d = *Vz; d.reinit(x); + VAp = memory.alloc(); Vector& Ad = *VAp; Ad.reinit(x); // Implementation taken from the DEAL // library