eigen_sparse_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 #include "libmesh/libmesh_common.h"
21 
22 #ifdef LIBMESH_HAVE_EIGEN
23 
24 
25 // Local Includes
28 #include "libmesh/string_to_enum.h"
32 
33 // GMRES is an "unsupported" iterative solver in Eigen.
35 #include <unsupported/Eigen/IterativeSolvers>
37 
38 namespace libMesh
39 {
40 
41 template <typename T>
44  LinearSolver<T>(comm_in),
45  _comp_info(Eigen::Success)
46 {
47  // The GMRES _solver_type can be used in EigenSparseLinearSolver,
48  // however, the GMRES iterative solver is currently in the Eigen
49  // "unsupported" directory, so we use BICGSTAB as our default.
50  this->_solver_type = BICGSTAB;
51 }
52 
53 
54 
55 template <typename T>
57 {
58  if (this->initialized())
59  {
60  this->_is_initialized = false;
61 
62  this->_solver_type = BICGSTAB;
63  this->_preconditioner_type = ILU_PRECOND;
64  }
65 }
66 
67 
68 
69 template <typename T>
70 void EigenSparseLinearSolver<T>::init (const char * /*name*/)
71 {
72  // Initialize the data structures if not done so already.
73  if (!this->initialized())
74  {
75  this->_is_initialized = true;
76  }
77 }
78 
79 
80 
81 template <typename T>
82 std::pair<unsigned int, Real>
84  NumericVector<T> & solution_in,
85  NumericVector<T> & rhs_in,
86  const double tol,
87  const unsigned int m_its)
88 {
89  LOG_SCOPE("solve()", "EigenSparseLinearSolver");
90  this->init ();
91 
92  // Make sure the data passed in are really Eigen types
93  EigenSparseMatrix<T> & matrix = cast_ref<EigenSparseMatrix<T> &>(matrix_in);
94  EigenSparseVector<T> & solution = cast_ref<EigenSparseVector<T> &>(solution_in);
95  EigenSparseVector<T> & rhs = cast_ref<EigenSparseVector<T> &>(rhs_in);
96 
97  // Close the matrix and vectors in case this wasn't already done.
98  matrix.close();
99  solution.close();
100  rhs.close();
101 
102  std::pair<unsigned int, Real> retval(0,0.);
103 
104  // Solve the linear system
105  switch (this->_solver_type)
106  {
107  // Conjugate-Gradient
108  case CG:
109  {
110  Eigen::ConjugateGradient<EigenSM> solver (matrix._mat);
111  solver.setMaxIterations(m_its);
112  solver.setTolerance(tol);
113  solution._vec = solver.solveWithGuess(rhs._vec,solution._vec);
114  libMesh::out << "#iterations: " << solver.iterations() << std::endl;
115  libMesh::out << "estimated error: " << solver.error() << std::endl;
116  retval = std::make_pair(solver.iterations(), solver.error());
117  _comp_info = solver.info();
118  break;
119  }
120 
121  // Bi-Conjugate Gradient Stabilized
122  case BICGSTAB:
123  {
124  Eigen::BiCGSTAB<EigenSM> solver (matrix._mat);
125  solver.setMaxIterations(m_its);
126  solver.setTolerance(tol);
127  solution._vec = solver.solveWithGuess(rhs._vec,solution._vec);
128  libMesh::out << "#iterations: " << solver.iterations() << std::endl;
129  libMesh::out << "estimated error: " << solver.error() << std::endl;
130  retval = std::make_pair(solver.iterations(), solver.error());
131  _comp_info = solver.info();
132  break;
133  }
134 
135  // Generalized Minimum Residual
136  case GMRES:
137  {
138  Eigen::GMRES<EigenSM> solver (matrix._mat);
139  solver.setMaxIterations(m_its);
140  solver.setTolerance(tol);
141 
142  // If there is an int parameter called "gmres_restart" in the
143  // SolverConfiguration object, pass it to the Eigen GMRES
144  // solver.
145  if (this->_solver_configuration)
146  {
147  auto it = this->_solver_configuration->int_valued_data.find("gmres_restart");
148 
149  if (it != this->_solver_configuration->int_valued_data.end())
150  solver.set_restart(it->second);
151  }
152 
153  libMesh::out << "Eigen GMRES solver, restart = " << solver.get_restart() << std::endl;
154  solution._vec = solver.solveWithGuess(rhs._vec, solution._vec);
155  libMesh::out << "#iterations: " << solver.iterations() << std::endl;
156  libMesh::out << "estimated error: " << solver.error() << std::endl;
157  retval = std::make_pair(solver.iterations(), solver.error());
158  _comp_info = solver.info();
159  break;
160  }
161 
162  case SPARSELU:
163  {
164  // SparseLU solver code adapted from:
165  // http://eigen.tuxfamily.org/dox-devel/classEigen_1_1SparseLU.html
166  //
167  // From Eigen docs:
168  // The input matrix A should be in a compressed and
169  // column-major form. Otherwise an expensive copy will be
170  // made. You can call the inexpensive makeCompressed() to get
171  // a compressed matrix.
172  //
173  // Note: we don't have a column-major storage format here, so
174  // I think a copy must be made in order to use SparseLU. It
175  // appears that we also have to call makeCompressed(),
176  // otherwise you get a segfault.
177  matrix._mat.makeCompressed();
178 
179  // Build the SparseLU solver object. Note, there is one other
180  // sparse direct solver available in Eigen:
181  //
182  // Eigen::SparseQR<EigenSM, Eigen::AMDOrdering<int>> solver;
183  //
184  // I've tested it, and it works, but it is much slower than
185  // SparseLU. The main benefit of SparseQR is that it can
186  // handle non-square matrices, but we don't allow non-square
187  // sparse matrices to be built in libmesh...
188  Eigen::SparseLU<EigenSM> solver;
189 
190  // Compute the ordering permutation vector from the structural pattern of the matrix.
191  solver.analyzePattern(matrix._mat);
192 
193  // Compute the numerical factorization
194  solver.factorize(matrix._mat);
195 
196  // Use the factors to solve the linear system
197  solution._vec = solver.solve(rhs._vec);
198 
199  // Set up the return value. The SparseLU solver doesn't
200  // support asking for the number of iterations or the final
201  // error, so we'll just report back 1 and 0, respectively.
202  retval = std::make_pair(/*n. iterations=*/1, /*error=*/0);
203 
204  // Store the success/failure reason and break out.
205  _comp_info = solver.info();
206  break;
207  }
208 
209  // Unknown solver, use BICGSTAB
210  default:
211  {
212  libMesh::err << "ERROR: Unsupported Eigen Solver: "
213  << Utility::enum_to_string(this->_solver_type) << std::endl
214  << "Continuing with BICGSTAB" << std::endl;
215 
216  this->_solver_type = BICGSTAB;
217 
218  return this->solve (matrix,
219  solution,
220  rhs,
221  tol,
222  m_its);
223  }
224  }
225 
226  return retval;
227 }
228 
229 
230 
231 template <typename T>
232 std::pair<unsigned int, Real>
234  NumericVector<T> & solution_in,
235  NumericVector<T> & rhs_in,
236  const double tol,
237  const unsigned int m_its)
238 {
239  LOG_SCOPE("adjoint_solve()", "EigenSparseLinearSolver");
240 
241  libmesh_experimental();
242  EigenSparseMatrix<T> mat_trans(this->comm());
243  matrix_in.get_transpose(mat_trans);
244 
245  std::pair<unsigned int, Real> retval = this->solve (mat_trans,
246  solution_in,
247  rhs_in,
248  tol,
249  m_its);
250 
251  return retval;
252 }
253 
254 
255 
256 
257 template <typename T>
258 std::pair<unsigned int, Real>
260  NumericVector<T> & /*solution_in*/,
261  NumericVector<T> & /*rhs_in*/,
262  const double /*tol*/,
263  const unsigned int /*m_its*/)
264 {
265  libmesh_not_implemented();
266  return std::make_pair(0,0.0);
267 }
268 
269 
270 
271 template <typename T>
272 std::pair<unsigned int, Real>
274  const SparseMatrix<T> & /*precond_matrix*/,
275  NumericVector<T> & /*solution_in*/,
276  NumericVector<T> & /*rhs_in*/,
277  const double /*tol*/,
278  const unsigned int /*m_its*/)
279 {
280  libmesh_not_implemented();
281  return std::make_pair(0,0.0);
282 }
283 
284 
285 
286 template <typename T>
288 {
289  libmesh_not_implemented();
290 
291  // switch (this->_preconditioner_type)
292  // {
293  // case IDENTITY_PRECOND:
294  // _precond_type = nullptr; return;
295 
296  // case ILU_PRECOND:
297  // _precond_type = ILUPrecond; return;
298 
299  // case JACOBI_PRECOND:
300  // _precond_type = JacobiPrecond; return;
301 
302  // case SSOR_PRECOND:
303  // _precond_type = SSORPrecond; return;
304 
305 
306  // default:
307  // libMesh::err << "ERROR: Unsupported LASPACK Preconditioner: "
308  // << this->_preconditioner_type << std::endl
309  // << "Continuing with ILU" << std::endl;
310  // this->_preconditioner_type = ILU_PRECOND;
311  // this->set_laspack_preconditioner_type();
312  // }
313 }
314 
315 
316 
317 template <typename T>
319 {
320  auto it = _convergence_reasons.find(_comp_info);
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
virtual std::pair< unsigned int, Real > solve(SparseMatrix< T > &matrix, NumericVector< T > &solution, NumericVector< T > &rhs, const double tol, const unsigned int m_its) override
virtual void close() override
EigenSparseLinearSolver(const libMesh::Parallel::Communicator &comm_in)
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) override
Provides a uniform interface to vector storage schemes for different linear algebra libraries...
Definition: diff_context.h:40
virtual void init(const char *name=nullptr) override
void init(triangulateio &t)
virtual LinearConvergenceReason get_converged_reason() const override
OStreamProxy err(std::cerr)
std::string enum_to_string(const T e)
virtual void get_transpose(SparseMatrix< T > &dest) const =0
bool initialized()
Definition: libmesh.C:258
OStreamProxy out(std::cout)