linear_solver.h

Go to the documentation of this file.

// The libMesh Finite Element Library.
// Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner

// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.

// This library 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
// Lesser General Public License for more details.

// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA



#ifndef LIBMESH_LINEAR_SOLVER_H
#define LIBMESH_LINEAR_SOLVER_H

// Local includes
#include "libmesh/libmesh_common.h"
#include "libmesh/enum_subset_solve_mode.h" // SUBSET_ZERO
#include "libmesh/reference_counted_object.h"
#include "libmesh/libmesh.h"
#include "libmesh/parallel_object.h"

// C++ includes
#include <cstddef>
#include <vector>
#include <memory>
#include <optional>

namespace libMesh
{

// forward declarations
template <typename T> class SparseMatrix;
template <typename T> class NumericVector;
template <typename T> class ShellMatrix;
template <typename T> class Preconditioner;
class System;
class SolverConfiguration;
enum SolverPackage : int;
enum PreconditionerType : int;
enum SolverType : int;
enum LinearConvergenceReason : int;

template <typename T>
class LinearSolver : public ReferenceCountedObject<LinearSolver<T>>,
                     public ParallelObject
{
public:

  LinearSolver (const libMesh::Parallel::Communicator & comm_in);

  virtual ~LinearSolver ();

  static std::unique_ptr<LinearSolver<T>> build(const libMesh::Parallel::Communicator & comm_in,
                                                const SolverPackage solver_package = libMesh::default_solver_package());

  bool initialized () const { return _is_initialized; }

  virtual void clear () {}

  virtual void init (const char * name = nullptr) = 0;

  virtual void init_names (const System &) {}

  SolverType solver_type () const { return _solver_type; }

  void set_solver_type (const SolverType st)
  { _solver_type = st; }

  PreconditionerType preconditioner_type () const;

  void set_preconditioner_type (const PreconditionerType pct);

  void attach_preconditioner(Preconditioner<T> * preconditioner);

  virtual void reuse_preconditioner(bool );

  bool get_same_preconditioner();

  virtual void restrict_solve_to (const std::vector<unsigned int> * const dofs,
                                  const SubsetSolveMode subset_solve_mode=SUBSET_ZERO);

  virtual std::pair<unsigned int, Real> solve (SparseMatrix<T> &,  // System Matrix
                                               NumericVector<T> &, // Solution vector
                                               NumericVector<T> &, // RHS vector
                                               const std::optional<double> tol = std::nullopt,
                                               const std::optional<unsigned int> m_its = std::nullopt) = 0; // N. Iterations

  virtual std::pair<unsigned int, Real> adjoint_solve (SparseMatrix<T> &,  // System Matrix
                                                       NumericVector<T> &, // Solution vector
                                                       NumericVector<T> &, // RHS vector
                                                       const std::optional<double> tol = std::nullopt,
                                                       const std::optional<unsigned int> m_its = std::nullopt); // N. Iterations

  virtual std::pair<unsigned int, Real> solve (SparseMatrix<T> &,  // System Matrix
                                               SparseMatrix<T> &,  // Preconditioning Matrix
                                               NumericVector<T> &, // Solution vector
                                               NumericVector<T> &, // RHS vector
                                               const std::optional<double> tol = std::nullopt,
                                               const std::optional<unsigned int> m_its = std::nullopt) = 0; // N. Iterations

  std::pair<unsigned int, Real> solve (SparseMatrix<T> & matrix,
                                       SparseMatrix<T> * precond_matrix,
                                       NumericVector<T> &, // Solution vector
                                       NumericVector<T> &, // RHS vector
                                       const std::optional<double> tol = std::nullopt,
                                       const std::optional<unsigned int> m_its = std::nullopt); // N. Iterations



  virtual std::pair<unsigned int, Real> solve (const ShellMatrix<T> & shell_matrix,
                                               NumericVector<T> &, // Solution vector
                                               NumericVector<T> &, // RHS vector
                                               const std::optional<double> tol = std::nullopt,
                                               const std::optional<unsigned int> m_its = std::nullopt) = 0; // N. Iterations


  virtual std::pair<unsigned int, Real> solve (const ShellMatrix<T> & shell_matrix,
                                               const SparseMatrix<T> & precond_matrix,
                                               NumericVector<T> &, // Solution vector
                                               NumericVector<T> &, // RHS vector
                                               const std::optional<double> tol = std::nullopt,
                                               const std::optional<unsigned int> m_its = std::nullopt) = 0; // N. Iterations


  std::pair<unsigned int, Real> solve (const ShellMatrix<T> & matrix,
                                       const SparseMatrix<T> * precond_matrix,
                                       NumericVector<T> &, // Solution vector
                                       NumericVector<T> &, // RHS vector
                                       const std::optional<double> tol = std::nullopt,
                                       const std::optional<unsigned int> m_its = std::nullopt); // N. Iterations


  virtual void print_converged_reason() const;

  virtual LinearConvergenceReason get_converged_reason() const = 0;

  void set_solver_configuration(SolverConfiguration & solver_configuration);

protected:
  SolverType _solver_type;

  PreconditionerType _preconditioner_type;

  bool _is_initialized;

  Preconditioner<T> * _preconditioner;

  bool same_preconditioner;

  SolverConfiguration * _solver_configuration;

  double get_real_solver_setting(const std::string & setting_name,
                                 const std::optional<double> & setting,
                                 const std::optional<double> default_value = std::nullopt);

  int get_int_solver_setting(const std::string & setting_name,
                             const std::optional<int> & setting,
                             const std::optional<int> default_value = std::nullopt);
};




/*----------------------- inline functions ----------------------------------*/
template <typename T>
inline
LinearSolver<T>::~LinearSolver ()
{
  this->LinearSolver::clear ();
}

template <typename T>
inline
bool LinearSolver<T>::get_same_preconditioner()
{
  return same_preconditioner;
}

template <typename T>
inline
std::pair<unsigned int, Real>
LinearSolver<T>::solve (SparseMatrix<T> & mat,
                        SparseMatrix<T> * pc_mat,
                        NumericVector<T> & sol,
                        NumericVector<T> & rhs,
                        const std::optional<double> tol,
                        const std::optional<unsigned int> n_iter)
{
  if (pc_mat)
    return this->solve(mat, *pc_mat, sol, rhs, tol, n_iter);
  else
    return this->solve(mat, sol, rhs, tol, n_iter);
}


template <typename T>
inline
std::pair<unsigned int, Real>
LinearSolver<T>::solve (const ShellMatrix<T> & mat,
                        const SparseMatrix<T> * pc_mat,
                        NumericVector<T> & sol,
                        NumericVector<T> & rhs,
                        const std::optional<double> tol,
                        const std::optional<unsigned int> n_iter)
{
  if (pc_mat)
    return this->solve(mat, *pc_mat, sol, rhs, tol, n_iter);
  else
    return this->solve(mat, sol, rhs, tol, n_iter);
}


} // namespace libMesh


#endif // LIBMESH_LINEAR_SOLVER_H