IncompleteCholeskiFactorization.hpp

Go to the documentation of this file.

//  This file is part of ff3d - http://www.freefem.org/ff3d
//  Copyright (C) 2001, 2002, 2003 Stéphane Del Pino

//  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; either version 2, or (at your option)
//  any later version.

//  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.  

//  $Id: IncompleteCholeskiFactorization.hpp,v 1.5 2007/05/20 23:02:49 delpinux Exp $


#ifndef INCOMPLETE_CHOLESKI_FACTORIZATION_HPP
#define INCOMPLETE_CHOLESKI_FACTORIZATION_HPP

#include <TinyVector.hpp>
#include <TinyMatrix.hpp>
#include <Vector.hpp>
#include <Preconditioner.hpp>

#include <ErrorHandler.hpp>

class IncompleteCholeskiFactorization
  : public Preconditioner
{
private:
  ReferenceCounting<SparseMatrix> __Lstar;

  const SparseMatrix& A;

public:
  std::string name() const
  {
    return "incomplete Choleski factorization";
  }

  void initializes();

  void computes(const Vector<real_t>& r ,
                Vector<real_t>& z) const;

  IncompleteCholeskiFactorization(const Problem& problem,
                                  const BaseMatrix& AA)
    : Preconditioner (problem,
                      incompleteCholeskiFactorization),
      A(static_cast<const SparseMatrix&>(AA))
  {
    if(AA.type() != BaseMatrix::sparseMatrix) {
      throw ErrorHandler(__FILE__, __LINE__,
                         "cannot use incomplete Choleski Preconditioner:\n"
                         "\tINCOMPATIBLE matrix type!",
                         ErrorHandler::normal);
    }
    __Lstar = new SparseMatrix();
    (*__Lstar).copyProfile(static_cast<const SparseMatrix&>(AA));
  }


  ~IncompleteCholeskiFactorization()
  {
    ;
  }
};


void IncompleteCholeskiFactorization::
computes(const Vector<real_t>& r ,
           Vector<real_t>& z) const
{
  const SparseMatrix& Lstar = *__Lstar;
  ASSERT ((Lstar.size() == r.size())&&(Lstar.size() == z.size()));

  z = r;

  for(size_t i=0; i<z.size(); ++i) {
    real_t sum = z[i];

    const SparseLine& line = Lstar.line(i);

    size_t j = 0;
    while (line[j].column()<(int)i) {
      sum -= line[j].value() * z[line[j].column()];
      ++j;
    }
    z[i] = sum/line[j].value();
  }

  for (int i=z.size()-1; i>=0; --i) {
    const SparseLine& line = Lstar.line(i);

    size_t j = 0;
    while (line[j].column()<i) // Find the diagonal term ...
      ++j;

    int ii = j; // The diagonal index in the SparseLine.
    j++;

    while ((j<line.size())&&(line[j].column() != -1)) {
      z[i] -= line[j].value()*z[line[j].column()];
      ++j;
    }
    z[i] /= line[ii].value();
  }
}

void IncompleteCholeskiFactorization::initializes()
{
  SparseMatrix& Lstar = *__Lstar;  
  const SparseMatrix& LstarRO = Lstar;
  
  for (size_t i=0; i < A.size(); ++i) {
    const SparseLine& Lline = Lstar.line(i);
    const SparseLine& Aline = A.line(i);
    size_t j = 0;

    real_t sumLik = 0;

    while ((Lline[j].column()<(int)i)&&(Lline[j].column()!=-1)) {
      sumLik += Lline[j].value() * Lline[j].value();
      j++;
    }

    // j contains the diagonal position in the SparseLine.
    const real_t& Lii
      = Lstar(i,i)
      = std::sqrt(real_t(std::abs(A(i,i) - sumLik)));

    if (Lstar(i,i) == 0.)
      Lstar(i,i) = 1.;

    j++;

    while(j<Aline.size()) {
      const size_t& I = i;
      const int& J = Aline[j].column();
      if(J == -1)
        break;

      real_t sumLkmLim = 0;

      size_t k=0;
      while (Aline[k].column()<(int)i) {
        const size_t& K = Aline[k].column();
        sumLkmLim
          += LstarRO(J,K) * Lline[k].value();
        k++;
      }
      
      Lstar(I,J)
        = Lstar(J,I)
        = (A(I,J) - sumLkmLim) / Lii;

      ++j;
    }
  }
}

#endif // INCOMPLETE_CHOLESKI_FACTORIZATION_HPP

---