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

#ifndef SYSMATML_H
#define SYSMATML_H

#include "sysmat.h"
#include "sysmatmlnodes.h"
#include "sysmatma28.h"

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


class MLSparseMatrix : public MLMatrix, public MA28Matrix
{
  protected:

    int   dimension;
    Bool  decomposed, ILUDecomposed;

    Vector<Num>     	    D;			// diagonal
    Vector<NodeNeighbour*>  L;			// lower triangle
    Vector<NodeNeighbour*>* symU;		// used for symmetric expansion
    Vector<SBool> 	    mark, dirichlet, smooth;

  public:

    MLSparseMatrix(symmetryType symmetry, int spaceDim, int dim=1);   
    ~MLSparseMatrix();

    virtual int  Dim()  const     { return dimension; }
    virtual int  SpaceDim() const { return spaceDim; }
    virtual void reset();

    virtual Num& operator() (int row, int col);
    virtual Num& Diag(int i) { return D[i]; }

    virtual void extend(int noOfNodes);

    virtual Bool isMarked(int node) { return Bool(mark[node]); }
    virtual void setMark(int node)  { mark[node] = True; }
    virtual void markNodes();
    virtual void unMarkNodes();

    virtual void Smooth(int node) 	{ smooth[node] = True; }
    virtual void dontSmooth(int node)	{ smooth[node] = False; }
    virtual Bool isSmoothed(int node)	{ return Bool(smooth[node]); }
    virtual void setDirichlet(int node) { dirichlet[node] = True; }

    virtual void resetRow (int node);
    virtual void removeRow(int node);

    virtual void removeNonSmoothedEntries() { }
    virtual void removeLU();

    virtual void GalerkinRestriction(MLMatrix& Ap1, Generation& generation,
				     const Vector<SBool>* critical=0);

    virtual void smoothNode(int no, Vector<Num>& e, const Vector<Num>& r);
    virtual Num  MultRow(int row, const Vector<Num>& e);

    virtual void symmetricExpansion();
    virtual void removeSymmetricExpansion();

    virtual void Mult  (Vector<Num>& lhs, Vector<Num>& rhs);
    virtual void ATMult(Vector<Num>& lhs, Vector<Num>& rhs);
    virtual void DiagMult(Vector<Num>& lhs, Vector<Num>& rhs, Real omega=1.0);
    virtual void DiagDiv (Vector<Num>& lhs, Vector<Num>& rhs, Real omega=1.0);
    virtual void LMult(Vector<Num>& lhs, Vector<Num>& rhs);
    virtual void UMult(Vector<Num>& lhs, Vector<Num>& rhs);
    virtual void F  (Vector<Num>& lhs, Vector<Num>& rhs, Real omega=1.0);
    virtual void FT (Vector<Num>& lhs, Vector<Num>& rhs, Real omega=1.0);
    virtual void Fm1(Vector<Num>& lhs, Vector<Num>& rhs, Real omega=1.0);
    virtual void FmT(Vector<Num>& lhs, Vector<Num>& rhs, Real omega=1.0);

    virtual Real omegaOpt(int mode) const;

    virtual void setDirichletBCs(DirichletBCs& dirichletBCs, Vector<Num>& b,
				 Vector<Num>& bS, Num& Funct, Num& E0);
    virtual void resetDirichletEntries(Vector<Num>& r) const;

    virtual MLSparseMatrix* castToMLSparseMatrix() { return this; }

    virtual int  memSpace  (int print=0) const;
    virtual int  memSpaceLU(int print=0) const;


    // --      routines for complete/incomplete factorization:

  public:

    virtual Bool DirectSolution() const { return decomposed; }

    virtual void Decompose(Bool removeOriginal=True);
    virtual void FBSubst(Vector<Num>& lhs, Vector<Num>& rhs);

    virtual void ILUDecompose(Real dropTol);
    virtual void ILUFBSubst(Vector<Num>& lhs, Vector<Num>& rhs);

  protected:

    virtual void setSymDirichletBCs(DirichletBCs& dirichletBCs, Vector<Num>& b,
				    Vector<Num>& bS, Num& Funct, Num& E0);
    virtual void setASymDirichletBCs(DirichletBCs& dirichletBCs, Vector<Num>& b,
				     Vector<Num>& bS, Num& Funct, Num& E0);

    // virtual void resetMarkedEntries();

    virtual void fillMA28Vectors();
    virtual void countEntries(int* dim, int* nEntries);


  public:

    virtual void print() const;
    friend class MLInt;
};
//-------------------------------------------------------------------------


class LocalMLSparseMatrix : public MLSparseMatrix
{
  public:

    LocalMLSparseMatrix(symmetryType symmetry0, int spaceDim0, int dim=1) 
			    	: MLSparseMatrix(symmetry0, spaceDim0, dim) { }

    virtual LocalMLSparseMatrix* castToLocalMLSparseMatrix() { return this; }

    virtual void DiagMult(Vector<Num>& lhs, Vector<Num>& rhs, Real omega=1.0);
    virtual void DiagDiv (Vector<Num>& lhs, Vector<Num>& rhs, Real omega=1.0);
    virtual void F  (Vector<Num>& lhs, Vector<Num>& rhs, Real omega=1.0);
    virtual void FT (Vector<Num>& lhs, Vector<Num>& rhs, Real omega=1.0);
    virtual void Fm1(Vector<Num>& lhs, Vector<Num>& rhs, Real omega=1.0);
    virtual void FmT(Vector<Num>& lhs, Vector<Num>& rhs, Real omega=1.0);

    virtual void smoothNode(int no, Vector<Num>& e, const Vector<Num>& r);
    virtual void removeNonSmoothedEntries();
};

#endif