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

#ifndef SYSMATBL_H
#define SYSMATBL_H

#include "block.h"
#include "sysmat.h"
#include "sysmatma28.h"

class Generation;

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


class MLBlockMatrix : public MLMatrix, public MA28Matrix
{
  protected:

    const int nComp;
    int   dimension, noOfBlocks;
    Bool  decomposed, ILUDecomposed, inverseDiag;

    Vector<BlockNode*> 	    D;
    Vector<NeighbourBlock*> L;
    Vector<NeighbourBlock*> *symU;
    Vector<SBool>  	    dirichlet, mark, smooth;

    Vector<Num*>   	    rp, lp;
    Vector<int>   	    Block;
    Num*		    aux;

    FixedSizeAllocator      dataAlloc;


  protected:

    virtual void invertDiagonal();
    void setBlockPointers(Vector<Num>& lhs, Vector<Num>& rhs);
    void addBlock(int row, int col, Vector<Num*>& Ap);
    void addBlock(int row, int col, Vector<Num>& data);
    void setBlockArray();

  public:

    MLBlockMatrix(symmetryType symmetry, int spaceDim, int dim, int nComp);
    ~MLBlockMatrix();

    virtual MLBlockMatrix* castToMLBlockMatrix() { return this; }

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

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

    virtual Num& operator()(int row, int col);
    virtual Num& Diag(int i);
    virtual Num  Diag(int i) const;

    virtual void extend(int noOfNodes);

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

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

    virtual void setDirichlet(int node) { dirichlet[node] = True; }

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

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

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

    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);

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

    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 void removeDirichletEntries(DirichletBCs& dirichletBCs);

    // --	  operations for nonlinear solvers:

    //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() { } 


  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 fillMA28Vectors();
    virtual void countEntries(int* dim, int* nEntries);

  public:
 
    virtual void print() const;

    friend class Interface;
};
//-------------------------------------------------------------------------


class LocalMLBlockMatrix : public MLBlockMatrix
{

  public:

    LocalMLBlockMatrix(symmetryType symmetry0, int spaceDim0, int dim0, 
		       int nComp0)
    		 : MLBlockMatrix(symmetry0, spaceDim0, dim0, nComp0) { }

    ~LocalMLBlockMatrix() { }

    virtual LocalMLBlockMatrix* castToLocalMLBlockMatrix() { return this; }

    virtual void invertDiagonal();

    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 removeNonSmoothedEntries();
};

#endif