#ifndef FILE_NGS_SPARSEMATRIX
#define FILE_NGS_SPARSEMATRIX

/* ************************************************************************/
/* File:   sparsematrix.hpp                                               */
/* Author: Joachim Schoeberl                                              */
/* Date:   01. Oct. 94, 15 Jan. 02                                        */
/* ************************************************************************/


template <class TM> class SparseMatrix;

class BaseJacobiPrecond;
template <class TM> class JacobiPrecond;
template <class TM> class JacobiPrecondSymmetric;

class BaseBlockJacobiPrecond;
template <class TM> class BlockJacobiPrecond;
template <class TM> class BlockJacobiPrecondSymmetric;
 
template <class TM> class SparseCholesky;

template <class TM> class PardisoInverse;
template <class TM> class SuperLUInverse;

/** 
  The graph of a sparse matrix.
*/
class MatrixGraph
{
protected:
  /// number of rows
  int size;

  /// non-zero elements
  int nze; 

  /// column numbers
  // int * colnr;
  MoveableMem<int> colnr;

  /// pointer to first in row
  // int * firsti;
  MoveableMem<int> firsti;
  
  /// position of diagonal element
  //   int * diagi;
  MoveableMem<int> diagi;

  /// owner of arrays ?
  bool owner;

public:
  /// matrix of hight as, uniform number of els/row
  MatrixGraph (int as, int max_elsperrow);
  /// arbitrary number of els/row
  MatrixGraph (const ARRAY<int> & elsperrow);
  /// shadow matrix graph
  MatrixGraph (const MatrixGraph & graph);

  virtual ~MatrixGraph ();

  /// eliminate unused columne indices
  void Compress();
  
   /// returns position of Element (i, j), exception for unused
  int GetPosition (int i, int j) const;

  /// returns position of Element (i, j), -1 for unused
  int GetPositionTest (int i, int j) const;

  /// find positions of n sorted elements, overwrite pos, exception for unused
  void GetPositionsSorted (int row, int n, int * pos) const;

  /// returns position of new element
  int CreatePosition (int i, int j);

  int Size() const { return size; }

  int NZE() const { return nze; }

  FlatArray<const int> GetRowIndices(int i) const
  { return FlatArray<const int> (firsti[i+1]-firsti[i], colnr+firsti[i]); }

  int * GetRowIndices (int i) 
    { return colnr+firsti[i]; }

  int First (int i) const { return firsti[i]; }

  ostream & Print (ostream & ost) const;

  virtual void MemoryUsage (ARRAY<MemoryUsageStruct*> & mu) const;
};





/// A virtual base class for all sparse matrices
class BaseSparseMatrix : virtual public BaseMatrix, 
			 public MatrixGraph
{
public:
  BaseSparseMatrix (int as, int max_elsperrow)
    : MatrixGraph (as, max_elsperrow)
  { ; }

  BaseSparseMatrix (const ARRAY<int> & elsperrow)
    : MatrixGraph (elsperrow)
  { ; }

  BaseSparseMatrix (const MatrixGraph & agraph)
    : MatrixGraph (agraph)
  { ; }   

  BaseSparseMatrix (const BaseSparseMatrix & amat)
    : BaseMatrix(), MatrixGraph (amat)
  { ; }   


  BaseSparseMatrix & operator= (double s)
  {
    AsVector() = s;
    return *this;
  }

  BaseSparseMatrix & Add (double s, const BaseSparseMatrix & m2)
  {
    AsVector() += s * m2.AsVector();
    return *this;
  }


  virtual BaseJacobiPrecond * CreateJacobiPrecond (const BitArray * inner = 0) const 
  {
    throw Exception ("BaseSparseMatrix::CreateJacobiPrecond");
  }

  virtual BaseBlockJacobiPrecond * 
  CreateBlockJacobiPrecond (Table<int> & blocks,
			    const BaseVector * constraint = 0) const
  { 
    throw Exception ("BaseSparseMatrix::CreateBlockJacobiPrecond");
  }

  virtual BaseMatrix * 
  InverseMatrix (BitArray * subset = 0) const
  { 
    throw Exception ("BaseSparseMatrix::CreateInverse called");
  }

  virtual BaseMatrix * 
  InverseMatrix (ARRAY<int> * clusters) const
  { 
    throw Exception ("BaseSparseMatrix::CreateInverse called");
  }

  virtual BaseSparseMatrix * Restrict (const SparseMatrix<double> & prol) const
  {
    throw Exception ("BaseSparseMatrix::Restrict");
  }
};




/// A general, sparse matrix
template<class TM>
class SparseMatrix : public BaseSparseMatrix, 
		     public S_BaseMatrix<typename mat_traits<TM>::TSCAL>
{
protected:
  // TM * data;
  MoveableMem<TM> data;
  VFlatVector<typename mat_traits<TM>::TSCAL> asvec;
  TM nul;
public:
  typedef typename mat_traits<TM>::TSCAL TSCAL;
  typedef typename mat_traits<TM>::TV_ROW TVX;
  typedef typename mat_traits<TM>::TV_COL TVY;

  ///
  SparseMatrix (int as, int max_elsperrow);

  SparseMatrix (const ARRAY<int> & elsperrow);

  SparseMatrix (const MatrixGraph & agraph);

  SparseMatrix (const SparseMatrix<TM> & amat);

  virtual ~SparseMatrix ();

  int Height() const { return size; }
  int Width() const { return size; }
  virtual int VHeight() const { return size; }
  virtual int VWidth() const { return size; }

  TM & operator[] (int i)  { return data[i]; }
  const TM & operator[] (int i) const { return data[i]; }

  FlatArray<const TM> GetRow(int i) const
  { 
    return FlatArray<const TM> (firsti[i+1]-firsti[i], data+firsti[i]); 
  }

  TM & operator() (int row, int col)
  {
    return data[GetPosition(row, col)];
  }

  const TM & operator() (int row, int col) const
  {
    int pos = GetPositionTest (row,col);
    if (pos != -1)
      return data[pos];
    else
      return nul;
  }

  virtual BaseVector & AsVector() 
  {
    asvec.AssignMemory (nze*sizeof(TM)/sizeof(TSCAL), (void*)&data[0]);
    return asvec; 
  }

  virtual const BaseVector & AsVector() const
  { 
    const_cast<VFlatVector<TSCAL>&> (asvec).
      AssignMemory (nze*sizeof(TM)/sizeof(TSCAL), (void*)&data[0]);
    return asvec; 
  }

  virtual BaseMatrix * CreateMatrix () const
  {
    return new SparseMatrix (*this);
  }

  ///
  virtual BaseVector * CreateVector () const
  {
    return new VVector<TVY> (size);
  }


  virtual BaseJacobiPrecond * 
  CreateJacobiPrecond (const BitArray * inner) const
  { 
    return new JacobiPrecond<TM> (*this, inner);
  }

  virtual BaseBlockJacobiPrecond * 
  CreateBlockJacobiPrecond (Table<int> & blocks,
			    const BaseVector * constraint = 0) const
  { 
    return new BlockJacobiPrecond<TM> (*this, blocks);
  }

  virtual BaseMatrix * InverseMatrix (BitArray * subset = 0) const;

  virtual BaseMatrix * InverseMatrix (ARRAY<int> * clusters) const;

  virtual BaseSparseMatrix * Restrict (const SparseMatrix<double> & prol) const
  { return 0; }

  ///
  virtual ostream & Print (ostream & ost) const;

  ///
  virtual void MemoryUsage (ARRAY<MemoryUsageStruct*> & mu) const;

  
  ///
  // template <class TVEC>
  TVY RowTimesVector (int row, const FlatVector<TVX> & vec) const
  {
    /*
    TVY sum = TSCAL(0);
    for (int j = firsti[row]; j < firsti[row+1]; j++)
      sum += data[j] * vec(colnr[j]);
    return sum;
    */

    int nj = firsti[row+1] - firsti[row];
    const int * colpi = colnr+firsti[row];
    const TM * datap = data+firsti[row];
    const TVX * vecp = &vec(0);
    TVY sum = TSCAL(0);
    for (int j = 0; j < nj; j++, colpi++, datap++)
      sum += *datap * vecp[*colpi];
    return sum;
  }

  ///
  // template <class TVEC, class TEL>
  // void AddRowTransToVector (int row, TEL el, TVEC & vec) const
  void AddRowTransToVector (int row, TVY el, FlatVector<TVX> & vec) const
  {
    /*
      // cannot be optimized by the compiler, since hvec could alias
      // the matrix --> keyword 'restrict' will hopefully solve this problem
      
      for (int j = firsti[row]; j < firsti[row+1]; j++)
        hvec(colnr[j]) += Trans(data[j]) * el;
    */

    int nj = firsti[row+1] - firsti[row];
    const int * colpi = colnr+firsti[row];
    const TM * datap = data+firsti[row];
    TVX * vecp = &vec(0);
    for (int j = 0; j < nj; j++, colpi++, datap++)
      vecp[*colpi] += Trans(*datap) * el;
  }

  virtual void MultAdd (double s, const BaseVector & x, BaseVector & y) const;
  virtual void MultTransAdd (double s, const BaseVector & x, BaseVector & y) const;
};






/// A symmetric sparse matrix
template<class TM>
class SparseMatrixSymmetric : public SparseMatrix<TM>
{
public:
  typedef typename mat_traits<TM>::TSCAL TSCAL;
  typedef typename mat_traits<TM>::TV_COL TV_COL;
  typedef typename mat_traits<TM>::TV_ROW TV_ROW;
  typedef typename mat_traits<TM>::TV_COL TVY;
  typedef typename mat_traits<TM>::TV_ROW TVX;

  ///
  SparseMatrixSymmetric (int as, int max_elsperrow);
  ///
  SparseMatrixSymmetric (const ARRAY<int> & elsperrow);
  ///
  SparseMatrixSymmetric (const MatrixGraph & agraph);
  ///
  SparseMatrixSymmetric (const SparseMatrixSymmetric & agraph);
  ///
  virtual ~SparseMatrixSymmetric ();

  SparseMatrixSymmetric & operator= (double s)
  {
    this->AsVector() = s;
    return *this;
  }

  virtual BaseMatrix * CreateMatrix () const
  {
    return new SparseMatrixSymmetric (*this);
  }

  virtual BaseJacobiPrecond * CreateJacobiPrecond (const BitArray * inner) const
  { 
    return new JacobiPrecondSymmetric<TM> (*this, inner);
  }

  virtual BaseBlockJacobiPrecond * 
  CreateBlockJacobiPrecond (Table<int> & blocks,
			    const BaseVector * constraint = 0) const
  { 
    if (!constraint)
      return new BlockJacobiPrecondSymmetric<TM> (*this, blocks);
    else
      return new BlockJacobiPrecondSymmetric<TM> 
	(*this, dynamic_cast<const T_BaseVector<TVX>&>(*constraint).FV(),
	 blocks);
  }

  virtual BaseSparseMatrix * Restrict (const SparseMatrix<double> & prol) const;

  ///
  virtual void MultAdd (double s, const BaseVector & x, BaseVector & y) const;

  virtual void MultTransAdd (double s, const BaseVector & x, BaseVector & y) const
  {
    MultAdd (s, x, y);
  }

  // template <class TVEC>
  TV_COL RowTimesVectorNoDiag (int row, const FlatVector<TVX> & vec) const
  {
    int last = this->firsti[row+1];
    if (this->colnr[last-1] == row) last--;

    /*
    TV_COL sum = TSCAL(0);
    for (int j = this->firsti[row]; j < last; j++)
      sum += this->data[j] * vec(this->colnr[j]);
    return sum;
    */

    // int nj = firsti[row+1] - firsti[row];
    const int * colpi = this->colnr+this->firsti[row];
    const TM * datap = this->data+this->firsti[row];
    const TVX * vecp = &vec(0);

    TVY sum = TSCAL(0);
    for (int j = this->firsti[row]; j < last; j++, colpi++, datap++)
      sum += *datap * vecp[*colpi];
    return sum;
  }

  // template <class TVEC, class TEL>
  void AddRowTransToVectorNoDiag (int row, TVY el, FlatVector<TVX> & vec) const
  {
    /*
    int last = firsti[row+1];
    if (colnr[last-1] == row) last--;
    
    for (int j = firsti[row]; j < last; j++)
      vec(colnr[j]) += Trans(data[j]) * el;
    */
    int last = this->firsti[row+1];
    if (this->colnr[last-1] == row) last--;

    int nj = last - this->firsti[row];
    const int * colpi = this->colnr+this->firsti[row];
    const TM * datap = this->data+this->firsti[row];
    TVX * vecp = &vec(0);
    for (int j = 0; j < nj; j++, colpi++, datap++)
      vecp[*colpi] += Trans(*datap) * el;
  }
  
  BaseSparseMatrix & AddMerge (double s, const SparseMatrixSymmetric  & m2);

  virtual BaseMatrix * InverseMatrix (BitArray * subset = 0) const;

  virtual BaseMatrix * InverseMatrix (ARRAY<int> * clusters) const;
};


template<class TM>
class VarBlockSparseMatrix : public BaseSparseMatrix, 
			     public S_BaseMatrix<typename mat_traits<TM>::TSCAL>
{
  ARRAY<int> block2linear;
  ARRAY<int> data_index;
  MoveableMem<TM> data;

public:

  typedef typename mat_traits<TM>::TSCAL TSCAL;
  typedef typename mat_traits<TM>::TV_COL TV_COL;
  typedef typename mat_traits<TM>::TV_ROW TV_ROW;
  typedef typename mat_traits<TM>::TV_COL TVY;
  typedef typename mat_traits<TM>::TV_ROW TVX;



  VarBlockSparseMatrix (ARRAY<int> & elsperrow, 
			ARRAY<int> & ablock2linear, 
			ARRAY<int> & linear2block, 
			const SparseMatrix<TM> & sm);
  static VarBlockSparseMatrix * Create (const SparseMatrix<TM> & sm);
  virtual ~VarBlockSparseMatrix ();
  virtual void MultAdd (double s, const BaseVector & x, BaseVector & y) const;
};




#endif