---

ILU0Preconditioner.h

/****************************************************************************
 *                                                                          *
 *  Program: FDDlib                                                         *
 *  Version: 1.1                                                            *
 *                                                                          *
 *  Copyright (C) 2002 - 2004 by Eric Miller and Dana Brooks                *
 *  All rights reserved.                                                    *
 *                                                                          *
 *  This software is Version 1.1 of the fddlib tomography toolbox.          *
 *  It is not to be redistributed or used for any commercial purpose        *
 *  without the prior written consent of the authors and Northeastern       *
 *  University.                                                             *
 *                                                                          *
 *  This software is provided as is, and any express or implied warranties, *
 *  including but not limited to the implied warranty of merchantability    *
 *  and the implied warranty of fitness for a particular purpose, are dis-  *
 *  claimed.  In no event shall the authors or Northeastern University be   *
 *  liable for any direct, indirect, incidental, special, exemplary, or     *
 *  consequential damages (including but not limited to procurement of      * 
 *  substitute goods or services; loss of use, data, or profits; or busi-   *
 *  ness interruption) however caused and on any theory of liability,       *
 *  whether in contract, strict liability, or tort (including negligence or *
 *  otherwise) arising in any way from the use of this software, even if    *
 *  advised of the possibility of such damage.                              *
 *                                                                          *
 *  Portions of this code benefit from ideas from Kyle Guilbert, Greg       *
 *  Boverman, Derek Uluski, David Kaeli, and Jennifer Black                 *
 *                                                                          *
 ****************************************************************************/

#ifndef _ILU0_PRECONDITIONER_H_
#define _ILU0_PRECONDITIONER_H_

#include "Preconditioner.h"
#include "SparseRowMatrix.h"
#include "RealVector.h"
#include <stdio.h>
#include <string>

namespace FDDlib {

class ILU0Preconditioner : public Preconditioner
{
private:
  SparseRowMatrix<double> SRMD_;
  RealVector<double>          invDiag_;

public:
  ILU0Preconditioner() {len_ = 0;}
  ~ILU0Preconditioner() {}

  void init(SparseRowMatrix<double>& srmd) throw(std::string)
  {
    int    i,j,k,l,m;
    double v;

    //
    // We store the inverse of the diagonal of the matrix
    // to avoid redoing costly divisions.
    //

    SRMD_ = srmd;

    len_ = SRMD_.getNumRows();
    invDiag_.init(len_);

    for (i = 1; i < len_; i++)
    {
      for (j = 0; j < SRMD_.getNumCols(i); j++)
      {
        k = SRMD_.getCol(i,j);
        if (k < i)
        {
          v = SRMD_.val(i,k)/SRMD_.val(k,k);
          SRMD_.set(i,k,v);
          for (l = 0; l < SRMD_.getNumCols(i); l++)
          {
            m = SRMD_.getCol(i,l);
            if (m > k)
            {
              v = SRMD_.val(i,m) - SRMD_.val(i,k)*SRMD_.val(k,m);
              SRMD_.set(i,m,v);
            }
          }
        }
      }
    }

    for (i = 0; i < len_; i++)
    {
      invDiag_(i) = 1/SRMD_.val(i,i);
    }

  }


  template<class Vector>
  Vector lowerSolve(const Vector& x) const
  {
    int            i, j;
    int            cols, col;
    Vector         r(len_);
    SparseRowEntry<double>* ep;
    double         v;

    for (i = 0; i < len_; i++)
    {
      if (i == 0)
      {
        r(i) = invDiag_(i)*x(i);
        r(i) = x(i);
      }
      else
      {
        cols = SRMD_.getNumCols(i);
        ep = SRMD_.getRowEntry(i);
        v = 0;
        for (j = 0; j < cols; j++)
        {
          col = ep->column;
          if (col < i)
          {
            v += (ep->val)*r(col);
          }
          ep++;
        }
        v = x(i) - v;
        r(i) = v;
      }
    }

    return r;
  }

  template<class Vector>
  Vector upperSolve(const Vector& x) const
  {
    int            i, j;
    int            cols, col;
    Vector         r(len_);
    SparseRowEntry<double>* ep;
    double         v;

    for (i = len_ - 1; i >= 0; i--)
    {
      if (i == len_ - 1)
      {
        r(i) = invDiag_(i)*x(i);
      }
      else
      {
        cols = SRMD_.getNumCols(i);
        v = 0;
        ep = SRMD_.getRowEntry(i);
        for (j = 0; j < cols; j++)
        {
          col = ep->column;;
          if (col > i)
          {
            v += (ep->val)*r(col);
          }
          ep++;
        }
        v = x(i) - v;
        v *= invDiag_(i);
        r(i) = v;
      }
    }

    return r;
  }


  RealVector<double> solve(const RealVector<double>& x) const
  {
    RealVector<double> u(len_);
    RealVector<double> l(len_);

    l = lowerSolve(x);
    u = upperSolve(l);


    return u;

  }


};
} // end namespace
#endif

---