#ifndef FILE_BLOCK_JACOBI
#define FILE_BLOCK_JACOBI

/* *************************************************************************/
/* File:   blockjacobi.hpp                                                  */
/* Author: Joachim Schoeberl                                               */
/* Date:   06. Oct. 96                                                     */
/* *************************************************************************/


/**
   Block - Jacobi and Gauss Seidel smoother
 */

class BaseBlockJacobiPrecond : virtual public BaseMatrix
{
protected:
  Table<int> & blocktable;
  int maxbs;
public:
  BaseBlockJacobiPrecond (Table<int> & ablocktable);
  virtual ~BaseBlockJacobiPrecond ();

  virtual void GSSmooth (BaseVector & x, const BaseVector & b,
			 int steps = 1) const = 0;
  virtual void GSSmoothBack (BaseVector & x, const BaseVector & b,
			     int steps = 1) const = 0;

  int Reorder (FlatArray<int> block, const MatrixGraph & graph,
	       LocalHeap & lh);
};

/**
   A block-Jacobi preconditioner.
   The blocks are specified by a table container
*/
template <class TM>
class BlockJacobiPrecond : virtual public BaseBlockJacobiPrecond,
			   virtual public S_BaseMatrix<typename mat_traits<TM>::TSCAL>
{
protected:
  const SparseMatrix<TM> & mat;
  ARRAY<Matrix<TM>*> invdiag;
public:
  typedef typename mat_traits<TM>::TV_ROW TVX;
  typedef typename mat_traits<TM>::TSCAL TSCAL;

  ///
  BlockJacobiPrecond (const SparseMatrix<TM> & amat, 
		      Table<int> & ablocktable);
  ///
  virtual ~BlockJacobiPrecond ();

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


  ///
  virtual void MultAdd (TSCAL s, const BaseVector & x, BaseVector & y) const 
  {
    int i, j;
    const FlatVector<TVX> fx = 
      dynamic_cast<const T_BaseVector<TVX> &> (x).FV();
    FlatVector<TVX> fy = 
      dynamic_cast<T_BaseVector<TVX> &> (y).FV();

    Vector<TVX> hxmax(maxbs);
    Vector<TVX> hymax(maxbs);

    for (i = 0; i < blocktable.Size(); i++)
      {
	int bs = blocktable[i].Size();
	if (!bs) continue;

	FlatVector<TVX> hx(bs, &hxmax(0));
	FlatVector<TVX> hy(bs, &hymax(0));

	for (j = 0; j < bs; j++)
	  hx(j) = fx(blocktable[i][j]);
	
	hy = (*invdiag[i]) * hx;

	for (j = 0; j < bs; j++)
	  fy(blocktable[i][j]) += s * hy(j);
      }
  }


  ///
  virtual void MultTransAdd (TSCAL s, const BaseVector & x, BaseVector & y) const 
  {
    int i, j;
    const FlatVector<TVX> fx = 
      dynamic_cast<const T_BaseVector<TVX> &> (x).FV();
    FlatVector<TVX> fy = 
      dynamic_cast<T_BaseVector<TVX> &> (y).FV();

    Vector<TVX> hxmax(maxbs);
    Vector<TVX> hymax(maxbs);

    for (i = 0; i < blocktable.Size(); i++)
      {
	int bs = blocktable[i].Size();
	if (!bs) continue;

	FlatVector<TVX> hx(bs, &hxmax(0));
	FlatVector<TVX> hy(bs, &hymax(0));

	for (j = 0; j < bs; j++)
	  hx(j) = fx(blocktable[i][j]);
	
	hy = Trans(*invdiag[i]) * hx;

	for (j = 0; j < bs; j++)
	  fy(blocktable[i][j]) += s * hy(j);
      }
  }


  ///
  virtual BaseVector * CreateVector () const 
  {
    return mat.CreateVector();
  }


  ///
  virtual void GSSmooth (BaseVector & x, const BaseVector & b,
			 int steps = 1) const 
  {
    int i, j, k;

    const FlatVector<TVX> fb = 
      dynamic_cast<const T_BaseVector<TVX> &> (b).FV();
    FlatVector<TVX> fx = 
      dynamic_cast<T_BaseVector<TVX> &> (x).FV();

    Vector<TVX> hxmax(maxbs);
    Vector<TVX> hymax(maxbs);
    for (k = 0; k < steps; k++)
      for (i = 0; i < blocktable.Size(); i++)
	{
	  int bs = blocktable[i].Size();
	  if (!bs) continue;
	  
	  FlatVector<TVX> hx(bs, &hxmax(0));
	  FlatVector<TVX> hy(bs, &hymax(0));
	  
	  for (j = 0; j < bs; j++)
	    {
	      int jj = blocktable[i][j];
	      hx(j) = fb(jj) - mat.RowTimesVector (jj, fx);
	    }
	  
	  hy = (*invdiag[i]) * hx;
	  
	  for (j = 0; j < bs; j++)
	    fx(blocktable[i][j]) += hy(j);
	}
  }
  
  
  ///
  virtual void GSSmoothBack (BaseVector & x, const BaseVector & b,
			     int steps = 1) const 
  {
    int i, j, k;

    const FlatVector<TVX> fb = 
      dynamic_cast<const T_BaseVector<TVX> &> (b).FV();
    FlatVector<TVX> fx = 
      dynamic_cast<T_BaseVector<TVX> &> (x).FV();

    Vector<TVX> hxmax(maxbs);
    Vector<TVX> hymax(maxbs);

    for (k = 0; k < steps; k++)
      for (i = blocktable.Size()-1; i >= 0; i--)
	{
	  int bs = blocktable[i].Size();
	  if (!bs) continue;
	  
	  FlatVector<TVX> hx(bs, &hxmax(0));
	  FlatVector<TVX> hy(bs, &hymax(0));
	  
	  for (j = 0; j < bs; j++)
	    {
	      int jj = blocktable[i][j];
	      hx(j) = fb(jj) - mat.RowTimesVector (jj, fx);
	    }

	  hy = (*invdiag[i]) * hx;
	  
	  for (j = 0; j < bs; j++)
	    fx(blocktable[i][j]) += hy(j);
	}  
  }
  
  ///
  virtual void GSSmoothNumbering (BaseVector & x, const BaseVector & b,
				  const ARRAY<int> & numbering, 
				  int forward = 1) const
  {
    ;
  }

};




#ifdef SYMCHOLESKY

///
template <class TM>
class BlockJacobiPrecondSymmetric : 
  virtual public BaseBlockJacobiPrecond,
  virtual public S_BaseMatrix<typename mat_traits<TM>::TSCAL>
{
protected:
  const SparseMatrixSymmetric<TM> & mat;
  ARRAY<CholeskyFactors<TM>*> invdiag;

public:

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

  ///
  BlockJacobiPrecondSymmetric (const SparseMatrixSymmetric<TM> & amat, 
			       Table<int> & ablocktable);

  BlockJacobiPrecondSymmetric (const SparseMatrixSymmetric<TM> & amat, 
			       const FlatVector<TVX> & constraint,
			       Table<int> & ablocktable);

  ///
  virtual ~BlockJacobiPrecondSymmetric ();

  ///
  virtual void MultAdd (TSCAL s, const BaseVector & x, BaseVector & y) const 
  {
    int i, j;

    const FlatVector<TVX> fx = 
      dynamic_cast<const T_BaseVector<TVX> &> (x).FV();
    FlatVector<TVX> fy = 
      dynamic_cast<T_BaseVector<TVX> &> (y).FV();

    Vector<TVX> hxmax(maxbs);
    Vector<TVX> hymax(maxbs);

    for (i = 0; i < blocktable.Size(); i++)
      {
	int bs = blocktable[i].Size();
	if (!bs) continue;

	FlatVector<TVX> hx(bs, &hxmax(0));
	FlatVector<TVX> hy(bs, &hymax(0));

	for (j = 0; j < bs; j++)
	  hx(j) = fx(blocktable[i][j]);
	
	invdiag[i]->Mult (hx, hy);

	for (j = 0; j < bs; j++)
	  fy(blocktable[i][j]) += s * hy(j);
      }
  }

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

  ///
  virtual BaseVector * CreateVector () const 
  {
    return mat.CreateVector();
  }


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


  ///
  virtual void GSSmooth (BaseVector & x, const BaseVector & b,
			 int steps = 1) const 
  {
    const FlatVector<TVX> fb = 
      dynamic_cast<const T_BaseVector<TVX> &> (b).FV();
    FlatVector<TVX> fx = 
      dynamic_cast<T_BaseVector<TVX> &> (x).FV();

    Vector<TVX> fy(fx.Size());

    // y = b - (D L^T) x

    fy = fb;
    for (int j = 0; j < mat.Height(); j++)
      mat.AddRowTransToVector (j, -fx(j), fy);

    for (int k = 1; k <= steps; k++)
      for (int i = 0; i < blocktable.Size(); i++)
	SmoothBlock (i, fx, fb, fy);

  }
  
  
  ///
  virtual void GSSmoothBack (BaseVector & x, const BaseVector & b,
			     int steps = 1) const 
  {
    const FlatVector<TVX> fb = 
      dynamic_cast<const T_BaseVector<TVX> &> (b).FV();
    FlatVector<TVX> fx = 
      dynamic_cast<T_BaseVector<TVX> &> (x).FV();

    Vector<TVX> fy(fx.Size());

    // y = b - (D L^T) x
    fy = fb;
    for (int j = 0; j < mat.Height(); j++)
      mat.AddRowTransToVector (j, -fx(j), fy);

    for (int k = 1; k <= steps; k++)
      for (int i = blocktable.Size()-1; i >= 0; i--)
	SmoothBlock (i, fx, fb, fy);
  }




  void SmoothBlock (int i, 
		    FlatVector<TVX> & x,
		    const FlatVector<TVX> & b,
		    FlatVector<TVX> & y) const
  {
    int j;

    FlatArray<int> row = blocktable[i];

    int bs = row.Size();
    if (bs == 0) return;
    
    ArrayMem<TVX,1000> dimem(bs);
    ArrayMem<TVX,1000> wimem(bs);
    FlatVector<TVX> di (bs, &dimem[0]);
    FlatVector<TVX> wi (bs, &wimem[0]);

    
    // di = P_i (y - L x)
    for (j = 0; j < bs; j++)
      di(j) = y(row[j]) - mat.RowTimesVectorNoDiag (row[j], x);

    invdiag[i]->Mult (di, wi);

    // x += P_i w
    // y -= (D L^t) P_i w
    for (j = 0; j < bs; j++)
      {
	x(row[j]) += wi(j);
	mat.AddRowTransToVector (row[j], -wi(j), y);
      }
  }

  
  ///
  virtual void GSSmoothNumbering (BaseVector & x, const BaseVector & b,
				  const ARRAY<int> & numbering, 
				  int forward = 1) const
  {
    ;
  }
  
  
  virtual void MemoryUsage (ARRAY<MemoryUsageStruct*> & mu) const
  {
    int nels = 0;
    for (int i = 0; i < blocktable.Size(); i++)
      {
	int bs = blocktable[i].Size();
	nels += (bs * (bs+1))/2;
      }
    mu.Append (new MemoryUsageStruct ("BlockJacSym", nels*sizeof(TM), blocktable.Size()));
  }
};


#else


///
template <class TM>
class BlockJacobiPrecondSymmetric : 
  virtual public BaseBlockJacobiPrecond,
  virtual public S_BaseMatrix<typename mat_traits<TM>::TSCAL>
{
protected:
  const SparseMatrixSymmetric<TM> & mat;
  // ARRAY<FlatBandCholeskyFactors<TM>*> invdiag;
  
  MoveableMem<int> blockstart, blocksize, blockbw;
  MoveableMem<TM> data;
  bool lowmem;
public:
  
  typedef typename mat_traits<TM>::TV_ROW TVX;
  typedef typename mat_traits<TM>::TSCAL TSCAL;
  
  ///
  BlockJacobiPrecondSymmetric (const SparseMatrixSymmetric<TM> & amat, 
			       Table<int> & ablocktable);
  ///
  BlockJacobiPrecondSymmetric (const SparseMatrixSymmetric<TM> & amat, 
			       const FlatVector<TVX> & constraint,
			       Table<int> & ablocktable);
  
  ///
  virtual ~BlockJacobiPrecondSymmetric ();
  
  FlatBandCholeskyFactors<TM> InvDiag (int i) const
  {
    return FlatBandCholeskyFactors<TM> (blocksize[i], 
					blockbw[i], 
					const_cast<TM*>(data+blockstart[i]));
  }

  void ComputeBlockFactor (FlatArray<int> block, int bw, FlatBandCholeskyFactors<TM> & inv) const;
  
  ///
  virtual void MultAdd (TSCAL s, const BaseVector & x, BaseVector & y) const 
  {
    int i, j;

    const FlatVector<TVX> fx = 
      dynamic_cast<const T_BaseVector<TVX> &> (x).FV();
    FlatVector<TVX> fy = 
      dynamic_cast<T_BaseVector<TVX> &> (y).FV();

    Vector<TVX> hxmax(maxbs);
    Vector<TVX> hymax(maxbs);

    for (i = 0; i < blocktable.Size(); i++)
      {
	int bs = blocktable[i].Size();
	if (!bs) continue;

	FlatVector<TVX> hx(bs, &hxmax(0));
	FlatVector<TVX> hy(bs, &hymax(0));

	for (j = 0; j < bs; j++)
	  hx(j) = fx(blocktable[i][j]);
	
	// invdiag[i]->Mult (hx, hy);
	InvDiag(i).Mult (hx, hy);

	/*
	TSCAL scal = InnerProduct (hx, hy);
	(*testout) << "i = " << i << endl;
	(*testout) << "blockjac, (hx,hy) = " << scal << endl;
	// if (scal < 0)
	  {
	    InvDiag(i).Print (*testout) << endl;
	  }
	*/

	for (j = 0; j < bs; j++)
	  fy(blocktable[i][j]) += s * hy(j);
      }
  }


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


  ///
  virtual BaseVector * CreateVector () const 
  {
    return mat.CreateVector();
  }


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


  ///
  virtual void GSSmooth (BaseVector & x, const BaseVector & b,
			 int steps = 1) const 
  {
    const FlatVector<TVX> fb = 
      dynamic_cast<const T_BaseVector<TVX> &> (b).FV();
    FlatVector<TVX> fx = 
      dynamic_cast<T_BaseVector<TVX> &> (x).FV();

    Vector<TVX> fy(fx.Size());

    // y = b - (D L^T) x

    fy = fb;
    for (int j = 0; j < mat.Height(); j++)
      mat.AddRowTransToVector (j, -fx(j), fy);

    for (int k = 1; k <= steps; k++)
      for (int i = 0; i < blocktable.Size(); i++)
	SmoothBlock (i, fx, fb, fy);
  }
  
  
  ///
  virtual void GSSmoothBack (BaseVector & x, const BaseVector & b,
			     int steps = 1) const 
  {
    const FlatVector<TVX> fb = 
      dynamic_cast<const T_BaseVector<TVX> &> (b).FV();
    FlatVector<TVX> fx = 
      dynamic_cast<T_BaseVector<TVX> &> (x).FV();

    Vector<TVX> fy(fx.Size());

    // y = b - (D L^T) x
    fy = fb;
    for (int j = 0; j < mat.Height(); j++)
      mat.AddRowTransToVector (j, -fx(j), fy);

    for (int k = 1; k <= steps; k++)
      for (int i = blocktable.Size()-1; i >= 0; i--)
	SmoothBlock (i, fx, fb, fy);
  }




  void SmoothBlock (int i, 
		    FlatVector<TVX> & x,
		    const FlatVector<TVX> & b,
		    FlatVector<TVX> & y) const
  {
    int j;

    FlatArray<int> row = blocktable[i];

    int bs = row.Size();
    if (bs == 0) return;
    
    ArrayMem<TVX,1000> dimem(bs);
    ArrayMem<TVX,1000> wimem(bs);
    FlatVector<TVX> di (bs, &dimem[0]);
    FlatVector<TVX> wi (bs, &wimem[0]);
    
    // di = P_i (y - L x)
    for (j = 0; j < bs; j++)
      di(j) = y(row[j]) - mat.RowTimesVectorNoDiag (row[j], x);

    // invdiag[i]->Mult (di, wi);
    if (!lowmem)
      InvDiag(i).Mult (di, wi);
    else
      {
	int bw = blockbw[i];
	int bs = blocktable[i].Size();
	DynamicMem<TM> mem(bs*bw);
	FlatBandCholeskyFactors<TM> inv(bs, bw, mem.Ptr());
	ComputeBlockFactor (blocktable[i], bw, inv);
	inv.Mult (di, wi);
      }

    // x += P_i w
    // y -= (D L^t) P_i w
    for (j = 0; j < bs; j++)
      {
	x(row[j]) += wi(j);
	mat.AddRowTransToVector (row[j], -wi(j), y);
      }
  }

  
  ///
  virtual void GSSmoothNumbering (BaseVector & x, const BaseVector & b,
				  const ARRAY<int> & numbering, 
				  int forward = 1) const
  {
    ;
  }
  
  /*  
  virtual void MemoryUsage (ARRAY<MemoryUsageStruct*> & mu) const
  {
    // ARRAY<BandCholeskyFactors<TM>*> invdiag;


    int nels = 0;
    for (int i = 0; i < blocktable.Size(); i++)
      if (invdiag[i])
	{
	  int s = invdiag[i]->Size();
	  int bw = invdiag[i]->BandWidth();
	  nels += s*bw - (bw * (bw-1)) / 2;
	  nels += s;
	}
    mu.Append (new MemoryUsageStruct ("BlockJacSym Table", blocktable.NElements()*sizeof(int), 1));
    mu.Append (new MemoryUsageStruct ("BlockJacSym", nels*sizeof(TM), 2*blocktable.Size()));
  }
  */
};

#endif





#endif