linear_solver.h

Go to the documentation of this file.

// The libMesh Finite Element Library.
// Copyright (C) 2002-2026 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_systems (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