linear_solver.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2018 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 
20 // Local Includes
21 #include "libmesh/auto_ptr.h" // libmesh_make_unique
23 #include "libmesh/linear_solver.h"
28 #include "libmesh/preconditioner.h"
29 #include "libmesh/sparse_matrix.h"
30 #include "libmesh/string_to_enum.h"
35 
36 namespace libMesh
37 {
38 
39 //------------------------------------------------------------------
40 // LinearSolver members
41 template <typename T>
43  ParallelObject (comm_in),
44  _solver_type (GMRES),
45  _preconditioner_type (ILU_PRECOND),
46  _is_initialized (false),
47  _preconditioner (nullptr),
48  same_preconditioner (false),
49  _solver_configuration(nullptr)
50 {
51 }
52 
53 
54 
55 template <typename T>
56 std::unique_ptr<LinearSolver<T>>
58  const SolverPackage solver_package)
59 {
60  // Avoid unused parameter warnings when no solver packages are enabled.
61  libmesh_ignore(comm);
62 
63  // Build the appropriate solver
64  switch (solver_package)
65  {
66 #ifdef LIBMESH_HAVE_LASPACK
67  case LASPACK_SOLVERS:
68  return libmesh_make_unique<LaspackLinearSolver<T>>(comm);
69 #endif
70 
71 
72 #ifdef LIBMESH_HAVE_PETSC
73  case PETSC_SOLVERS:
74  return libmesh_make_unique<PetscLinearSolver<T>>(comm);
75 #endif
76 
77 
78 #ifdef LIBMESH_TRILINOS_HAVE_AZTECOO
79  case TRILINOS_SOLVERS:
80  return libmesh_make_unique<AztecLinearSolver<T>>(comm);
81 #endif
82 
83 
84 #ifdef LIBMESH_HAVE_EIGEN
85  case EIGEN_SOLVERS:
86  return libmesh_make_unique<EigenSparseLinearSolver<T>>(comm);
87 #endif
88 
89  default:
90  libmesh_error_msg("ERROR: Unrecognized solver package: " << solver_package);
91  }
92 
93  return std::unique_ptr<LinearSolver<T>>();
94 }
95 
96 template <typename T>
99 {
100  if (_preconditioner)
101  return _preconditioner->type();
102 
103  return _preconditioner_type;
104 }
105 
106 template <typename T>
107 void
109 {
110  if (_preconditioner)
111  _preconditioner->set_type(pct);
112  else
113  _preconditioner_type = pct;
114 }
115 
116 template <typename T>
117 void
119 {
120  if (this->_is_initialized)
121  libmesh_error_msg("Preconditioner must be attached before the solver is initialized!");
122 
123  _preconditioner_type = SHELL_PRECOND;
124  _preconditioner = preconditioner;
125 }
126 
127 template <typename T>
128 void
130 {
131  same_preconditioner = reuse_flag;
132 }
133 
134 template <typename T>
135 void
136 LinearSolver<T>::restrict_solve_to(const std::vector<unsigned int> * const dofs,
137  const SubsetSolveMode /*subset_solve_mode*/)
138 {
139  if (dofs != nullptr)
140  libmesh_not_implemented();
141 }
142 
143 
144 template <typename T>
145 std::pair<unsigned int, Real> LinearSolver<T>::adjoint_solve (SparseMatrix<T> & mat,
146  NumericVector<T> & sol,
147  NumericVector<T> & rhs,
148  const double tol,
149  const unsigned int n_iter)
150 {
151  // Log how long the linear solve takes.
152  LOG_SCOPE("adjoint_solve()", "LinearSolver");
153 
154  // Take the discrete adjoint
155  mat.close();
156  mat.get_transpose(mat);
157 
158  // Call the solve function for the relevant linear algebra library and
159  // solve the transpose matrix
160  const std::pair<unsigned int, Real> totalrval = this->solve (mat, sol, rhs, tol, n_iter);
161 
162  // Now transpose back and restore the original matrix
163  // by taking the discrete adjoint
164  mat.get_transpose(mat);
165 
166  return totalrval;
167 }
168 
169 template <typename T>
171 {
172  LinearConvergenceReason reason = this->get_converged_reason();
173  libMesh::out << "Linear solver convergence/divergence reason: " << Utility::enum_to_string(reason) << std::endl;
174 }
175 
176 template <typename T>
178 {
179  _solver_configuration = &solver_configuration;
180 }
181 
182 //------------------------------------------------------------------
183 // Explicit instantiations
184 template class LinearSolver<Number>;
185 
186 
187 
188 } // namespace libMesh
EIGEN_SOLVERS
Definition: libmesh.C:246
TRILINOS_SOLVERS
Definition: libmesh.C:244
Provides a uniform interface to vector storage schemes for different linear algebra libraries...
Definition: diff_context.h:40
void libmesh_ignore(const Args &...)
LinearSolver(const libMesh::Parallel::Communicator &comm_in)
Definition: linear_solver.C:42
LASPACK_SOLVERS
Definition: libmesh.C:248
An object whose state is distributed along a set of processors.
std::string enum_to_string(const T e)
virtual void close()=0
virtual void get_transpose(SparseMatrix< T > &dest) const =0
OStreamProxy out(std::cout)