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

#ifndef SYSMATSP_H
#define SYSMATSP_H

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

class ConnectionPattern;
class DirichletBCs;

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


class SparseMatrix : public SystemMatrix, public MA28Matrix
{
  public:

    int 	 dimension;
    Bool	 decomposed, ILUDecomposed;

    Vector<Num>  D, *L, *U;			// sparse matrix vectors
    Vector<int>  col, end;
    Vector<char> dirichlet;

    Vector<Num>  *symU;				// used in symmetric expansion
    Vector<int>  symCol, symEnd;
 

    SparseMatrix (ConnectionPattern& cPattern, symmetryType symmetry, 
		  int spaceDim);
    SparseMatrix (int spaceDim, SystemMatrix* A, SystemMatrix* B=0, 
		  Num lambda=0.0);
    virtual ~SparseMatrix();

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

    virtual void removeLU();

    virtual void reset();
    virtual int  memSpace(int print=0) const;
    virtual Num& operator() (int i, int k);
    virtual Num& Diag(int i) { return D[i]; }
    virtual Num  Diag(int i) const { return D[i]; }

    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 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 setDirichletBCs (DirichletBCs& dirichletBCs, Vector<Num>& b,
				  Vector<Num>& bS, Num& Fct0, Num& E0); 
    virtual void resetDirichletEntries(Vector<Num>& r) const;

    virtual SparseMatrix*  castToSparseMatrix() { return this; }

    virtual void print() const;
    virtual void printMatLabFormat() const;

  protected:

    virtual void allocate(ConnectionPattern& cPattern);

    void setSymDirichletBCs (DirichletBCs& dirichletBCs, Vector<Num>& b, 
			     Vector<Num>& bSave, Num& Fct0, Num& E0); 
    void setASymDirichletBCs(DirichletBCs& dirichletBCs, Vector<Num>& b,
			     Vector<Num>& bSave, Num& Fct0, Num& E0); 


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


    friend class LinSystem;
};


#endif