/* $Id: linsystem.h,v 1.1.1.1 1996/10/02 10:35:51 roitzsch Exp $ */

#ifndef LINSYSTEM_H
#define LINSYSTEM_H

#include "general.h"
#include "utils.h"
#include "numerics.h"
#include "kvector.h"
#include "matrix.h"

class SystemMatrix;
class Preconditioner;
class Interface;
class DirichletBCs;
class Statistic;

//-------------------------------------------------------------------------


class LinSystem
{
  public:

    const symmetryType  symmetry;
    const int 		spaceDim;

  public:

    SystemMatrix*	A;
    Vector<Num>   	b;		// right hand side of lin. system
    Vector<Num>   	bSave;		// to store right hand side before
    					//     dirichlet BCs are applied
    Preconditioner* 	preCond;
    const DirichletBCs* dirichletBCs;

    SystemMatrix* 	AA;			// used in EigenSystems
    SystemMatrix* 	BB;			// used in EigenSystems

    Statistic*          statistic;     // Statistische Daten

  public:

    LinSystem(int spaceDim, symmetryType sym, Preconditioner* preCond,
	      const DirichletBCs* dirichletBCs);

    virtual ~LinSystem();
    
    int Dim() const;
    virtual void update(SystemMatrix* ANew);
    virtual void reset();

    virtual void store(const Matrix<Num>& AElem,const  Vector<Num>& bElem, 
		       const Vector<int>& globalNodes);

    virtual void setDirichletBCs(DirichletBCs& dirichletBCs);

    virtual Bool solve(Vector<Num>& x, Real eNorm, Real error, 
		       Real globalPrecision, int step, int maxIter = -1);

    Real normRhs() const;

    void setRhs(Vector<Num>& rhs);
    void compEnergy(Vector<Num>& x, Real* eNorm, Num* fct, Bool print=False);
    void Residual(Vector<Num>& res, Vector<Num>& sol) const;


    virtual void formNewEigenSystem(Num lambda);

    virtual void memSpaceInfo();
    virtual void print();


  protected:

    int 	  	directSolverLimit, level0DirectSolverLimit;

    Real		ENorm, E0;	// 1/2xAx (energy norm of solution)
    Num			Fct, Fct0;	// 1/2xAx-xb

    Bool 		status;		// solution status
    Real		Error;
    Real		globalPrecision; 
    Real		extPrecFactor;	// tol *= extPrecFactor  (in conv. test)
    int 		maxOrthoGMRes;  // max no. of orthogonalized vectors


    enum LinSolverType 	{ cgOmin, cgOdir, crOmin, crOdir, relax, nonLinRelax,
			  gmRes, biCG, biCGStab, cgs, lsqCG, ddCG, 
			  testSolver1, testSolver2 };
    enum ConvTestType   { ccgDd, ccgDB, ci, residual, decayOfResidual,
			  vectorIteration };

    LinSolverType	solver;
    ConvTestType 	convTest, convTest0;

    int infoLinSystem;
    int timeLinSystem, accTimeLinSystem;

    void notImplemented(const char* s);


  protected:

    Bool CG      (Vector<Num>& x, int maxIter);		// 2-term recurrence
    Bool CGOdir  (Vector<Num>& x, int maxIter);		// 3-term recurrence

    Bool CROmin  (Vector<Num>& x, int maxIter);		// conjugate residuals(2)
    Bool CROdir  (Vector<Num>& x, int maxIter);		// conjugate residuals(3)

    Bool BiCG    (Vector<Num>& x, int maxIter);
    Bool CGS     (Vector<Num>& x, int maxIter);
    Bool BiCGStab(Vector<Num>& x, int maxIter);

    Bool LSQCG   (Vector<Num>& x, int maxIter);		// CG on Normal equations
    Bool DdCG    (Vector<Num>& x, int maxIter);

    Bool GMRes   (Vector<Num>& x, int maxRestarts, int maxOrtho);
    Bool Relax   (Vector<Num>& x, int maxIter);

    Bool NonLinRelax(Vector<Num>& x, int maxIter);

    Bool TestSolver1(Vector<Num>& x, int maxIter);
    Bool TestSolver2(Vector<Num>& x, int maxIter);


    Bool GMResInnerIteration(Vector<Num>& x, int maxOrtho, Real bNorm,
			     Vector<Num>& g, Matrix<Num>& R, 
			     Vector<Vector<Num>*> &v, Vector<Num>& aux, 
			     int nRestart, int* kMax);
    void GMResGivensMult(int j, Num* v, Vector<Num>& c, Vector<Num>& s);
    void GMResQMult     (int k, Num* v, Vector<Num>& c, Vector<Num>& s);
    void GMResRSolve(Vector<Num>& lhs, Matrix<Num>& R, Vector<Num>& rhs, int k);


    Bool convergenceTest(int iter, Vector<Num>& x, Vector<Num>& r, 
			 Real bNorm, Real deltaE=0, Real eAe=0);

    Bool residualConvergence(int iter, Vector<Num>& x, Vector<Num>& r, 
			     Real bNorm);
    Bool residualDecay   (int iter, Vector<Num>& x, Vector<Num>& r, Real bNorm);
    Bool ccgDdConvergence(int iter, Vector<Num>& x, Vector<Num>& r, Real bNorm, 
			  Real deltaE, Real eAe);
    Bool ccgDBConvergence(int iter, Vector<Num>& x, Vector<Num>& r, Real bNorm, 
			  Real deltaE, Real eAe);
    Bool ciConvergence   (int iter, Vector<Num>& x, Vector<Num>& r, Real bNorm, 
			  Real deltaE, Real eAe);
    Bool vectorIterConvergence(int iter, Vector<Num>& x, Vector<Num>& r, 
			       Real bNorm);
    Bool nonLinEnergyConvergence(int iter, Real ENormX, Real ENormE); 

    friend class Problem;
    friend class StaticPotential;
    friend class StaticHeatConduction;
    friend class TransientHeatConduction;
    friend class RK1A;
    friend class RK2;
};
//-------------------------------------------------------------------------

/*

class SparseLinSystem : public LinSystem
{
  public:

    SparseLinSystem (int spaceDim, symmetryType sym, 
		     Preconditioner* 	 preCond,
		     const DirichletBCs* dirichletBCs);
};
//-------------------------------------------------------------------------

class SparseEigenSystem : public SparseLinSystem
{

  public:

    SparseEigenSystem (int spaceDim, symmetryType sym, 
		       const Interface*     interface,
		       const DirichletBCs*  dirichletBCs);
    
    SparseEigenSystem:: ~SparseEigenSystem();
    
    virtual void reset();
    virtual void update();
    virtual void updateAssembledSystem();

    virtual void formNewEigenSystem(Num lambda);
};
*/

//-------------------------------------------------------------------------

#endif