PETScKSP.cpp

Go to the documentation of this file.

/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   This file is part of CEPS.
 
   CEPS 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 3 of the License, or
   (at your option) any later version.
 
   CEPS 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 CEPS (see file LICENSE at root of project).
   If not, see <https://www.gnu.org/licenses/>.
 
 
   Copyright 2019-2024 Inria, Universite de Bordeaux
 
   Authors, in alphabetical order:
 
     Pierre-Elliott BECUE, Florian CARO, Yves COUDIERE(*), Andjela DAVIDOVIC, 
     Charlie DOUANLA-LONTSI, Marc FUENTES, Mehdi JUHOOR, Michael LEGUEBE(*), 
     Pauline MIGERDITICHAN, Valentin PANNETIER(*), Nejib ZEMZEMI.
     * : currently active authors
 
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
#include "linearAlgebra/LinearSystem.hpp"
 
LinearSystem::LinearSystem()
{
  this->initialize();
}
 
LinearSystem::LinearSystem(InputParameters *params)
{
  this->initialize();
 
  if (ceps::isValidPtr(params))
    setupWithParameters(params);
}
 
LinearSystem::~LinearSystem()
{
  if (ceps::isProfiling())
    ceps::profilingLog() << "#   Average number of CG iterations: " << m_avIterNb/(1.*m_nSolv) << std::endl;
 
  KSPDestroy(&m_solver);
}
 
void
LinearSystem::solve(DVecPtr x)
{
  if (!m_isSetUp)
    setUp();
 
  CEPS_ABORT_IF(not x, "Provided solution vector is nullptr");
  CEPS_ABORT_IF(not m_rhsVector, "right-hand side vector has not been set.");
 
  // Global size
  CepsIndex rhsSize, xSize;
  m_rhsVector->getSize(&rhsSize);
  x->getSize(&xSize);
  CEPS_ABORT_IF(rhsSize != xSize,
    "right hand side and solution vectors " << std::endl
    << "have incompatible sizes. RHS has global size " << rhsSize << "\n"
    << "while solution has size " << xSize
  );
  // Local range
  CepsIndex rhsLo, rhsHi, xLo, xHi;
  m_rhsVector->getLocalRange(&rhsLo, &rhsHi);
  x->getLocalRange(&xLo, &xHi);
  CEPS_ABORT_IF(rhsLo != xLo or rhsHi != xHi,
    "right hand side and solution vectors" << std::endl
    << "have incompatible rangees. RHS has local range (" << rhsLo << "," << rhsHi << ")\n"
    << "while solution has range (" << xLo << "," << xHi << ")"
    );
 
  if (ceps::isProfiling ())
    getProfiler()->start("LSsolve", "solving linear system(s)");
 
  // solve system
  KSPSolve (m_solver,m_rhsVector->getVector(),x->getVector());
 
  if (ceps::isProfiling())
    getProfiler()->stop("LSsolve");
 
  CepsInt nbIterations;
  KSPGetIterationNumber(m_solver, &nbIterations);
 
  m_avIterNb += nbIterations;
  m_nSolv++;
 
  // check if converged
  KSPConvergedReason reason;
  KSPGetConvergedReason (m_solver, &reason);
  // PetscViewer viewer;
  // PetscViewerCreate(PETSC_COMM_WORLD,&viewer);
  // KSPConvergedReasonView(m_solver,viewer);
 
  static constexpr const char* description[] = {
    "",
    "",
    "KSP_DIVERGED_NULL: breakdown when solving the Hessenberg system within GMRES", 
    "KSP_DIVERGED_ITS: requested number of iterations", 
    "KSP_DIVERGED_DTOL: large increase in the residual norm",
    "KSP_DIVERGED_BREAKDOWN: breakdown in the Krylov method",
    "KSP_DIVERGED_BREAKDOWN_BICG: breakdown in the KSPBGCS Krylov method",
    "KSP_DIVERGED_NONSYMMETRIC: the operator or preconditioner was not symmetric for a KSPType that requires symmetry",
    "KSP_DIVERGED_INDEFINITE_PC: the preconditioner was indefinite for a KSPType that requires it be definite",
    "KSP_DIVERGED_NANORINF: a not a number of infinity was detected in a vector during the computation",
    "KSP_DIVERGED_INDEFINITE_MAT: the operator was indefinite for a KSPType that requires it be definite",
    "KSP_DIVERGED_PC_FAILED: the action of the preconditioner failed for some reason"};
  CEPS_ABORT_IF(reason<=-1,
    "KSP diverged: " << description[-reason]
    // << "Please refer to the KSPConvergedReason page of PETSc\n"
    // << "documentation for more explanations"
  );
}
 
void LinearSystem::
setLhsMatrix (DMatPtr lhs)
{
  CEPS_ABORT_IF(not lhs,"nullptr passed as lhs matrix");
  m_lhsMatrix = lhs;
  // forces call of setOperators function
  m_isSetUp = false;
}
 
void
LinearSystem::setPreconditioningMatrix(DMatPtr precond)
{
  // nullptr matrix can be given: LHS matrix will be used when solving
  m_precondMatrix = precond;
  m_isSetUp       = false;
}
 
void
LinearSystem::setRhsVector(DVecPtr rhs)
{
  CEPS_ABORT_IF(not rhs,"nullptr passed as rhs vector");
  m_rhsVector = rhs;
  return;
}
 
void
LinearSystem::setType(const CepsString &stype)
{
  CepsString S = ceps::toUpper(stype);
  if (S == "CG" or S == "CONJUGATE GRADIENT")
  {
    PC pc;
    KSPGetPC (m_solver, &pc);
    PCSetType (pc, PCBJACOBI);
    KSPSetType (m_solver, KSPCG);
  }
  else if (S == "BICGSTAB")
  {
    PC pc;
    KSPGetPC (m_solver, &pc);
    PCSetType (pc, PCBJACOBI);
    KSPSetType (m_solver, KSPBCGS);
  }
  else if (S == "GMRES")
  {
    PC pc;
    KSPGetPC (m_solver, &pc);
    PCSetType (pc, PCBJACOBI);
    KSPSetType (m_solver, KSPGMRES);
  }
  else if (S == "LU")
  {
    CEPS_ABORT_IF(ceps::isParallel(),"Impossible to use a direct solver in parallel.");
    PC pc;
    KSPGetPC (m_solver, &pc);
    PCSetType (pc, PCLU);
  }
  else
  {
    // unsupported solverType
    CEPS_WARNS (
      "unsupported solver type. Will use default value GMRES.\n" <<
      "   Valid types are CG, BICGSTAB, GMRES, and LU (LU for sequential runs only)"
    );
  }
}
 
void
LinearSystem::setUp()
{
  // check correct lhs
  CEPS_ABORT_IF(not m_lhsMatrix,"no LHS matrix was provided to the solver");
 
  // if precond is nullptr, default is to give LHS matrix as preconditioning matrix
  if (not m_precondMatrix)
    m_precondMatrix = m_lhsMatrix;
 
  KSPSetOperators (m_solver, m_lhsMatrix->getMatrix (), m_precondMatrix->getMatrix ());
  KSPSetUp (m_solver);
  m_isSetUp = true;
}
 
void
LinearSystem::setTolerances(CepsReal rTol, CepsReal aTol, CepsReal dTol, CepsUInt iter)
{
  m_relativeTolerance   = rTol;
  m_absoluteTolerance   = aTol;
  m_divergenceTolerance = dTol;
  m_maxIterations       = iter;
  KSPSetTolerances (m_solver, rTol, aTol, dTol, iter);
}
 
void
LinearSystem::setMaxNbIterations(CepsUInt maxIterations)
{
  m_maxIterations = maxIterations;
  CepsReal rtol, abstol, dtol;
 
  // get the previous tolerances, since PETSc API does not
  // have a function that changes only the number of iterations
  KSPGetTolerances (m_solver, &rtol, &abstol, &dtol, nullptr);
  KSPSetTolerances (m_solver, rtol, abstol, dtol, maxIterations);
}
 
void
LinearSystem::setRelativeTolerance(CepsReal rTol)
{
  m_relativeTolerance = rTol;
  CepsReal absTol, dTol;
  CepsInt  iter;
  KSPGetTolerances (m_solver, nullptr, &absTol, &dTol, &iter);
  KSPSetTolerances (m_solver, rTol, absTol, dTol, iter);
}
 
void
LinearSystem::setAbsoluteTolerance(CepsReal aTol)
{
  m_absoluteTolerance = aTol;
  CepsReal relTol, dtol;
  CepsInt  iter;
  KSPGetTolerances(m_solver, &relTol, nullptr, &dtol, &iter);
  KSPSetTolerances(m_solver, relTol, aTol, dtol, iter);
}
 
void
LinearSystem::setDivergenceTolerance(CepsReal dTol)
{
  m_divergenceTolerance = dTol;
  CepsReal absTol, relTol;
  CepsInt  iter;
  KSPGetTolerances (m_solver, &relTol, &absTol, nullptr, &iter);
  KSPSetTolerances (m_solver, relTol, absTol, dTol, iter);
}
 
CepsUInt
LinearSystem::getMaxNbIterations() const
{
  return m_maxIterations;
}
 
void
LinearSystem::initialize()
{
  KSPCreate(PETSC_COMM_WORLD, &m_solver);
  m_isSetUp       = false;
  m_lhsMatrix     = nullptr;
  m_rhsVector     = nullptr;
  m_precondMatrix = nullptr;
 
  // default solver is GMRES
  setType("GMRES");
 
  // Get PetscDefaults for tolerances
  CepsReal rTol, aTol, dTol;
  CepsInt  iter;
  KSPGetTolerances(m_solver, &rTol, &aTol, &dTol, &iter);
  m_relativeTolerance   = rTol;
  m_absoluteTolerance   = aTol;
  m_divergenceTolerance = dTol;
  m_maxIterations       = iter;
 
  // Initialize iteration counters
  m_nSolv    = 0;
  m_avIterNb = 0;
}