KrylovSolver.cpp

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: KrylovSolver.cpp,v 1.13 2008/03/02 22:28:23 delpinux Exp $


#include <BaseVector.hpp>
#include <BaseMatrix.hpp>

#include <Vector.hpp>
#include <SparseMatrix.hpp>

#include <KrylovSolver.hpp>

#include <Preconditioner.hpp>

#include <DiagPrecond.hpp>
#include <IncompleteCholeskiFactorization.hpp>
#include <IdentityPrecond.hpp>

#include <SpectralFEMPreconditioner.hpp>

#include <PDESystem.hpp>

#include <ConjugateGradient.hpp>
#include <BiConjugateGradient.hpp>
#include <BiConjugateGradientStabilized.hpp>

#include <FGMRES.hpp>
#include <GMRES.hpp>

#include <PETScKrylovSolver.hpp>

#include <MultiGrid.hpp>

#include <Timer.hpp>

void KrylovSolverDim(ReferenceCounting<Vector<real_t> > u,
                     const BaseMatrix& A, const BaseVector& b,
                     const Problem& problem,
                     const KrylovSolverOptions::Type& type,
                     const KrylovSolverOptions::PreconditionerType& pType,
                     const DegreeOfFreedomSet& degreeOfFreedomSet)
{
#warning temporay implementation
#ifdef    HAVE_PETSC
  ASSERT(A.type() == BaseMatrix::petscMatrix);
  PETScKrylovSolver pks(static_cast<const Vector<real_t>&>(b),
                        static_cast<const PETScMatrix&>(A),
                        static_cast<Vector<real_t>&>(u));
  return;
#endif // HAVE_PETSC
  ReferenceCounting<Preconditioner> P = 0;

  switch (pType) {
  case KrylovSolverOptions::none: {
    P = new IdentityPrecond(problem);
    break;
  }
  case KrylovSolverOptions::diagonal: {
    P = new DiagPrecond(problem, A);
    break;
  }
  case KrylovSolverOptions::incompleteCholeski: {
    P = new IncompleteCholeskiFactorization(problem, A);
    break;
  }
  case KrylovSolverOptions::multiGrid: {
    P = new MultiGrid(problem,
                      static_cast<const SparseMatrix&>(A),
                      degreeOfFreedomSet);
    break;
  }
  case KrylovSolverOptions::spectralFEM: {
    P = new SpectralFEMPreconditioner(problem);
    break;
  }
  default: {
    throw ErrorHandler(__FILE__,__LINE__,
                       "unexpected preconditioner type",
                       ErrorHandler::unexpected);
  }
  }
  ffout(2) << "- preconditioner: " << (*P).name() << '\n';
  ffout(2) << "- preconditioner initialization\n";
  (*P).initializes();
  ffout(2) << "  preconditioner initialization: done\n";

  switch (type) {
  case (KrylovSolverOptions::conjugateGradient): {
    ConjugateGradient cg(static_cast<const Vector<real_t>&>(b),
                         A, *P, *u);
    break;
  }
  case (KrylovSolverOptions::biConjugateGradient): {
    BiConjugateGradient bicg(static_cast<const Vector<real_t>&>(b),
                             A, *P, *u);
    break;
  }
  case (KrylovSolverOptions::biConjugateGradientStabilized): {
    BiConjugateGradientStabilized bicgstab(static_cast<const Vector<real_t>&>(b),
                                           A, *P, *u);
    break;
  }
  case (KrylovSolverOptions::fgmres): {
    FGMRES fgmres(static_cast<const Vector<real_t>&>(b),
                  A, *P, *u);
    break;
  }
  case (KrylovSolverOptions::gmres): {
    GMRES gmres(static_cast<const Vector<real_t>&>(b),
                A, *P, *u);
    break;
  }
//   case (KrylovSolverOptions::multiGrid): {
//     MultiGrid<1, real, real_t> mg(static_cast<const Vector<real_t>&>(b),
//                                  static_cast<SparseMatrix<real_t>&>(A),
//                                  static_cast<Vector<real_t>&>(u),
//                                  os);
//     break;
//   }
  case (KrylovSolverOptions::iterativeLUFactorization): {
//     IterativeLUFactorization<1,real,real_t>(static_cast<const Vector<real_t>&>(b),
//                                            static_cast<SparseMatrix<real_t>&>(A),
//                                            static_cast<Vector<real_t>&>(u));
    throw ErrorHandler(__FILE__,__LINE__,
                       "not implemented",
                       ErrorHandler::unexpected);
    break;
  }
  default: {
    throw ErrorHandler(__FILE__,__LINE__,
                       "unexpected linear solver type",
                       ErrorHandler::unexpected);
  }
  }
}

void KrylovSolver::solve(const Problem& problem,
                         ReferenceCounting<Vector<real_t> > u)
{
  // initializing the timer.
  Timer t;
  t.start();
  
  ffout(2) << "Krylov solver\n";
  ffout(2) << "- number of unknowns: " << u->size() << '\n';
  
  KrylovSolverDim(u, __A, __b, problem, __type, __pType,
                  __degreeOfFreedomSet);
  // measure the time
  t.stop();
  ffout(2) << "Krylov solver: done";
  ffout(3) << " [cost: "  << t << ']';
  ffout(2) << '\n';
}

---