/**************************************************************************/
/* File:   blockjacobi.cpp                                                */
/* Author: Joachim Schoeberl                                              */
/* Date:   14. Aug. 2002                                                  */
/**************************************************************************/


#include <la.hpp>

namespace ngla
{
  using namespace ngla;


  
  template <class TM>
  JacobiPrecond<TM> ::
  JacobiPrecond (const SparseMatrix<TM> & amat, 
		 const BitArray * ainner)
    : mat(amat), inner(ainner)
  { 
    height = mat.Height();
    invdiag = new TM[height];
    for (int i = 0; i < height; i++)
      ngbla::CalcInverse (mat(i,i), invdiag[i]);
  }

  ///
  template <class TM>
  JacobiPrecond<TM> :: ~JacobiPrecond () 
  {
    delete invdiag; 
  }

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

    for (int i = 0; i < height; i++)
      fy(i) += s * (invdiag[i] * fx(i));
  }


  ///
  template <class TM>
  BaseVector * JacobiPrecond<TM> :: CreateVector () const 
  {
    return mat.CreateVector();
  }


  ///
  template <class TM>
  void JacobiPrecond<TM> ::
  GSSmooth (BaseVector & x, const BaseVector & b) const 
  {
    
    FlatVector<TVX> fx = 
      dynamic_cast<T_BaseVector<TVX> &> (x).FV();
    const FlatVector<TVX> fb = 
      dynamic_cast<const T_BaseVector<TVX> &> (b).FV();

    for (int i = 0; i < height; i++)
      {
	TVX ax = mat.RowTimesVector (i, fx);
	fx(i) += invdiag[i] * (fb(i) - ax);
      }
  }


  ///
  template <class TM>
  void JacobiPrecond<TM> ::
  GSSmoothBack (BaseVector & x, const BaseVector & b) const 
  {
    FlatVector<TVX> fx = 
      dynamic_cast<T_BaseVector<TVX> &> (x).FV();
    const FlatVector<TVX> fb = 
      dynamic_cast<const T_BaseVector<TVX> &> (b).FV();


    for (int i = height-1; i >= 0; i--)
      {
	TVX ax = mat.RowTimesVector (i, fx);
	fx(i) += invdiag[i] * (fb(i) - ax);
      }
  }

  ///
  template <class TM>
  void JacobiPrecond<TM> ::
  GSSmoothNumbering (BaseVector & x, const BaseVector & b,
		     const ARRAY<int> & numbering, 
		     int forward) const
  {
    ;
  }




  template <class TM>
  JacobiPrecondSymmetric<TM> ::
  JacobiPrecondSymmetric (const SparseMatrixSymmetric<TM> & amat, 
			  const BitArray * ainner)
    : JacobiPrecond<TM> (amat, ainner)
  { 
    ;
  }

  ///
  template <class TM>
  void JacobiPrecondSymmetric<TM> ::
  GSSmooth (BaseVector & x, const BaseVector & b) const 
  {
    int i;

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

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

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

  ///
  template <class TM>
  void JacobiPrecondSymmetric<TM> ::
  GSSmoothBack (BaseVector & x, const BaseVector & b) const 
  {
    int i;
    FlatVector<TVX> fx = 
      dynamic_cast<T_BaseVector<TVX> &> (x).FV();
    const FlatVector<TVX> fb = 
      dynamic_cast<const T_BaseVector<TVX> &> (b).FV();

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

  ///
  template <class TM>
  void JacobiPrecondSymmetric<TM> ::
  GSSmoothNumbering (BaseVector & x, const BaseVector & b,
		     const ARRAY<int> & numbering, 
		     int forward) const
  {
    ;
  }



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


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


}