PETScLinearSolver.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright (c) OpenGeoSys Community (opengeosys.org)
// SPDX-License-Identifier: BSD-3-Clause
 
#pragma once
 
#include <petscksp.h>
 
#include <algorithm>
#include <array>
#include <string>
 
#include "MathLib/LinAlg/LinearSolverBehaviour.h"
#include "PETScMatrix.h"
#include "PETScVector.h"
 
namespace MathLib
{
class PETScLinearSolver
{
public:
    PETScLinearSolver(std::string const& prefix,
                      std::string const& petsc_options);
 
    ~PETScLinearSolver()
    {
        PetscCallAbort(PETSC_COMM_WORLD, KSPDestroy(&solver_));
    }
    // TODO check if some args in LinearSolver interface can be made const&.
    bool solve(PETScMatrix& A, PETScVector& b, PETScVector& x);

    PetscInt getNumberOfIterations() const
    {
        PetscInt its = 0;
        PetscCallAbort(PETSC_COMM_WORLD, KSPGetIterationNumber(solver_, &its));
        return its;
    }

    double getElapsedTime() const { return elapsed_ctime_; }

    bool canSolveRectangular() const { return can_solve_rectangular_; }

    bool willCompute(
        MathLib::LinearSolverBehaviour const /*linear_solver_behaviour*/) const
    {
        return true;
    }
 
private:
    bool canSolverHandleRectangular(std::string_view ksp_type)
    {
        // List of KSP types that can handle rectangular matrices
        constexpr auto rectangular_solvers = std::to_array<std::string_view>(
            {KSPGMRES, KSPFGMRES, KSPBCGS, KSPBCGSL, KSPTFQMR, KSPCGS});
 
        return std::ranges::any_of(rectangular_solvers,
                                   [&](const auto& solver)
                                   { return solver == ksp_type; });
    }
 
    KSP solver_;  
    PC pc_;       
 
    bool can_solve_rectangular_ = false;
 
    double elapsed_ctime_ = 0.0;  
};
 
}  // namespace MathLib