17 #ifndef __deal2__solver_bicgstab_h
18 #define __deal2__solver_bicgstab_h
21 #include <deal.II/base/config.h>
22 #include <deal.II/base/logstream.h>
23 #include <deal.II/lac/solver.h>
24 #include <deal.II/lac/solver_control.h>
26 #include <deal.II/base/subscriptor.h>
28 DEAL_II_NAMESPACE_OPEN
66 template <
class VECTOR = Vector<
double> >
125 template<
class MATRIX,
class PRECONDITIONER>
130 const PRECONDITIONER &precondition);
136 template <
class MATRIX>
147 const VECTOR &d)
const;
226 template <
class MATRIX>
232 template<
class MATRIX,
class PRECONDITIONER>
234 iterate(
const MATRIX &A,
const PRECONDITIONER &precondition);
243 template<
class VECTOR>
246 const AdditionalData &data)
249 additional_data(data)
254 template<
class VECTOR>
256 const AdditionalData &data)
259 additional_data(data)
264 template<
class VECTOR>
270 template <
class VECTOR>
271 template <
class MATRIX>
284 template <
class VECTOR >
285 template <
class MATRIX>
290 Vr->sadd(-1.,1.,*Vb);
293 return this->control().check(step, res);
298 template<
class VECTOR>
303 const VECTOR &)
const
308 template<
class VECTOR>
309 template<
class MATRIX,
class PRECONDITIONER>
312 const PRECONDITIONER &precondition)
316 alpha = omega = rho = 1.;
319 VECTOR &rbar = *Vrbar;
334 beta = rhobar * alpha / (rho * omega);
342 p.sadd(beta, 1., r, -beta*omega, v);
344 precondition.vmult(y,p);
352 if (std::fabs(alpha) > 1.e10)
362 print_vectors(step, *Vx, r, y);
366 precondition.vmult(z,r);
369 omega = rhobar/(t*t);
370 Vx->add(alpha, y, omega, z);
373 if (additional_data.exact_residual)
374 res = criterion(A, *Vx, *Vb);
378 state = this->control().check(step, res);
379 print_vectors(step, *Vx, r, y);
386 template<
class VECTOR>
387 template<
class MATRIX,
class PRECONDITIONER>
392 const PRECONDITIONER &precondition)
394 deallog.
push(
"Bicgstab");
395 Vr = this->memory.alloc();
397 Vrbar = this->memory.alloc();
398 Vrbar->reinit(x,
true);
399 Vp = this->memory.alloc();
401 Vy = this->memory.alloc();
403 Vz = this->memory.alloc();
405 Vt = this->memory.alloc();
407 Vv = this->memory.alloc();
420 deallog <<
"Restart step " << step << std::endl;
423 state = iterate(A, precondition);
427 this->memory.free(Vr);
428 this->memory.free(Vrbar);
429 this->memory.free(Vp);
430 this->memory.free(Vy);
431 this->memory.free(Vz);
432 this->memory.free(Vt);
433 this->memory.free(Vv);
440 this->control().last_value()));
446 DEAL_II_NAMESPACE_CLOSE
AdditionalData(const bool exact_residual=true, const double breakdown=1.e-10)
void solve(const MATRIX &A, VECTOR &x, const VECTOR &b, const PRECONDITIONER &precondition)
void vmult(VECTOR &u, const VECTOR &v) const
SolverBicgstab(SolverControl &cn, VectorMemory< VECTOR > &mem, const AdditionalData &data=AdditionalData())
#define AssertThrow(cond, exc)
double criterion(const MATRIX &A, const VECTOR &x, const VECTOR &b)
Stop iteration, goal reached.
virtual void print_vectors(const unsigned int step, const VECTOR &x, const VECTOR &r, const VECTOR &d) const
bool iterate(const MATRIX &A, const PRECONDITIONER &precondition)
void push(const std::string &text)
AdditionalData additional_data
SolverControl::State start(const MATRIX &A)
virtual ~SolverBicgstab()