operator/gmres.hh

Go to the documentation of this file.

/* Iterative linear solver: GMRes
 */

#ifndef gmres_hh
#define gmres_hh

#include <string>

#include "basics/noConvException.hh"
#include "basics/typedefs.hh"
#include "compositions.hh"

namespace concepts {

  // ***************************************************************** GMRes **

  template<class F>
  class GMRes : public VecOperator<F> {
  public:
    GMRes(Operator<F>& A, Real maxeps, int maxit = 0, uint rs = 0,
          uint relres = 0, Operator<F>* W = 0)
      : VecOperator<F>(A.dimY(), A.dimX()), W_(W), A_(A), maxeps_(maxeps), maxit_(maxit)
      , eps_(1.0), it_(0), rs_(rs ? rs : maxit), relres_(relres) {}
    virtual ~GMRes();

    uint iterations() const { return it_; }

    Real epsilon() const { return eps_; }
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    struct R {
      R* right;
      R* left;
      Vector<F>* v;
      F sin;
      Real cos;
      F* h;

      R(uint szh, R* lnk = 0) : right(0), v(0) {
        if ((left = lnk)) lnk->right = this;
        h = szh ? new F[szh] : 0;
        std::memset(h, 0, szh * sizeof(F));
      }
      ~R() { delete v;  delete[] h;}
    };

    Operator<F>*  W_;
    Operator<F>&  A_;

    Real maxeps_;
    uint maxit_;
    Real eps_;
    uint it_;
    uint rs_;
    uint relres_;

    virtual void apply_(const Vector<F>& fncY, Vector<F>& fncX)
      throw(NoConvergence, MissingFeature);
  };

} // namespace concepts

#endif // gmres_hh