operator/bicgstab.hh

Go to the documentation of this file.

/* Iterative linear solver: BICGSTAB
 */

#ifndef bicgstab_hh
#define bicgstab_hh

#include <string>
#include "basics/typedefs.hh"
#include "toolbox/sequence.hh"
#include "compositions.hh"

namespace concepts {

  // ******************************************************************* CG **

  template<class F, class G = F>
  class BiCGStab : public VecOperator<F> {
  public:
    BiCGStab(Operator<F>& A, Real maxeps, int maxit = 0,
             uint relres = false, bool throwing = true)
      : VecOperator<F>(A.dimY(), A.dimX()), W_(0), A_(A), M_(0)
      , maxeps_(maxeps), maxit_(maxit), eps_(1.0), it_(0)
      , relres_(relres), throwing_(throwing), stag_(false) {}

    BiCGStab(Operator<F>& A, Operator<G>& Minv, Real maxeps, int maxit = 0,
             bool relres = 0, bool throwing = true, Operator<G>* M = 0)
      : VecOperator<F>(A.dimY(), A.dimX()), W_(&Minv), A_(A), M_(M)
      , maxeps_(maxeps), maxit_(maxit), eps_(1.0), it_(0)
      , relres_(relres), throwing_(throwing), stag_(false) {}

//     virtual void operator()(const Function<F,G>& fncY, Function<F,G>& fncX);
//     void operator()(const Vector<F,G>& fncY, Vector<F,G>& fncX);

    uint iterations() const { return it_; }

    Real epsilon() const { return eps_; }
    bool stagnated() const { return stag_; }
  protected:
    std::ostream& info(std::ostream& os) const;
  private:
    Operator<G>*  W_;
    Operator<F>&  A_;
    Operator<G>*  M_;

    Real maxeps_;
    uint maxit_;
    Real eps_;
    concepts::Sequence<Real> epsVec_;
    uint it_;
    bool relres_;
    bool throwing_;
    bool stag_;

    virtual void apply_(const Vector<F>& fncY, Vector<F>& fncX);

    template<class I>
    bool scalarTooSmall_(I& scalar, const std::string name, const uint& it,
                         const Real& eps, const Real& maxeps) const;
    bool converged_(const Vector<F>& fncY, Vector<F>& fncX,
                    const Vector<F>& x, const uint& it, Real& eps,
                    const Real& maxeps, const Real& l2Y);
    bool stagnated_(const Vector<F>& x, const Vector<F>& s_p,
                    F& alpha_omega);
  };

} // namespace concepts

#endif // bicgstab_hh