#ifndef FILE_JACOBI
#define FILE_JACOBI

/* *************************************************************************/
/* File:   jacobi.hh                                                      */
/* Author: Joachim Schoeberl                                              */
/* Date:   06. Oct. 96                                                    */
/* *************************************************************************/


/**
   Jacobi and Gauss Seidel smoother
   for scalar, block and system matrices
 */

class BaseJacobiPrecond : virtual public BaseMatrix
{
public:
  virtual void GSSmooth (BaseVector & x, const BaseVector & b) const = 0;
  virtual void GSSmoothBack (BaseVector & x, const BaseVector & b) const = 0;
};

/// A Jaboci preconditioner for general sparse matrices
template <class TM>
class JacobiPrecond : virtual public BaseJacobiPrecond,
		      virtual public S_BaseMatrix<typename mat_traits<TM>::TSCAL>
{
protected:
  const SparseMatrix<TM> & mat;
  ///
  const BitArray * inner;
  ///
  int height;
  ///
  TM * invdiag;
public:
  typedef typename mat_traits<TM>::TV_ROW TVX;

  ///
  JacobiPrecond (const SparseMatrix<TM> & amat, 
		 const BitArray * ainner = NULL);

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

  ///
  virtual BaseVector * CreateVector () const;
  ///
  virtual void GSSmooth (BaseVector & x, const BaseVector & b) const;

  ///
  virtual void GSSmoothBack (BaseVector & x, const BaseVector & b) const;

  ///
  virtual void GSSmoothNumbering (BaseVector & x, const BaseVector & b,
				  const ARRAY<int> & numbering, 
				  int forward = 1) const;
};









/// A Jaboci preconditioner for symmetric sparse matrices
template <class TM>
class JacobiPrecondSymmetric : public JacobiPrecond<TM>
{
public:
  typedef typename mat_traits<TM>::TV_ROW TVX;

  ///
  JacobiPrecondSymmetric (const SparseMatrixSymmetric<TM> & amat, 
			  const BitArray * ainner = NULL);
    /*
    : JacobiPrecond<TM> (amat, ainner);
  { 
    ;
  }
    */

  ///
  virtual void GSSmooth (BaseVector & x, const BaseVector & b) const;
  /*
  {
    int i;

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

    const SparseMatrixSymmetric<TM> & smat =
      dynamic_cast<const SparseMatrixSymmetric<TM>&> (mat);

    // x := b - L^t x
    for (i = 0; i < height; i++)
      {
	smat.AddRowTransToVectorNoDiag (i, -fx(i), fx);
	fx(i) = fb(i);
      }
    
    // x := (L+D)^{-1} x
    for (i = 0; i < height; i++)
      {
	TVX hv = fx(i) - smat.RowTimesVectorNoDiag (i, fx);
	fx(i) = invdiag[i] * hv;
      }
  }
  */

  ///
  virtual void GSSmoothBack (BaseVector & x, const BaseVector & b) const;
  /*
  {
    int i;
    FlatVector<TVX> & fx = 
      dynamic_cast<VFlatVector<TVX> &> (x).FV();
    const FlatVector<TVX> & fb = 
      dynamic_cast<const VFlatVector<TVX> &> (b).FV();

    const SparseMatrixSymmetric<TM> & smat =
      dynamic_cast<const SparseMatrixSymmetric<TM>&> (mat);
    
    for (i = height-1; i >= 0; i--)
      {
	fx(i) = fb(i) - smat.RowTimesVectorNoDiag (i, fx);
      }
    
    for (i = height-1; i >= 0; i--)
      {
	TVX val = invdiag[i] * fx(i);
	fx(i) = val;
	smat.AddRowTransToVectorNoDiag (i, -val, fx);
      }	 
  }
  */

  ///
  virtual void GSSmoothNumbering (BaseVector & x, const BaseVector & b,
				  const ARRAY<int> & numbering, 
				  int forward = 1) const;
};







#endif