/*********************************************************************/
/* File:   sparsematrix.cpp                                          */
/* Author: Joachim Schoeberl                                         */
/* Date:   25. Mar. 2000                                             */
/*********************************************************************/

/* 
   bilinear-form and linear-form integrators
*/

#include <la.hpp>
namespace ngla
{
  using namespace ngla;

  MatrixGraph :: MatrixGraph (int as, int max_elsperrow)
  {
    size = as;
    nze = as * max_elsperrow;

    // colnr = new int[as*max_elsperrow+1];
    colnr.Alloc (as*max_elsperrow+1);
    colnr.SetName ("matrix graph");

    // firsti = new int[as+1];
    firsti.Alloc (as+1);
    firsti.SetName ("matrix graph, table 1");

    // diagi = new int[as+1];
    diagi.Alloc (as+1);
    diagi.SetName ("matrix graph, table 2");

    owner = true;
    
    for (int i = 0; i < as*max_elsperrow; i++)
      colnr[i] = -1;
    colnr[as*max_elsperrow] = 0;
    
    for (int i = 0; i < as+1; i++)
      {
	firsti[i] = i*max_elsperrow;
	diagi[i] = i*max_elsperrow;
      }
  }
  
  MatrixGraph :: MatrixGraph (const ARRAY<int> & elsperrow)
  {
    size = elsperrow.Size();
    owner = true;
    
    //    firsti = new int[size+1];
    firsti.Alloc (size+1);
    firsti.SetName ("matrix graph, table 1");

    // diagi = new int[size+1];
    diagi.Alloc (size+1);
    diagi.SetName ("matrix graph, table 2");
    
    nze = 0;
    for (int i = 0; i < size; i++)
      {
	firsti[i] = diagi[i] = nze;
	nze += elsperrow[i];
      }
    firsti[size] = diagi[size] = nze;
    
    // colnr = new int[nze+1];
    colnr.Alloc (nze+1);
    colnr.SetName ("matrix graph");
    
    for (int i = 0; i < nze; i++)
      colnr[i] = -1;
    colnr[nze] = 0;
  }
    
  MatrixGraph :: MatrixGraph (const MatrixGraph & agraph)
  {
    int i;
    MatrixGraph & graph = const_cast<MatrixGraph&> (agraph);
    size = graph.size;
    nze = graph.nze;
    owner = false;

    /*
    //  firsti = graph.firsti;
    firsti.Swap (graph.firsti);
    //     diagi = graph.diagi;
    diagi.Swap (graph.diagi);
    // colnr = graph.colnr;
    colnr.Swap (graph.colnr);
    */
    firsti.Alloc (size+1);
    diagi.Alloc (size);
    colnr.Alloc (nze);
    for (i = 0; i < size+1; i++)
      firsti[i] = graph.firsti[i];
    for (i = 0; i < size; i++)
      diagi[i] = graph.diagi[i];
    for (i = 0; i < nze; i++)
      colnr[i] = graph.colnr[i];
  }
  
  MatrixGraph :: ~MatrixGraph ()
  {
    if (owner)
      {
	// delete colnr;
	// delete firsti;
	// delete diagi;
      }
  }
  
  void MatrixGraph :: Compress()
  {
    cout << "compress" << endl;
  }
  

  /// returns position of Element (i, j), exception for unused
  int MatrixGraph :: GetPosition (int i, int j) const
  {
    /*
      for (int k = firsti[i]; k < firsti[i+1]; k++)
      if (colnr[k] == j) return k;
    */
    
    int first = firsti[i];
    int last = firsti[i+1];
    while (last > first + 5)
      {
	unsigned int mid = (first+last) / 2;
	if (colnr[mid] > j)
	  last = mid;
	else
	  {
	    if (colnr[mid] == j) return mid;
	    first = mid+1;
	  }
      }
    for (int k = first; k < last; k++)
      if (colnr[k] == j) return k;
    
    stringstream err;
    err << "illegal position: " << i << ", " << j << endl;
    throw Exception (err.str());
  }
  
  
  /// returns position of Element (i, j), -1 for unused
  int MatrixGraph :: GetPositionTest (int i, int j) const
  {
    for (int k = firsti[i]; k < firsti[i+1]; k++)
      if (colnr[k] == j) return k;
    
    return -1;
  }
  
  int MatrixGraph :: CreatePosition (int i, int j)
  {
    for (int k = firsti[i]; k < firsti[i+1]; k++)
      {
	if (colnr[k] == -1)
	  {
	    colnr[k] = j;
	    return k;
	  }
	
	if (colnr[k] == j) return k;
	
	if (colnr[k] > j)
	  {
	    if (colnr[firsti[i+1]-1] != -1)
	      throw Exception ("sparse matrix row full 1 !");
	    
	    for (int l = firsti[i+1]-1; l > k; l--)
	      colnr[l] = colnr[l-1];

	    colnr[k] = j;
	    return k;
	  }
      }
    throw Exception ("sparse matrix row full 2 !");
  }
  


  void MatrixGraph :: 
  GetPositionsSorted (int row, int n, int * pos) const
  {
    if (n == 1)
      {
	pos[0] = GetPosition (row, pos[0]);
	return;
      }
    
    int i = 0;
    int posi = pos[i];
    int endk = firsti[row+1];
    for (int k = firsti[row]; k < endk; k++)
      {
	if (colnr[k] == posi)
	  {
	    pos[i] = k;
	    i++;
	    if (i == n) return;
	    posi = pos[i];
	  }
      }

    throw Exception ("GetPositionSorted: not matching");
  }



  
  
  ostream & MatrixGraph :: Print (ostream & ost) const
  {
    for (int i = 0; i < size; i++)
      {
	ost << "Row " << i << ":";
	
	for (int j = firsti[i]; j < firsti[i+1]; j++)
	  ost << " " << colnr[j];
	ost << "\n";
      }
    return ost;
  }


  void MatrixGraph :: MemoryUsage (ARRAY<MemoryUsageStruct*> & mu) const
  {
    mu.Append (new MemoryUsageStruct ("MatrixGraph", (nze+size)*sizeof(int), 1));
  }







  template <class TM>
  SparseMatrix<TM> ::
  SparseMatrix (int as, int max_elsperrow)
    : BaseMatrix(),
      BaseSparseMatrix (as, max_elsperrow), 
      S_BaseMatrix<typename mat_traits<TM>::TSCAL> (),
      /* asvec(nze, data), */ nul(TSCAL(0))
  {
    // data = new TM[nze];
    data.Alloc (nze);
    data.SetName ("sparse matrix");
  }

  template <class TM>
  SparseMatrix<TM> ::
  SparseMatrix (const ARRAY<int> & elsperrow)
    : BaseMatrix(),
      BaseSparseMatrix (elsperrow), 
      S_BaseMatrix<typename mat_traits<TM>::TSCAL> (),
      /* asvec(nze, data), */ nul(TSCAL(0))
  { 
    // data = new TM[nze];
    data.Alloc (nze);
    data.SetName ("sparse matrix");
  }

  template <class TM>
  SparseMatrix<TM> ::
  SparseMatrix (const MatrixGraph & agraph)
    : BaseMatrix(),
      BaseSparseMatrix (agraph), 
      S_BaseMatrix<typename mat_traits<TM>::TSCAL> (),
      /* data(new TM[nze]), */
      /* asvec(nze, data), */ nul(TSCAL(0))
  { 
    // data = new TM[nze];
    data.Alloc (nze);
    data.SetName ("sparse matrix");
  }
  
  template <class TM>
  SparseMatrix<TM> ::
  SparseMatrix (const SparseMatrix<TM> & amat)
    : BaseMatrix(),
      BaseSparseMatrix (amat), 
      S_BaseMatrix<typename mat_traits<TM>::TSCAL> (),
      /* data(new TM[nze]), */
      /* asvec(nze, data), */ nul(TSCAL(0)) 
  { 
    // BaseMoveableMem::Print ();
    // data = new TM[nze];
    data.Alloc (nze);
    data.SetName ("sparse matrix");
    AsVector() = amat.AsVector(); 
  }
  

  template <class TM>
  SparseMatrix<TM> :: ~SparseMatrix ()
  { 
    // delete data; 
    data.Free();
  }





  template <class TM>
  void SparseMatrix<TM> :: 
  MultAdd (double s, const BaseVector & x, BaseVector & y) const
  {
    /*
    FlatVector<TVX> fx = 
      dynamic_cast<const T_BaseVector<TVX> &> (x).FV();
    FlatVector<TVY> fy = 
      dynamic_cast<T_BaseVector<TVY> &> (y).FV();
    */

    FlatVector<TVX> fx (x.Size(), x.Memory());
    FlatVector<TVY> fy (y.Size(), y.Memory());

    for (int i = 0; i < Height(); i++)
      fy(i) += s * RowTimesVector (i, fx);
  }
  
  template <class TM>
  void SparseMatrix<TM> :: 
  MultTransAdd (double s, const BaseVector & x, BaseVector & y) const
  {
    /*
    FlatVector<TVX> fx = 
      dynamic_cast<const T_BaseVector<TVX> &> (x).FV();
    FlatVector<TVY> fy = 
      dynamic_cast<T_BaseVector<TVY> &> (y).FV();
    */

    FlatVector<TVX> fx (x.Size(), x.Memory());
    FlatVector<TVX> fy (y.Size(), y.Memory());
    
    for (int i = 0; i < Height(); i++)
      AddRowTransToVector (i, fx(i), fy);
  }



  template <class TM>
  BaseMatrix * SparseMatrix<TM> :: InverseMatrix (BitArray * subset) const
  {
#ifdef USE_SUPERLU
    return new SuperLUInverse<TM> (*this, subset);
#else
#ifdef USE_PARDISO
    return new PardisoInverse<TM> (*this, subset);
#else
    return new SparseCholesky<TM> (*this, subset);
#endif
#endif
  }

  template <class TM>
  BaseMatrix * SparseMatrix<TM> :: InverseMatrix (ARRAY<int> * clusters) const
  {
#ifdef USE_SUPERLU
    return new SuperLUInverse<TM> (*this, 0, clusters);
#else
#ifdef USE_PARDISO
    return new PardisoInverse<TM> (*this, 0, clusters);
#else
    return new SparseCholesky<TM> (*this, 0, clusters);
#endif
#endif
  }




  ///
  template <class TM>
  ostream & SparseMatrix<TM> :: Print (ostream & ost) const
  {
    ost << "data = " << data << endl;
    // ost << "&asvec = " << &asvec(0) << endl;
    for (int i = 0; i < size; i++)
      {
	ost << "Row " << i << ":";
	
	for (int j = firsti[i]; j < firsti[i+1]; j++)
	  ost << " " << colnr[j] << ", " << data[j];
	ost << endl;
      }
    return ost;
  }


  template <class TM>
  void SparseMatrix<TM> :: MemoryUsage (ARRAY<MemoryUsageStruct*> & mu) const
  {
    mu.Append (new MemoryUsageStruct ("SparseMatrix", nze*sizeof(TM), 1));
    if (owner) MatrixGraph::MemoryUsage (mu);
  }


















    // : BaseSparseMatrix (as, max_elsperrow), 

  template <class TM>
  SparseMatrixSymmetric<TM> ::
  SparseMatrixSymmetric (int as, int max_elsperrow)
    : SparseMatrix<TM> (as, max_elsperrow)
  { ; }
  
  template <class TM>
  SparseMatrixSymmetric<TM> ::
  SparseMatrixSymmetric (const ARRAY<int> & elsperrow)
    //    : BaseSparseMatrix (elsperrow)
    : SparseMatrix<TM> (elsperrow)
  { ; }


  template <class TM>  
  SparseMatrixSymmetric<TM> ::
  SparseMatrixSymmetric (const MatrixGraph & agraph)
    //    : BaseSparseMatrix (agraph)
    : SparseMatrix<TM> (agraph)
  { ; }

  template <class TM>
  SparseMatrixSymmetric<TM> ::
  SparseMatrixSymmetric (const SparseMatrixSymmetric & amat)
    : BaseMatrix(), SparseMatrix<TM> (amat)
    //    : BaseSparseMatrix (amat), data(new TM[nze]), 
    //      asvec(nze, data), nul(TSCAL(0))
  { 
    this->AsVector() = amat.AsVector(); 
    // *this = amat; 
  }
  
  
  template <class TM>
  SparseMatrixSymmetric<TM> :: ~SparseMatrixSymmetric ()
  {
    ; 
  }

  template <class TM>
  void SparseMatrixSymmetric<TM> :: 
  MultAdd (double s, const BaseVector & x, BaseVector & y) const
  {
    const FlatVector<TV_ROW> fx = 
      dynamic_cast<const T_BaseVector<TV_ROW> &> (x).FV();
    FlatVector<TV_COL> fy = 
      dynamic_cast<T_BaseVector<TV_COL> &> (y).FV();


    for (int i = 0; i < this->Height(); i++)
      {
	fy(i) += s * RowTimesVector (i, fx);
	AddRowTransToVectorNoDiag (i, s * fx(i), fy);
      }
  }

  template <class TM>
  BaseSparseMatrix & SparseMatrixSymmetric<TM> :: 
  AddMerge (double s, const SparseMatrixSymmetric  & m2)
  {
    for (int i = 0; i < m2.Height(); i++)
      for (int j = 0; j < m2.GetRowIndices(i).Size(); j++)
	(*this)(i, m2.GetRowIndices(i)[j]) += s * m2(i, m2.GetRowIndices(i)[j]);
    return *this;
  }


  template <class TM>
  BaseSparseMatrix * 
  SparseMatrixSymmetric<TM> :: Restrict (const SparseMatrix<double> & prol) const
  {
    int n = this->Height();
    IntTable cols(n);
    IntTable hcols(n);
    ARRAY<int> marks(n);
    /*
    cout << "1" << flush;
    int i, j, k, l;
    for (i = 0; i < n; i++)
      for (j = 0; j < this->GetRowIndices(i).Size(); j++)
	{
	  int jj = this->GetRowIndices(i)[j];

	  for (k = 0; k < prol.GetRowIndices(i).Size(); k++)
	    for (l = 0; l < prol.GetRowIndices(jj).Size(); l++)
	      {
		int kk = prol.GetRowIndices(i)[k];
		int ll = prol.GetRowIndices(jj)[l];

		if (kk >= ll)
		  hcols.AddUnique (kk, ll);
		else
		  hcols.AddUnique (ll, kk);
	      }
	}
    */


    cout << "1" << flush;
    marks = -1;
    int i, j, k, l;
    int mcnt = -1;
    for (i = 0; i < n; i++)
      for (k = 0; k < prol.GetRowIndices(i).Size(); k++)
	{
	  int kk = prol.GetRowIndices(i)[k];
	  
	  mcnt++;
	  for (j = 0; j < hcols[kk].Size(); j++)
	    marks[hcols[kk][j]] = mcnt;

	  for (j = 0; j < this->GetRowIndices(i).Size(); j++)
	    {
	      int jj = this->GetRowIndices(i)[j];
	      
	      for (l = 0; l < prol.GetRowIndices(jj).Size(); l++)
		{
		  int ll = prol.GetRowIndices(jj)[l];

		  if (marks[ll] != mcnt)
		    {
		      hcols.Add (kk, ll);
		      marks[ll] = mcnt;
		    }
		}
	    }
	}




    cout << "." << flush;

    for (i = 0; i < hcols.Size(); i++)
      for (j = 0; j < hcols[i].Size(); j++)
	if (hcols[i][j] < i)
	  cols.AddUnique (i, hcols[i][j]);
	else
	  cols.AddUnique (hcols[i][j], i);


    cout << "2" << flush;

    int nc = 0;
    for (i = 0; i < n; i++)
      if (cols[i].Size()) nc = i;
    nc++;

    int nze = 0;
    ARRAY<int> cnts(nc);
    for (i = 0; i < nc; i++)
      {
	cnts[i] = cols[i].Size();
	nze += cnts[i];
      }
   
    SparseMatrixSymmetric<TM> * cmat = 
      new SparseMatrixSymmetric<TM> (cnts);

    for (i = 0; i < nc; i++)
      for (j = 0; j < cols[i].Size(); j++)
	cmat->CreatePosition (i, cols[i][j]);
    (*cmat) = 0;

    cout << "3" << flush;
    
    for (i = 0; i < n; i++)
      for (j = 0; j < this->GetRowIndices(i).Size(); j++)
	{
	  int jj = this->GetRowIndices(i)[j];
	  
	  for (k = 0; k < prol.GetRowIndices(i).Size(); k++)
	    for (l = 0; l < prol.GetRowIndices(jj).Size(); l++)
	      {
		int kk = prol.GetRowIndices(i)[k];
		int ll = prol.GetRowIndices(jj)[l];

		if (kk >= ll)
		  (*cmat)(kk,ll) += 
		    prol[prol.First(i)+k] * 
		    prol[prol.First(jj)+l] *
		    (*this)[this->firsti[i]+j];

		if (ll >= kk && i != jj)
		  (*cmat)(ll,kk) += 
		    prol[prol.First(jj)+l] *
		    prol[prol.First(i)+k] * 
		    (*this)[this->firsti[i]+j];

		/*
		if (kk >= ll)
		  (*cmat)(kk,ll) += prol[First(i)+k] * prol[First(jj)+l]
		    * (*this)[firsti[i]+j];

		if (ll >= kk && i != jj)
		  (*cmat)(ll,kk) += prol[First(jj)+l] * prol[First(i)+k]
		    * (*this)[firsti[i]+j];
		    */

		/*
		if (kk >= ll)
		  (*cmat)(kk,ll) += (prol(i,kk) * prol(jj,ll)) * (*this)(i,jj);

		if (ll >= kk && i != jj)
		  (*cmat)(ll,kk) += (prol(jj,ll) * prol(i,kk)) * (*this)(i,jj);
		*/
	      }
	}
    cout << "4" << flush << endl;
    return cmat;
  }
  
  

  template <class TM>
  BaseMatrix * SparseMatrixSymmetric<TM> :: InverseMatrix (BitArray * subset) const
  {
#ifdef USE_SUPERLU
    return new SuperLUInverse<TM> (*this, subset, 0, 1);
#else
#ifdef USE_PARDISO
    return new PardisoInverse<TM> (*this, subset, 0, 1);
#else
    return new SparseCholesky<TM> (*this, subset);
#endif
#endif
  }

  template <class TM>
  BaseMatrix * SparseMatrixSymmetric<TM> :: InverseMatrix (ARRAY<int> * clusters) const
  {
#ifdef USE_SUPERLU
    return new SuperLUInverse<TM> (*this, 0, clusters, 1);
#else
#ifdef USE_PARDISO
    return new PardisoInverse<TM> (*this, 0, clusters, 1);
#else
    return new SparseCholesky<TM> (*this, 0, clusters);
#endif
#endif
  }











  template <class TM>
  VarBlockSparseMatrix<TM> ::
  VarBlockSparseMatrix (ARRAY<int> & elsperrow, 
			ARRAY<int> & ablock2linear, 
			ARRAY<int> & linear2block,
			const SparseMatrix<TM> & sm)
    : BaseSparseMatrix (elsperrow), block2linear(ablock2linear), data_index(nze+1)
  {
    cout << "constr" << endl;

    cout << "nze = " << nze << endl;

    
    for (int i = 0; i < block2linear.Size()-1; i++)
      {
	FlatArray<const int> rlin = sm.GetRowIndices(block2linear[i]);
	ARRAY<int> rblock(rlin.Size());
	for (int j = 0; j < rlin.Size(); j++)
	  rblock[j] = linear2block[rlin[j]];
	BubbleSort (rblock.Size(), &rblock[0]);
	(*testout) << "rblock = " << rblock << endl;

	if (rblock.Size() > 0) 
	  CreatePosition (i, rblock[0]);

	for (int j = 1; j < rlin.Size(); j++)
	  if (rblock[j] != rblock[j-1])
	    CreatePosition (i, rblock[j]);
      }

    //    (*testout) << "graph = ";
    //    MatrixGraph::Print (*testout);

    cout << "has graph" << endl;



    int ii = 0, cnt = 0;
    for (int i = 0; i < size; i++)
      for (int j = firsti[i]; j < firsti[i+1]; j++)
	{
	  int col = colnr[j];
	  int h = block2linear[i+1]-block2linear[i];
	  int w = block2linear[col+1]-block2linear[col];

	  data_index[ii] = cnt;
	  ii++;
	  cnt += h*w;
	}

    data_index[ii] = cnt;
    cout << "has data_index" << endl;
    cout << "cnt = " << cnt << endl;
    cout << "sm.NZE = " << sm.NZE() << endl;
    data.Alloc (cnt);
    
    ii = 0;
    for (int i = 0; i < size; i++)
      for (int j = firsti[i]; j < firsti[i+1]; j++)
	{
	  int col = colnr[j];
	  int h = block2linear[i+1]-block2linear[i];
	  int w = block2linear[col+1]-block2linear[col];

	  FlatMatrix<TM> fm(h, w, &data[data_index[ii]]);

	  for (int k = 0; k < h; k++)
	    for (int l = 0; l < w; l++)
	      fm(k,l) = sm(block2linear[i]+k, block2linear[col]+l);

	  ii++;
	}
  }


  
  template <class TM>
  VarBlockSparseMatrix<TM> * VarBlockSparseMatrix<TM> ::
  Create (const SparseMatrix<TM> & sm)
  {
    ARRAY<int> b2l;
    ARRAY<int> l2b(sm.Height());
    ARRAY<int> elsperrow;
    
    b2l.Append (0);
    l2b[0] = 0;
    for (int i = 1; i < sm.Height(); i++)
      {
	FlatArray<const int> rold = sm.GetRowIndices(i-1);
	FlatArray<const int> rnew = sm.GetRowIndices(i);


	if (rold == rnew)
	  {
	    l2b[i] = l2b[i-1];
	  }
	else
	  {
	    l2b[i] = l2b[i-1]+1;
	    b2l.Append (i);
	  }
      }
    b2l.Append (sm.Height());
    //    (*testout) << "b2l = " << b2l << endl;
    //    (*testout) << "l2b = " << l2b << endl;
    
    elsperrow.SetSize (b2l.Size()-1);
    elsperrow = 0;

    for (int i = 0; i < b2l.Size()-1; i++)
      {
	FlatArray<const int> rlin = sm.GetRowIndices(b2l[i]);
	ARRAY<int> rblock(rlin.Size());
	for (int j = 0; j < rlin.Size(); j++)
	  rblock[j] = l2b[rlin[j]];
	BubbleSort (rblock.Size(), &rblock[0]);

	if (rblock.Size() > 0) elsperrow[i] = 1;
	for (int j = 1; j < rlin.Size(); j++)
	  if (rblock[j] != rblock[j-1])
	    elsperrow[i]++;
      }

    return new VarBlockSparseMatrix (elsperrow, b2l, l2b, sm);
  }
  

  
  template <class TM>
  VarBlockSparseMatrix<TM> ::
  ~VarBlockSparseMatrix ()
  {
    ;
  }


  template <class TM>
  void VarBlockSparseMatrix<TM> ::
  MultAdd (double s, const BaseVector & x, BaseVector & y) const
  {
    FlatVector<TVX> fx (x.Size(), x.Memory());
    FlatVector<TVY> fy (y.Size(), y.Memory());
    
    for (int i = 0; i < size; i++)
      {
	FlatVector<TVY> fyi (block2linear[i+1]-block2linear[i], &fy[block2linear[i]]);
	for (int j = firsti[i]; j < firsti[i+1]; j++)
	  {
	    int col = colnr[j];
	    FlatVector<TVX> fxi (block2linear[col+1]-block2linear[col], &fy[block2linear[col]]);
	    FlatMatrix<TM> fm(fyi.Size(), fxi.Size(), const_cast<TM*> (&data[data_index[j]]));
	    
	    fyi += s * (fm * fxi);
	  }
      }
  }
  






 
  template class SparseMatrix<double>;
  template class SparseMatrix<Complex>;
#if MAX_SYS_DIM >= 1
  template class SparseMatrix<Mat<1,1,double> >;
  template class SparseMatrix<Mat<1,1,Complex> >;
#endif
#if MAX_SYS_DIM >= 2
  template class SparseMatrix<Mat<2,2,double> >;
  template class SparseMatrix<Mat<2,2,Complex> >;
#endif
#if MAX_SYS_DIM >= 3
  template class SparseMatrix<Mat<3,3,double> >;
  template class SparseMatrix<Mat<3,3,Complex> >;
#endif
#if MAX_SYS_DIM >= 4
  template class SparseMatrix<Mat<4,4,double> >;
  template class SparseMatrix<Mat<4,4,Complex> >;
#endif
#if MAX_SYS_DIM >= 5
  template class SparseMatrix<Mat<5,5,double> >;
  template class SparseMatrix<Mat<5,5,Complex> >;
#endif
#if MAX_SYS_DIM >= 6
  template class SparseMatrix<Mat<6,6,double> >;
  template class SparseMatrix<Mat<6,6,Complex> >;
#endif
#if MAX_SYS_DIM >= 7
  template class SparseMatrix<Mat<7,7,double> >;
  template class SparseMatrix<Mat<7,7,Complex> >;
#endif
#if MAX_SYS_DIM >= 8
  template class SparseMatrix<Mat<8,8,double> >;
  template class SparseMatrix<Mat<8,8,Complex> >;
#endif




  template class SparseMatrixSymmetric<double>;
  template class SparseMatrixSymmetric<Complex>;
#if MAX_SYS_DIM >= 1
  template class SparseMatrixSymmetric<Mat<1,1,double> >;
  template class SparseMatrixSymmetric<Mat<1,1,Complex> >;
#endif
#if MAX_SYS_DIM >= 2
  template class SparseMatrixSymmetric<Mat<2,2,double> >;
  template class SparseMatrixSymmetric<Mat<2,2,Complex> >;
#endif
#if MAX_SYS_DIM >= 3
  template class SparseMatrixSymmetric<Mat<3,3,double> >;
  template class SparseMatrixSymmetric<Mat<3,3,Complex> >;
#endif
#if MAX_SYS_DIM >= 4
  template class SparseMatrixSymmetric<Mat<4,4,double> >;
  template class SparseMatrixSymmetric<Mat<4,4,Complex> >;
#endif
#if MAX_SYS_DIM >= 5
  template class SparseMatrixSymmetric<Mat<5,5,double> >;
  template class SparseMatrixSymmetric<Mat<5,5,Complex> >;
#endif
#if MAX_SYS_DIM >= 6
  template class SparseMatrixSymmetric<Mat<6,6,double> >;
  template class SparseMatrixSymmetric<Mat<6,6,Complex> >;
#endif
#if MAX_SYS_DIM >= 7
  template class SparseMatrixSymmetric<Mat<7,7,double> >;
  template class SparseMatrixSymmetric<Mat<7,7,Complex> >;
#endif
#if MAX_SYS_DIM >= 8
  template class SparseMatrixSymmetric<Mat<8,8,double> >;
  template class SparseMatrixSymmetric<Mat<8,8,Complex> >;
#endif




  template class VarBlockSparseMatrix<double>;
  template class VarBlockSparseMatrix<Complex>;
}