#ifndef FILE_CG_NEW
#define FILE_CG_NEW

/**************************************************************************/
/* File:   cg.hh                                                          */
/* Author: Joachim Schoeberl                                              */
/* Date:   5. Jul. 96                                                     */
/**************************************************************************/

/**
   Krylov Space Solver
  */ 
class KrylovSpaceSolver : public BaseMatrix
{
protected:
  ///
  const BaseMatrix *a, *c;
  ///
  double prec;
  ///
  int maxsteps;
  ///
  int steps;
  ///
  int initialize;
  ///
  bool stop_absolute;
  ///
  int printrates;
  ///
  int absoluteRes;
public:
  ///
  KrylovSpaceSolver ();
  ///
  KrylovSpaceSolver (const BaseMatrix & aa);
  ///
  KrylovSpaceSolver (const BaseMatrix & aa, const BaseMatrix & ac);
  ///
  void SetMatrix (const BaseMatrix & aa)
  { a = &aa; }
  ///
  void SetPrecond (const BaseMatrix & ac)
  { c = ∾ }
  ///
  void SetMaxSteps (int amaxsteps)
  { maxsteps = amaxsteps; }

  ///
  void SetPrecision (double aprec)
  { prec = aprec; stop_absolute = 0; }
  ///
  void SetAbsolutePrecision (double aprec)
  {  prec = aprec; stop_absolute = 1; }
  ///
  void SetRelativePrecision (double aprec)
  {  prec = aprec; stop_absolute = 0; }


  void SetPrintRates (int pr = 1)
  { printrates = pr; }
  ///
  void SetInitialize (int ai)
  { initialize = ai; }
  ///
  int GetSteps () const
  { return steps; }
  ///
  virtual void Mult (const BaseVector & v, BaseVector & prod) const = 0;
  ///
  virtual BaseVector * CreateVector () const;
};
  

/// The conjugate gradient solver
template <class SCAL>
class CGSolver : public KrylovSpaceSolver
{
public:
  ///
  CGSolver () 
    : KrylovSpaceSolver () { ; }
  ///
  CGSolver (const BaseMatrix & aa)
    : KrylovSpaceSolver (aa) { ; }

  ///
  CGSolver (const BaseMatrix & aa, const BaseMatrix & ac)
    : KrylovSpaceSolver (aa, ac) { ; }

  ///
  virtual void Mult (const BaseVector & v, BaseVector & prod) const;
};
  

/// The quasi-minimal residual (QMR) solver
template <class SCAL>
class QMRSolver : public KrylovSpaceSolver
{
  int status;
  const BaseMatrix * c2;
public:
  ///
  QMRSolver () 
    : KrylovSpaceSolver (), c2(0) { ; }
  ///
  QMRSolver (const BaseMatrix & aa)
    : KrylovSpaceSolver (aa), c2(0) { ; }

  ///
  QMRSolver (const BaseMatrix & aa, const BaseMatrix & ac)
    : KrylovSpaceSolver (aa, ac), c2(0) { ; }

  ///
  virtual void Mult (const BaseVector & v, BaseVector & prod) const;
};
  




/*

// Conjugate residual solver for symmetric, indefinite matrices


class CRSolver : public CGSolver
{
public:
  ///
  CRSolver ();
  ///
  CRSolver (const BaseMatrix & aa);
  ///
  CRSolver (const BaseMatrix & aa, const BaseMatrix & ac);
  ///
  ~CRSolver ();
  ///
  virtual void Mult (const BaseVector & v, BaseVector & prod) const;
};



//  CG Solver for normal problem 


class NormalCGSolver : public CGSolver
{
public:
  ///
  NormalCGSolver ();
  ///
  NormalCGSolver (const BaseMatrix & aa);
  ///
  NormalCGSolver (const BaseMatrix & aa, const BaseMatrix & ac);
  ///
  ~NormalCGSolver ();
  ///
  virtual void Mult (const BaseVector & v, BaseVector & prod) const;
};



/// gradient method for non-symmetric, indefinite matrices
class MinResGradientSolver : public CGSolver
{
public:
  ///
  MinResGradientSolver ();
  ///
  MinResGradientSolver (const BaseMatrix & aa);
  ///
  MinResGradientSolver (const BaseMatrix & aa, const BaseMatrix & ac);
  ///
  ~MinResGradientSolver ();
  ///
  virtual void Mult (const BaseVector & v, BaseVector & prod) const;
};



//  Bi-CG Stab solver for non-symmetric matrices

class BiCGStabSolver : public CGSolver
{
public:
  ///
  BiCGStabSolver ();
  ///
  BiCGStabSolver (const BaseMatrix & aa);
  ///
  BiCGStabSolver (const BaseMatrix & aa, const BaseMatrix & ac);
  ///
  ~BiCGStabSolver ();
  ///
  virtual void Mult (const BaseVector & v, BaseVector & prod) const;
};









//  QMR solver for non-symmetric matrices

class QMRSolver : public CGSolver
{
  int status;
  ///
  const BaseMatrix * c2;
public:
  ///
  QMRSolver ();
  ///
  QMRSolver (const BaseMatrix & aa);
  ///
  QMRSolver (const BaseMatrix & aa, const BaseMatrix & ac);
  ///
  ~QMRSolver ();
  ///
  virtual void Mult (const BaseVector & v, BaseVector & prod) const;
};






//  complex QMR solver for non-symmetric matrices

class ComplexQMRSolver : public CGSolver
{
public:
  ///
  int status;
  ///
  const BaseMatrix * c2;
public:
  ///
  ComplexQMRSolver ();
  ///
  ComplexQMRSolver (const BaseMatrix & aa);
  ///
  ComplexQMRSolver (const BaseMatrix & aa, const BaseMatrix & ac);
  ///
  ~ComplexQMRSolver ();
  ///
  virtual void Mult (const BaseVector & v, BaseVector & prod) const;
};






//   Simple iteration solver

class SimpleIterationSolver : public CGSolver
{
public:
  ///
  SimpleIterationSolver ();
  ///
  SimpleIterationSolver (const BaseMatrix & aa);
  ///
  SimpleIterationSolver (const BaseMatrix & aa, const BaseMatrix & ac);
  ///
  ~SimpleIterationSolver ();
  ///
  virtual void Mult (const BaseVector & v, BaseVector & prod) const;
};
*/

#endif