#ifndef FILE_SPARSECHOLESKY
#define FILE_SPARSECHOLESKY

/* *************************************************************************/
/* File:   sparsecholesky.hpp                                              */
/* Author: Joachim Schoeberl                                               */
/* Date:   18. Jun. 97                                                     */
/* *************************************************************************/

/*
  sparse cholesky factorization
*/



/**
   A sparse cholesky factorization.
   The unknowns are reordered by the minimum degree
   ordering algorithm
*/
template<class TM>
class SparseCholesky : public BaseMatrix
{
  int height, nze;
  // int *order, *lastinrow, *rowindex;
  MoveableMem<int> order, lastinrow, rowindex;
  MoveableMem<int> rowindexblock, rowindexbs, lastinrowblock;
  ///
  //  TM *lfact;
  MoveableMem<TM> lfact;
  MoveableMem<TM> diag;
  ///
  BitArray * inner;
  ARRAY<int> * cluster;
  ///
  MinimumDegreeOrdering * mdo;
  
public:
  typedef typename mat_traits<TM>::TV_COL TV;
  typedef typename mat_traits<TM>::TV_ROW TVX;

  ///
  SparseCholesky (const SparseMatrix<TM> & a, 
		  BitArray * ainner = NULL,
		  ARRAY<int> * acluster = NULL);
  
  ///
  SparseCholesky (const ARRAY<int> & aorder, 
		  const ARRAY<CliqueEl*> & cliques,
		  const ARRAY<MDOVertex> & vertices);		  
  ///
  ~SparseCholesky ();
  ///
  int VHeight() const { return height; }
  ///
  int VWidth() const { return height; }
  ///
  void Allocate (const ARRAY<int> & aorder, 
		 const ARRAY<CliqueEl*> & cliques,
		 const ARRAY<MDOVertex> & vertices);
  ///
  void Factor (const int * blocknr);
  ///
  void FactorNew (const SparseMatrix<TM> & a);
  ///
  virtual void Mult (const BaseVector & x, BaseVector & y) const;

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

  virtual void MemoryUsage (ARRAY<MemoryUsageStruct*> & mu) const
  {
    mu.Append (new MemoryUsageStruct ("SparseChol", nze*sizeof(TM), 1));
  }


  ///
  void Set (int i, int j, const TM & val);
  ///
  const TM & Get (int i, int j) const;
  ///
  void SetOrig (int i, int j, const TM & val)
  { Set (order[i], order[j], val); }

  virtual BaseVector * CreateVector () const
  {
    return new VVector<TV> (height);
  }
};




#endif