libMesh
eigen_sparse_linear_solver.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 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 #include "libmesh/libmesh_common.h"
21 
22 #ifdef LIBMESH_HAVE_EIGEN
23 
24 
25 // C++ includes
26 
27 // Local Includes
28 #include "libmesh/eigen_sparse_linear_solver.h"
29 #include "libmesh/libmesh_logging.h"
30 #include "libmesh/string_to_enum.h"
31 #include "libmesh/solver_configuration.h"
32 
33 // GMRES is an "unsupported" iterative solver in Eigen.
34 #include <unsupported/Eigen/IterativeSolvers>
35 
36 namespace libMesh
37 {
38 
39 template <typename T>
42  LinearSolver<T>(comm_in),
43  _comp_info(Eigen::Success)
44 {
45  // The GMRES _solver_type can be used in EigenSparseLinearSolver,
46  // however, the GMRES iterative solver is currently in the Eigen
47  // "unsupported" directory, so we use BICGSTAB as our default.
48  this->_solver_type = BICGSTAB;
49 }
50 
51 
52 
53 template <typename T>
55 {
56  if (this->initialized())
57  {
58  this->_is_initialized = false;
59 
60  this->_solver_type = BICGSTAB;
62  }
63 }
64 
65 
66 
67 template <typename T>
68 void EigenSparseLinearSolver<T>::init (const char * /*name*/)
69 {
70  // Initialize the data structures if not done so already.
71  if (!this->initialized())
72  {
73  this->_is_initialized = true;
74  }
75 }
76 
77 
78 
79 template <typename T>
80 std::pair<unsigned int, Real>
82  NumericVector<T> & solution_in,
83  NumericVector<T> & rhs_in,
84  const double tol,
85  const unsigned int m_its)
86 {
87  LOG_SCOPE("solve()", "EigenSparseLinearSolver");
88  this->init ();
89 
90  // Make sure the data passed in are really Eigen types
91  EigenSparseMatrix<T> & matrix = cast_ref<EigenSparseMatrix<T> &>(matrix_in);
92  EigenSparseVector<T> & solution = cast_ref<EigenSparseVector<T> &>(solution_in);
93  EigenSparseVector<T> & rhs = cast_ref<EigenSparseVector<T> &>(rhs_in);
94 
95  // Close the matrix and vectors in case this wasn't already done.
96  matrix.close();
97  solution.close();
98  rhs.close();
99 
100  std::pair<unsigned int, Real> retval(0,0.);
101 
102  // Solve the linear system
103  switch (this->_solver_type)
104  {
105  // Conjugate-Gradient
106  case CG:
107  {
108  Eigen::ConjugateGradient<EigenSM> solver (matrix._mat);
109  solver.setMaxIterations(m_its);
110  solver.setTolerance(tol);
111  solution._vec = solver.solveWithGuess(rhs._vec,solution._vec);
112  libMesh::out << "#iterations: " << solver.iterations() << std::endl;
113  libMesh::out << "estimated error: " << solver.error() << std::endl;
114  retval = std::make_pair(solver.iterations(), solver.error());
115  _comp_info = solver.info();
116  break;
117  }
118 
119  // Bi-Conjugate Gradient Stabilized
120  case BICGSTAB:
121  {
122  Eigen::BiCGSTAB<EigenSM> solver (matrix._mat);
123  solver.setMaxIterations(m_its);
124  solver.setTolerance(tol);
125  solution._vec = solver.solveWithGuess(rhs._vec,solution._vec);
126  libMesh::out << "#iterations: " << solver.iterations() << std::endl;
127  libMesh::out << "estimated error: " << solver.error() << std::endl;
128  retval = std::make_pair(solver.iterations(), solver.error());
129  _comp_info = solver.info();
130  break;
131  }
132 
133  // Generalized Minimum Residual
134  case GMRES:
135  {
136  Eigen::GMRES<EigenSM> solver (matrix._mat);
137  solver.setMaxIterations(m_its);
138  solver.setTolerance(tol);
139 
140  // If there is an int parameter called "gmres_restart" in the
141  // SolverConfiguration object, pass it to the Eigen GMRES
142  // solver.
143  if (this->_solver_configuration)
144  {
145  std::map<std::string, int>::iterator it =
146  this->_solver_configuration->int_valued_data.find("gmres_restart");
147 
148  if (it != this->_solver_configuration->int_valued_data.end())
149  solver.set_restart(it->second);
150  }
151 
152  libMesh::out << "Eigen GMRES solver, restart = " << solver.get_restart() << std::endl;
153  solution._vec = solver.solveWithGuess(rhs._vec, solution._vec);
154  libMesh::out << "#iterations: " << solver.iterations() << std::endl;
155  libMesh::out << "estimated error: " << solver.error() << std::endl;
156  retval = std::make_pair(solver.iterations(), solver.error());
157  _comp_info = solver.info();
158  break;
159  }
160 
161  case SPARSELU:
162  {
163  // SparseLU solver code adapted from:
164  // http://eigen.tuxfamily.org/dox-devel/classEigen_1_1SparseLU.html
165  //
166  // From Eigen docs:
167  // The input matrix A should be in a compressed and
168  // column-major form. Otherwise an expensive copy will be
169  // made. You can call the inexpensive makeCompressed() to get
170  // a compressed matrix.
171  //
172  // Note: we don't have a column-major storage format here, so
173  // I think a copy must be made in order to use SparseLU. It
174  // appears that we also have to call makeCompressed(),
175  // otherwise you get a segfault.
176  matrix._mat.makeCompressed();
177 
178  // Build the SparseLU solver object. Note, there is one other
179  // sparse direct solver available in Eigen:
180  //
181  // Eigen::SparseQR<EigenSM, Eigen::AMDOrdering<int>> solver;
182  //
183  // I've tested it, and it works, but it is much slower than
184  // SparseLU. The main benefit of SparseQR is that it can
185  // handle non-square matrices, but we don't allow non-square
186  // sparse matrices to be built in libmesh...
187  Eigen::SparseLU<EigenSM> solver;
188 
189  // Compute the ordering permutation vector from the structural pattern of the matrix.
190  solver.analyzePattern(matrix._mat);
191 
192  // Compute the numerical factorization
193  solver.factorize(matrix._mat);
194 
195  // Use the factors to solve the linear system
196  solution._vec = solver.solve(rhs._vec);
197 
198  // Set up the return value. The SparseLU solver doesn't
199  // support asking for the number of iterations or the final
200  // error, so we'll just report back 1 and 0, respectively.
201  retval = std::make_pair(/*n. iterations=*/1, /*error=*/0);
202 
203  // Store the success/failure reason and break out.
204  _comp_info = solver.info();
205  break;
206  }
207 
208  // Unknown solver, use BICGSTAB
209  default:
210  {
211  libMesh::err << "ERROR: Unsupported Eigen Solver: "
212  << Utility::enum_to_string(this->_solver_type) << std::endl
213  << "Continuing with BICGSTAB" << std::endl;
214 
215  this->_solver_type = BICGSTAB;
216 
217  return this->solve (matrix,
218  solution,
219  rhs,
220  tol,
221  m_its);
222  }
223  }
224 
225  return retval;
226 }
227 
228 
229 
230 template <typename T>
231 std::pair<unsigned int, Real>
233  NumericVector<T> & solution_in,
234  NumericVector<T> & rhs_in,
235  const double tol,
236  const unsigned int m_its)
237 {
238  LOG_SCOPE("adjoint_solve()", "EigenSparseLinearSolver");
239 
240  libmesh_experimental();
241  EigenSparseMatrix<T> mat_trans(this->comm());
242  matrix_in.get_transpose(mat_trans);
243 
244  std::pair<unsigned int, Real> retval = this->solve (mat_trans,
245  solution_in,
246  rhs_in,
247  tol,
248  m_its);
249 
250  return retval;
251 }
252 
253 
254 
255 
256 template <typename T>
257 std::pair<unsigned int, Real>
259  NumericVector<T> & /*solution_in*/,
260  NumericVector<T> & /*rhs_in*/,
261  const double /*tol*/,
262  const unsigned int /*m_its*/)
263 {
264  libmesh_not_implemented();
265  return std::make_pair(0,0.0);
266 }
267 
268 
269 
270 template <typename T>
271 std::pair<unsigned int, Real>
273  const SparseMatrix<T> & /*precond_matrix*/,
274  NumericVector<T> & /*solution_in*/,
275  NumericVector<T> & /*rhs_in*/,
276  const double /*tol*/,
277  const unsigned int /*m_its*/)
278 {
279  libmesh_not_implemented();
280  return std::make_pair(0,0.0);
281 }
282 
283 
284 
285 template <typename T>
287 {
288  libmesh_not_implemented();
289 
290  // switch (this->_preconditioner_type)
291  // {
292  // case IDENTITY_PRECOND:
293  // _precond_type = libmesh_nullptr; return;
294 
295  // case ILU_PRECOND:
296  // _precond_type = ILUPrecond; return;
297 
298  // case JACOBI_PRECOND:
299  // _precond_type = JacobiPrecond; return;
300 
301  // case SSOR_PRECOND:
302  // _precond_type = SSORPrecond; return;
303 
304 
305  // default:
306  // libMesh::err << "ERROR: Unsupported LASPACK Preconditioner: "
307  // << this->_preconditioner_type << std::endl
308  // << "Continuing with ILU" << std::endl;
309  // this->_preconditioner_type = ILU_PRECOND;
310  // this->set_laspack_preconditioner_type();
311  // }
312 }
313 
314 
315 
316 template <typename T>
318 {
319  std::map<Eigen::ComputationInfo, LinearConvergenceReason>::iterator it =
321 
322  // If later versions of Eigen start returning new enumerations,
323  // we'll need to add them to the map...
324  if (it == _convergence_reasons.end())
325  {
326  libmesh_warning("Warning: unknown Eigen::ComputationInfo: " \
327  << _comp_info \
328  << " returning CONVERGED_ITS." \
329  << std::endl);
330  return CONVERGED_ITS;
331  }
332  else
333  return it->second;
334 }
335 
336 
337 
338 //------------------------------------------------------------------
339 // Explicit instantiations
340 template class EigenSparseLinearSolver<Number>;
341 
342 } // namespace libMesh
343 
344 
345 #endif // #ifdef LIBMESH_HAVE_EIGEN
OStreamProxy err
SolverType _solver_type
Enum stating which type of iterative solver to use.
Eigen::ComputationInfo _comp_info
Store the result of the last solve.
virtual void get_transpose(SparseMatrix< T > &dest) const =0
Copies the transpose of the matrix into dest, which may be *this.
Encapsulates the MPI_Comm object.
Definition: parallel.h:657
bool initialized() const
Definition: linear_solver.h:85
This class provides an interface to Eigen iterative solvers that is compatible with the libMesh Linea...
Numeric vector.
Definition: dof_map.h:66
The libMesh namespace provides an interface to certain functionality in the library.
EigenSparseLinearSolver(const libMesh::Parallel::Communicator &comm_in LIBMESH_CAN_DEFAULT_TO_COMMWORLD)
Constructor.
PetscDiffSolver & solver
Generic sparse matrix.
Definition: dof_map.h:65
virtual void clear() libmesh_override
Release all memory and clear data structures.
This base class can be inherited from to provide interfaces to linear solvers from different packages...
Definition: linear_solver.h:58
virtual std::pair< unsigned int, Real > adjoint_solve(SparseMatrix< T > &matrix, NumericVector< T > &solution, NumericVector< T > &rhs, const double tol, const unsigned int m_its) libmesh_override
Call the Eigen solver to solve A^T x = b.
std::map< std::string, int > int_valued_data
Store integer solver parameters in this map, e.g.
SolverConfiguration * _solver_configuration
Optionally store a SolverOptions object that can be used to set parameters like solver type...
static std::map< Eigen::ComputationInfo, LinearConvergenceReason > _convergence_reasons
Static map between Eigen ComputationInfo enumerations and libMesh LinearConvergenceReason enumeration...
std::string enum_to_string(const T e)
virtual LinearConvergenceReason get_converged_reason() const libmesh_override
virtual unsigned int solve() libmesh_override
This method performs a solve.
void set_eigen_preconditioner_type()
Tells Eigen to use the user-specified preconditioner stored in _preconditioner_type.
const Parallel::Communicator & comm() const
OStreamProxy out
PreconditionerType _preconditioner_type
Enum stating with type of preconditioner to use.
LinearConvergenceReason
Linear solver convergence flags (taken from the PETSc flags)
virtual std::pair< unsigned int, Real > solve(SparseMatrix< T > &matrix, NumericVector< T > &solution, NumericVector< T > &rhs, const double tol, const unsigned int m_its) libmesh_override
Call the Eigen solver.
Generic shell matrix, i.e.
virtual void init(const char *name=libmesh_nullptr) libmesh_override
Initialize data structures if not done so already.
bool _is_initialized
Flag indicating if the data structures have been initialized.