---

SSORPreconditioner.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 _SSOR_PRECONDITIONER_H_
#define _SSOR_PRECONDITIONER_H_

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

namespace FDDlib {

class SSORPreconditioner : public Preconditioner
{
private:

  SparseRowMatrix<double>* SRMDp;

  RealVector<double>           invDiag_;

  RealVector<double>           diagonal_;

  double                       omega_;


public:
  SSORPreconditioner(double omega)
  {
    len_ = 0;
    omega_ = omega;
    SRMDp = (SparseRowMatrix<double>*) NULL;
  }

  ~SSORPreconditioner()
  {
  }

  void init(SparseRowMatrix<double>& SRMD) throw(std::string)
  {
    int i;

    SRMDp = &SRMD;
    len_ = SRMD.getNumRows();

    invDiag_.init(len_);
    diagonal_.init(len_);

    for (i = 0; i < len_; i++)
    {
      if (SRMD.val(i,i) == 0) {
        throw std::string("SSORPreconditioner::init: zero on matrix diagonal.");
      }
      else
      {
        invDiag_(i) = 1 / SRMD.val(i,i);
        diagonal_(i) = 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);
      }
      else
      {
        cols = SRMDp->getNumCols(i);
        ep = SRMDp->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) - omega_*v;
        v *= invDiag_(i);
        r(i) = v;
      }
    }

    return r;
  }

  template<class Vector>
  Vector diagonalMultiply(const Vector& x) const
  {
    Vector r(len_);
    int    i;

    for (i = 0; i < len_; i++)
    {
      r(i) = x(i)*diagonal_(i);
    }

    return r;
  }

  RealVector<double> upperSolve(const RealVector<double>& x) const
  {
    int            i, j;
    int            cols, col;
    RealVector<double>         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 = SRMDp->getNumCols(i);
        v = 0;
        ep = SRMDp->getRowEntry(i);
        for (j = 0; j < cols; j++)
        {
          col = ep->column;;
          if (col > i)
          {
            v += (ep->val)*r(col);
          }
          ep++;
        }
        v = x(i) - omega_*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_);
    RealVector<double> d(len_);

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

    return u;

  }

};

} // end namespace
#endif

---