linear_solver.h

Go to the documentation of this file.

// The libMesh Finite Element Library.
// Copyright (C) 2002-2014 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_convergence_flags.h"
#include "libmesh/enum_solver_package.h"
#include "libmesh/enum_solver_type.h"
#include "libmesh/enum_preconditioner_type.h"
#include "libmesh/enum_subset_solve_mode.h"
#include "libmesh/reference_counted_object.h"
#include "libmesh/libmesh.h"
#include "libmesh/parallel_object.h"
#include "libmesh/auto_ptr.h"

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

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;

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

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

  virtual ~LinearSolver ();

  static UniquePtr<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 = NULL) = 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 double,      // Stopping tolerance
                                               const unsigned int) = 0; // N. Iterations

  virtual std::pair<unsigned int, Real> adjoint_solve (SparseMatrix<T>&,  // System Matrix
                                                       NumericVector<T>&, // Solution vector
                                                       NumericVector<T>&, // RHS vector
                                                       const double,      // Stopping tolerance
                                                       const unsigned int); // 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 double,      // Stopping tolerance
                                               const unsigned int) = 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 double,      // Stopping tolerance
                                       const unsigned int); // N. Iterations



  virtual std::pair<unsigned int, Real> solve (const ShellMatrix<T>& shell_matrix,
                                               NumericVector<T>&, // Solution vector
                                               NumericVector<T>&, // RHS vector
                                               const double,      // Stopping tolerance
                                               const unsigned int) = 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 double,      // Stopping tolerance
                                               const unsigned int) = 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 double,      // Stopping tolerance
                                       const unsigned int); // N. Iterations


  virtual void print_converged_reason() const;

  virtual LinearConvergenceReason get_converged_reason() const = 0;

protected:


  SolverType _solver_type;

  PreconditionerType _preconditioner_type;

  bool _is_initialized;

  Preconditioner<T> * _preconditioner;

  bool same_preconditioner;

};




/*----------------------- inline functions ----------------------------------*/
template <typename T>
inline
LinearSolver<T>::LinearSolver (const libMesh::Parallel::Communicator &comm_in) :
  ParallelObject       (comm_in),
  _solver_type         (GMRES),
  _preconditioner_type (ILU_PRECOND),
  _is_initialized      (false),
  _preconditioner      (NULL),
  same_preconditioner  (false)
{
}



template <typename T>
inline
LinearSolver<T>::~LinearSolver ()
{
  this->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 double       tol,
                        const 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 double           tol,
                        const 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