An Object Finite Element Library

Precond.h

Go to the documentation of this file.

/*==============================================================================

                                    O  F  E  L  I

                          Object  Finite  Element  Library

  ==============================================================================

   Copyright (C) 1998 - 2008 Rachid Touzani

   This program is free software; you can redistribute it and/or modify it under
   the terms of the GNU General Public License as published by the Free
   Software Foundation; Version 2 of the License.

   This program is distributed in the hope that it will be useful, but WITHOUT
   ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
   FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
   details.

   You should have received a copy of the GNU General Public License
   along with this program; if not, write to the :

   Free Software Foundation
   Inc., 59 Temple Place - Suite 330
   Boston, MA  02111-1307, USA

  ==============================================================================

                             Preconditioner classes

  ==============================================================================*/

#ifndef __PRECOND_H
#define __PRECOND_H

#include <iostream>
using std::ostream;
using std::endl;

#include <iomanip>
using std::setw;

#include "OFELI_Config.h"
#include "Vect.h"
#include "util.h"
#include "Matrix.h"

namespace OFELI {

template<class T_> class SpMatrix;

#ifndef DOXYGEN_SHOULD_SKIP_THIS
template<class T_,class M_>
int inv_diag(size_t n, const M_ &A, valarray<T_> &diag)
{
   for (int i=1; i<=int(n); i++) {
      if (A(i,i) == T_(0.))
         return i;
      if ((diag[i-1]=T_(1.)/A(i,i)) == T_(0.))
         return -i;
   }
   return 0;
}

template<class T_>
int inv_diag(size_t n, const Matrix<T_> *A, valarray<T_> &diag)
{
   for (size_t i=1; i<=n; i++) {
      T_ d = (*A)(i,i);
      if (d == T_(0.))
         return i;
      diag[i-1] = T_(1)/d;
      if (diag[i-1] == T_(0.))
         return -i;
   }
   return 0;
}
#endif /* DOXYGEN_SHOULD_SKIP_THIS */

template<class T_> class Prec
{

 public:

    Prec() {
       _type = -1;
#ifdef _OFELI_DEBUG
       std::clog << "An instance of class Prec is constructed.\n";
       std::clog << "File: " << __FILE__ << ", Line: " << __LINE__ << endl;
#endif
    }

    Prec(int type) {
       _type = type;
#ifdef _OFELI_DEBUG
       std::clog << "An instance of class Prec is constructed.\n";
       std::clog << "File: " << __FILE__ << ", Line: " << __LINE__ << endl;
#endif
    }

   Prec(const SpMatrix<T_> &A, int type=IDENT_PREC) {
       _type = type;
       setMatrix(A);
#ifdef _OFELI_DEBUG
       std::clog << "An instance of class Prec is constructed.\n";
       std::clog << "File: " << __FILE__ << ", Line: " << __LINE__ << endl;
#endif
}

    Prec(const Matrix<T_> *A, int type=IDENT_PREC) {
       _type = type;
       setMatrix(A);
#ifdef _OFELI_DEBUG
       std::clog << "An instance of class Prec is constructed.\n";
       std::clog << "File: " << __FILE__ << ", Line: " << __LINE__ << endl;
#endif
    }

    ~Prec()
    {
#ifdef _OFELI_DEBUG
       std::clog << "An instance of class Prec is destructed.\n";
       std::clog << "File: " << __FILE__ << ", Line: " << __LINE__ << endl;
#endif
    }

    void setType(int type) { _type = type; }

    void setMatrix(const Matrix<T_> *A)
    {
       if (_type==-1) {
          cerr << "Error in OFELI::Prec" << endl;
          cerr << "In setMatrix(A) : Choose preconditioner before setting matrix" << endl;
          exit(1);
       }
       _a = (SpMatrix<T_> *)(A);
       _size = _a->getSize();
       if (_type==IDENT_PREC) { }
       else if (_type==DIAG_PREC) {
          _pivot.resize(_size);
          int i;
          if ((i = inv_diag()) != 0) {
             cerr << "\n*** Fatal Error in OFELI ***\nIn Prec::setMatrix(A) : ";
             cerr << "Zero detected in element " << i << endl;
             cerr << "Program stops\n";
             exit(1);
          }
       }
       else if (_type==ILU_PREC) {
          size_t i, j, k=1;
          _pivot.resize(_size);
          _id.resize(_size);
          for (i=0; i<_size; i++) {
             for (j=_a->getRowPtr(i+1); j<_a->getRowPtr(i+2); j++,k++) {
                if (_a->getColInd(j)==i+1) {
                   _id[i] = k;
                   if ((*_a)(k)==0) {
                      cerr << "\n*** Fatal Error in OFELI ***\nIn Prec::setMatrix(a) : ";
                      cerr << "Zero detected in element " << i+1 << endl;
                      cerr << "Program stops\n";
                      exit(1);
                   }
                   _pivot[i] = (*_a)(k);
                }
             }
          }
          bool found;
          T_ c = 0;
          for (i=0; i<_size; i++) {
             _pivot[i] = 1./_pivot[i];
             for (j=_id[i]+1; j<=_a->getRowPtr(i+2)-1; j++) {
                found = false;
                size_t l = _a->getColInd(j);
                for (k=_a->getRowPtr(l); k<=_id[l-1]-1; k++) {
                   if (_a->getColInd(k)==i+1) {
                      found = true;
                      c = (*_a)(k);
                   }
                }
                if (found)
                   _pivot[l-1] -= c * _pivot[i] * (*_a)(j);
             }
          }
       }
       else
          ;
    }

    void setMatrix(const SpMatrix<T_> &A)
    {
       if (_type==-1) {
          cerr << "Error in OFELI::Prec" << endl;
          cerr << "In setMatrix(A) : Choose preconditioner before setting matrix" << endl;
          exit(1);
       }
       _a = &A;
       _size = _a->getSize();
       if (_type==IDENT_PREC) {
       }
       else if (_type==DIAG_PREC) {
          _pivot.resize(_size);
          int i;
          if ((i=inv_diag()) != 0) {
             cerr << "\n*** Fatal Error in OFELI ***\nIn Prec::setMatrix(a) : ";
             cerr << "Zero detected in element " << i << endl;
             cerr << "Program stops\n";
             exit(1);
          }
       }
       else if (_type==ILU_PREC) {
          size_t i, j, k=1;
          _pivot.resize(_size);
          _id.resize(_size);
          for (i=0; i<_size; i++) {
             for (j=_a->getRowPtr(i+1); j<_a->getRowPtr(i+2); j++, k++) {
                if (_a->getColInd(j)==i+1) {
                   _id[i] = k;
                   if ((*_a)(k)==T_(0.)) {
                      cerr << "\n*** Fatal Error in OFELI ***\nIn Prec::setMatrix(a) : ";
                      cerr << "Zero detected in element " << i+1 << endl;
                      cerr << "Program stops\n";
                      exit(1);
                   }
                   _pivot[i] = (*_a)(k);
                }
             }
          }
          bool found;
          T_ c = T_(0.);
          for (i=0; i<_size; i++) {
             _pivot[i] = T_(1.)/_pivot[i];
             for (j=_id[i]+1; j<=_a->getRowPtr(i+2)-1; j++) {
                found = false;
                size_t l = _a->getColInd(j);
                for (k=_a->getRowPtr(l); k<=_id[l-1]-1; k++) {
                   if (_a->getColInd(k)==i+1) {
                      found = true;
                      c = (*_a)(k);
                   }
                }
                if (found)
                   _pivot[l-1] -= c * _pivot[i] * (*_a)(j);
             }
          }
       }
       else
          ;
    }

    void Solve(Vect<T_> &x) const
    {
       if (_type==DIAG_PREC) {
          for (size_t i=0; i<_size; i++)
             x[i] *= _pivot[i];
       }
       else if (_type==ILU_PREC) {
          int i, j;
          Vect<T_> z(_size);
          for (i=0; i<int(_size); i++) {
             T_ s = 0;
             for (j=_a->getRowPtr(i+1); j<=int(_id[i])-1; j++)
                s += (*_a)[j-1] * z(_a->getColInd(j));
             z[i] = _pivot[i] * (x[i]-s);
          }
          for (i=_size-1; i>=0; i--) {
             T_ s = 0;
             for (j=_id[i]+1; j<=int(_a->getRowPtr(i+2))-1; j++) {
                s += (*_a)[j-1] * x(_a->getColInd(j));
                x[i] = z[i] - _pivot[i] * s;
             }
          }
       }
       else
          ;
    }

    void Solve(const Vect<T_> &b, Vect<T_> &x) const
    {
       if (_type==IDENT_PREC) {
          x = b;
       }
       else if (_type==DIAG_PREC) {
          for (size_t i=0; i<_size; i++)
             x[i] = b[i] * _pivot[i];
       }
       else if (_type==ILU_PREC) {
          int i, j;
          Vect<T_> z(_size);
          for (i=0; i<int(_size); i++) {
             T_ s = 0;
             for (j=int(_a->getRowPtr(i+1)); j<=int(_id[i])-1; j++)
                s += (*_a)[j-1] * z(_a->getColInd(j));
             z[i] = _pivot[i] * (b[i]-s);
          }
          for (i=int(_size-1); i>=0; i--) {
             T_ s = 0;
             for (j=int(_id[i]+1); j<=int(_a->getRowPtr(i+2))-1; j++)
                s += (*_a)[j-1] * x(_a->getColInd(j));
             x[i] = z[i] - _pivot[i] * s;
          }
       }
       else
          ;
    }

    void TransSolve(Vect<T_> &x) const
    {
       if (_type==DIAG_PREC) {
          Solve(x);
       }
       else if (_type==ILU_PREC) {
          int i, j;
          Vect<T_> z(_size), temp(x);
          T_ tmp;
          for (i=0; i<int(_size); i++) {
             z[i] = temp[i];
             tmp = _pivot[i] * z[i];
             for (j=_id[i]+1; j<=int(_a->getRowPtr(i+2))-1; j++)
                temp(_a->getColInd(j)) -= tmp * (*_a)(j);
          }
          for (i=_size-1; i>=0; i--) {
             x[i] = _pivot[i] * z[i];
             for (j=_a->getRowPtr(i+1); j<=int(_id[i])-1; j++)
                z(_a->getColInd(j)) -= (*_a)(j) * x[i];
          }
       }
       else
          ;
    }


    void TransSolve(const Vect<T_> &b, Vect<T_> &x) const { x = b; TransSolve(x); }

    T_ & getPivot(size_t i) const { return _pivot[i-1]; }

 private:
   valarray<T_>       _pivot;
   valarray<size_t>   _id;
   const SpMatrix<T_> *_a;
   size_t             _size;
   int                _type;

   int inv_diag()
   {
      for (size_t i=1; i<=_a->getNbRows(); i++) {
         if ((*_a)(i,i) == T_(0.))
            return i;
         if ((_pivot[i-1]=T_(1.)/(*_a)(i,i)) == T_(0.))
            return -int(i);
      }
      return 0;
   }
};

//==============================================================================
//                                CLASS Precond
//==============================================================================

template<class T_> class Precond
{

 public :

    Precond() { }

    Precond (const Matrix<T_> *A) { }

    virtual ~Precond() { }

#ifndef DOXYGEN_SHOULD_SKIP_THIS
    virtual void Solve(Vect<T_> &x) const = 0;
    virtual void Solve(const Vect<T_> &b, Vect<T_> &x) const = 0;
    virtual void TransSolve(Vect<T_> &x) const = 0;
    virtual void TransSolve(const Vect<T_> &b, Vect<T_> &x) const = 0;
#endif /* DOXYGEN_SHOULD_SKIP_THIS */

 protected:
#ifndef DOXYGEN_SHOULD_SKIP_THIS
   valarray<T_>     _pivot;
   valarray<size_t> _id;
   const Matrix<T_> *_a;
   size_t _size;
#endif /* DOXYGEN_SHOULD_SKIP_THIS */

};

/*-----------------------------------------------------------------------------
                             Identity Preconditioner
  -----------------------------------------------------------------------------*/

template<class T_,class M_> class IdentPrec
{

 public:

    IdentPrec(const M_ &A);

    ~IdentPrec();

    void Solve(Vect<T_> &x) const;

    void Solve(const Vect<T_> &b, Vect<T_> &x) const;

    void TransSolve(Vect<T_> &x) const { }

    void TransSolve(const Vect<T_> &b, Vect<T_> &x) const;

};


//==============================================================================
//                                CLASS IdentPrecond
//==============================================================================


#ifndef DOXYGEN_SHOULD_SKIP_THIS

template<class T_> class IdentPrecond : virtual public Precond<T_>
{

 public:

    IdentPrecond(const Matrix<T_> *A);

    ~IdentPrecond();

    void Solve(Vect<T_> &x) const { }

    void Solve(const Vect<T_> &b, Vect<T_> &x) const;

    void TransSolve(Vect<T_> &x) const { }

    void TransSolve(const Vect<T_> &b, Vect<T_> &x) const;

};
#endif /* DOXYGEN_SHOULD_SKIP_THIS */


//==============================================================================
//                                CLASS DiagPrec
//==============================================================================

template<class T_, class M_> class DiagPrec
{

 public:

    DiagPrec(const M_ &a)
    {
       _size = a.getSize();
       _pivot.resize(_size);
       _a = &a;
       int i;
       if ((i = inv_diag(_size, *_a, _pivot)) != 0) {
          cerr << "\n*** Fatal Error in OFELI ***\nIn DiagPrec::DiagPrec(M_ &) : ";
          cerr << "Zero detected in element " << i << endl;
          cerr << "Program stops\n";
          exit(1);
       }
#ifdef _OFELI_DEBUG
       std::clog << "An instance of class DiagPrec is constructed.\n";
       std::clog << "File : " << __FILE__ << ", Line : " << __LINE__ << endl;
#endif
    }

    ~DiagPrec()
    {
#ifdef _OFELI_DEBUG
       std::clog << "An instance of class DiagPrec is constructed.\n";
       std::clog << "File : " << __FILE__ << ", Line : " << __LINE__ << endl;
#endif
    }

    void Solve(Vect<T_> &x) const
    {
       for (size_t i=0; i<_size; i++)
          x[i] *= _pivot[i];
    }

    void Solve(const Vect<T_> &b, Vect<T_> &x) const
    {
       for (size_t i=0; i<_size; i++)
          x[i] = b[i] * _pivot[i];
    }

    void TransSolve(Vect<T_> &x) const { Solve(x); }

    void TransSolve(const Vect<T_> &b, Vect<T_> &x) const { Solve(b,x); }

    const T_ & getDiag(size_t i) const { return _pivot[i-1]; }

    T_ & getPivot(size_t i) const { return _pivot[i-1]; }

 private:
   valarray<T_> _pivot;
   const M_     *_a;
   size_t       _size;
};


//==============================================================================
//                                CLASS DiagPrecond
//==============================================================================

#ifndef DOXYGEN_SHOULD_SKIP_THIS
template<class T_> class DiagPrecond : virtual public Precond<T_>
{

 public:

    DiagPrecond(const Matrix<T_> *a)
    {
      _size = a->getSize();
      _pivot.resize(_size);
      _a = a;
      int i;
      if ((i = inv_diag(_size, _a, _pivot)) != 0) {
        cerr << "\n*** Fatal Error in OFELI ***\nIn DiagPrecond::DiagPrecond(...) : ";
        cerr << "Zero detected in element " << i << endl;
        cerr << "Program stops\n";
        exit(1);
      }
#ifdef _OFELI_DEBUG
      std::clog << "An instance of class DiagPrecond is constructed.\n";
      std::clog << "File : " << __FILE__ << ", Line : " << __LINE__ << endl;
#endif
    }

    ~DiagPrecond()
    {
#ifdef _OFELI_DEBUG
       std::clog << "An instance of class DiagPrecond is constructed.\n";
       std::clog << "File : " << __FILE__ << ", Line : " << __LINE__ << endl;
#endif
    }

    void Solve(Vect<T_> &x) const
    {
       for (size_t i=0; i<_size; i++)
          x[i] *= _pivot[i];
    }

    void Solve(const Vect<T_> &b, Vect<T_> &x) const
    {
       for (size_t i=0; i<_size; i++)
          x[i] = b[i] * _pivot[i];
    }

    void TransSolve(Vect<T_> &x) const { Solve(x); }

    void TransSolve(const Vect<T_> &b, Vect<T_> &x) const { Solve(b,x); }

    const T_ & getDiag(size_t i) const { return _pivot[i-1]; }

    T_ & getPivot(size_t i) const { return _pivot[i-1]; }

   using Precond<T_>::_size;
   using Precond<T_>::_pivot;
   using Precond<T_>::_a;
};
#endif /* DOXYGEN_SHOULD_SKIP_THIS */


//==============================================================================
//                                CLASS ILUPrec
//==============================================================================

template<class T_, class M_> class ILUPrec
{

 public:

    ILUPrec() {
#ifdef _OFELI_DEBUG
       std::clog << "An instance of class ILUPrec is constructed.\n";
       std::clog << "File: " << __FILE__ << ", Line: " << __LINE__ << endl;
#endif
    }

    ILUPrec(const M_ &a) {
       setMatrix(a);
#ifdef _OFELI_DEBUG
       std::clog << "An instance of class ILUPrec is constructed.\n";
       std::clog << "File: " << __FILE__ << ", Line: " << __LINE__ << endl;
#endif
    }

    ~ILUPrec()
    {
#ifdef _OFELI_DEBUG
       std::clog << "An instance of class ILUPrec is destructed.\n";
       std::clog << "File : " << __FILE__ << ", Line : " << __LINE__ << endl;
#endif
    }

    void setMatrix(const M_ &a)
    {
       size_t i, j, k=1;
       _size = a.getSize();
       _pivot.resize(_size);
       _id.resize(_size);
       _a = &a;
       for (i=0; i<_size; i++) {
          for (j=_a->getRowPtr(i+1); j<_a->getRowPtr(i+2); j++,k++) {
             if (_a->getColInd(j)==i+1) {
                _id[i] = k;
                if ((*_a)(k)==0) {
                   cerr << "\n*** Fatal Error in OFELI ***\nIn ILUPrec::setMatrix(M_ &) : ";
                   cerr << "Zero detected in element " << i+1 << endl;
                   cerr << "Program stops\n";
                   exit(1);
                }
                _pivot[i] = (*_a)(k);
             }
          }
       }

       bool found;
       T_ c = 0;
       for (i=0; i<_size; i++) {
          _pivot[i] = 1./_pivot[i];
          for (j=_id[i]+1; j<=_a->getRowPtr(i+2)-1; j++) {
             found = false;
             size_t l = _a->getColInd(j);
             for (k=_a->getRowPtr(l); k<=_id[l-1]-1; k++) {
                if (_a->getColInd(k)==i+1) {
                   found = true;
                   c = (*_a)(k);
                }
             }
             if (found)
                _pivot[l-1] -= c * _pivot[i] * (*_a)(j);
          }
       }
    }

    void Solve(Vect<T_> &x) const
    {
       int i, j;
       Vect<T_> z(_size);
       for (i=0; i<_size; i++) {
          T_ s = 0;
          for (j=_a->getRowPtr(i+1); j<=_id[i]-1; j++)
             s += (*_a)[j-1] * z(_a->getColInd(j));
          z[i] = _pivot[i] * (x[i]-s);
       }
       for (i=_size-1; i>=0; i--) {
          T_ s = 0;
          for (j=_id[i]+1; j<=_a->getRowPtr(i+2)-1; j++) {
             s += (*_a)[j-1] * x(_a->getColInd(j));
             x[i] = z[i] - _pivot[i] * s;
          }
       }
    }

    void Solve(const Vect<T_> &b, Vect<T_> &x) const
    {
       int i, j;
       Vect<T_> z(_size);
       for (i=0; i<int(_size); i++) {
          T_ s = 0;
          for (j=int(_a->getRowPtr(i+1)); j<=int(_id[i])-1; j++)
             s += (*_a)[j-1] * z(_a->getColInd(j));
          z[i] = _pivot[i] * (b[i]-s);
       }
       for (i=int(_size-1); i>=0; i--) {
          T_ s = 0;
          for (j=int(_id[i]+1); j<=int(_a->getRowPtr(i+2))-1; j++)
             s += (*_a)[j-1] * x(_a->getColInd(j));
          x[i] = z[i] - _pivot[i] * s;
       }
    }

    void TransSolve(Vect<T_> &x) const
    {
       int i, j;
       Vect<T_> z(_size), temp(x);
       T_ tmp;

       for (i=0; i<_size; i++) {
          z[i] = temp[i];
          tmp = _pivot[i] * z[i];
          for (j=_id[i]+1; j<=_a->getRowPtr(i+2)-1; j++)
             temp(_a->getColInd(j)) -= tmp * (*_a)(j);
       }

       for (i=_size-1; i>=0; i--) {
          x[i] = _pivot[i] * z[i];
          for (j=_a->getRowPtr(i+1); j<=_id[i]-1; j++)
             z(_a->getColInd(j)) -= (*_a)(j) * x[i];
       }
    }

    void TransSolve(const Vect<T_> &b, Vect<T_> &x) const { x = b; TransSolve(x); }

    T_ & getPivot(size_t i) const { return _pivot[i-1]; }

 private:
  valarray<T_>     _pivot;
  valarray<size_t> _id;
  const M_         *_a;
  size_t           _size;
};

//==============================================================================
//                                CLASS ILUPrecond
//==============================================================================

#ifndef DOXYGEN_SHOULD_SKIP_THIS
template<class T_> class ILUPrecond : virtual public Precond<T_>
{

 public:

    ILUPrecond(const Matrix<T_> *A)
    {
       size_t i, j, k=1;
       _size = A->getSize();
       _pivot.resize(_size);
       _id.resize(_size);
       _a = A;

       for (i=0; i<_size; i++) {
          for (j=_a->getRowPtr(i+1); j<_a->getRowPtr(i+2); j++) {
             size_t jj = _a->getColInd(j);
             if (jj==i+1) {
                _id[i] = k;
                if ((*_a)(i+1,i+1)==0) {
                   cerr << "\n*** Fatal Error in OFELI ***\nIn ILUPrecond::ILUPrecond(...) : ";
                   cerr << "Zero detected in element " << i+1 << endl;
                   cerr << "Program stops\n";
                   exit(1);
                }
                _pivot[i] = (*_a)(i+1,i+1);
             }
             k++;
          }
       }
       T_ c = 0;
       for (i=0; i<_size; i++) {
          _pivot[i] = T_(1)/_pivot[i];
          for (j=_id[i]+1; j<=_a->getRowPtr(i+2)-1; j++) {
             bool found = false;
             size_t l = _a->getColInd(j);
             for (k=_a->getRowPtr(l); k<=_id[l-1]-1; k++) {
                if (_a->getColInd(k)==i+1) {
                   found = true;
                   c = (*_a)(i+1,i+1);
                }
             }
             if (found)
                _pivot[l-1] -= c * _pivot[i] * (*_a)(i+1,l);
          }
       }
#ifdef _OFELI_DEBUG
       std::clog << "An instance of class ILUPrecond is constructed.\n";
       std::clog << "File : " << __FILE__ << ", Line : " << __LINE__ << endl;
#endif
    }

    ~ILUPrecond()
    {
#ifdef _OFELI_DEBUG
       std::clog << "An instance of class ILUPrecond is destructed.\n";
       std::clog << "File : " << __FILE__ << ", Line : " << __LINE__ << endl;
#endif
    }

    void Solve(Vect<T_> &x) const
    {
       size_t i, j;
       Vect<T_> z(_size);
       for (i=0; i<_size; i++) {
          T_ s = 0;
          for (j=_a->getRowPtr(i+1); j<=_id[i]-1; j++) {
             size_t k = _a->getColInd(j);
             s += (*_a)(i+1,k) * z[k-1];
          }
          z[i] = _pivot[i] * (x[i]-s);
       }

       for (int ii=int(_size-1); ii>=0; ii--) {
          T_ s = 0;
          for (j=_id[ii]+1; j<=_a->getRowPtr(ii+2)-1; j++) {
             size_t k = _a->getColInd(j);
             s += (*_a)(ii+1,k) * x[k-1];
          }
          x[ii] = z[ii] - _pivot[ii] * s;
       }
    }

    void Solve(const Vect<T_> &b, Vect<T_> &x) const { x = b; Solve(x); }

    void TransSolve(Vect<T_> &x) const
    {
//     This is just to avoid warning on unused variable
       if (x.size()==0) { }
    }

    void TransSolve(const Vect<T_> &b, Vect<T_> &x) const { x = b; TransSolve(x); }

    T_ & getPivot(size_t i) const { return _pivot[i-1]; }


    using Precond<T_>::_size;
    using Precond<T_>::_pivot;
    using Precond<T_>::_id;
    using Precond<T_>::_a;
};
#endif /* DOXYGEN_SHOULD_SKIP_THIS */


//                         I M P L E M E N T A T I O N                       //

template<class T_,class M_>
IdentPrec<T_,M_>::IdentPrec(const M_ &A)
{
   A.getNbRows();
#ifdef _OFELI_DEBUG
   std::clog << "An instance of class IdentPrec is constructed.\n";
   std::clog << "File : " << __FILE__ << ", Line : " << __LINE__ << endl;
#endif
}


template<class T_,class M_>
IdentPrec<T_,M_>::~IdentPrec()
{
#ifdef _OFELI_DEBUG
   std::clog << "An instance of class IdentPrec is destructed.\n";
   std::clog << "File: " << __FILE__ << ", Line: " << __LINE__ << endl;
#endif
}


template<class T_,class M_>
void IdentPrec<T_,M_>::Solve(Vect<T_> &x) const
{
   x.size();
}


template<class T_,class M_>
void IdentPrec<T_,M_>::Solve(const Vect<T_> &b, Vect<T_> &x) const
{
   x = b;
}

template<class T_,class M_>
void IdentPrec<T_,M_>::TransSolve(const Vect<T_> &b, Vect<T_> &x) const
{
   x = b;
}


#ifndef DOXYGEN_SHOULD_SKIP_THIS
template<class T_>
IdentPrecond<T_>::IdentPrecond(const Matrix<T_> *A)
{
#ifdef _OFELI_DEBUG
   std::clog << "An instance of class IdentPrecond is constructed.\n";
   std::clog << "File: " << __FILE__ << ", Line: " << __LINE__ << endl;
#endif
}
#endif /* DOXYGEN_SHOULD_SKIP_THIS */


#ifndef DOXYGEN_SHOULD_SKIP_THIS
template<class T_>
IdentPrecond<T_>::~IdentPrecond()
{
#ifdef _OFELI_DEBUG
   std::clog << "An instance of class IdentPrecond is destructed.\n";
   std::clog << "File: " << __FILE__ << ", Line: " << __LINE__ << endl;
#endif
}
#endif /* DOXYGEN_SHOULD_SKIP_THIS */


#ifndef DOXYGEN_SHOULD_SKIP_THIS
template<class T_>
void IdentPrecond<T_>::Solve(const Vect<T_> &b, Vect<T_> &x) const
{
   x = b;
}
#endif /* DOXYGEN_SHOULD_SKIP_THIS */


#ifndef DOXYGEN_SHOULD_SKIP_THIS
template<class T_>
void IdentPrecond<T_>::TransSolve(const Vect<T_> &b, Vect<T_> &x) const
{
   x = b;
}
#endif /* DOXYGEN_SHOULD_SKIP_THIS */

} /* namespace OFELI */

#endif

Copyright © 1998-2007 Rachid Touzani