Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 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
26 : #include "libmesh/eigen_sparse_linear_solver.h"
27 : #include "libmesh/libmesh_logging.h"
28 : #include "libmesh/enum_to_string.h"
29 : #include "libmesh/solver_configuration.h"
30 : #include "libmesh/enum_preconditioner_type.h"
31 : #include "libmesh/enum_solver_type.h"
32 :
33 : // GMRES is an "unsupported" iterative solver in Eigen.
34 : #include "libmesh/ignore_warnings.h"
35 : #include <unsupported/Eigen/IterativeSolvers>
36 : #include "libmesh/restore_warnings.h"
37 :
38 : namespace libMesh
39 : {
40 :
41 : template <typename T>
42 276 : EigenSparseLinearSolver<T>::
43 : EigenSparseLinearSolver(const Parallel::Communicator & comm_in) :
44 : LinearSolver<T>(comm_in),
45 276 : _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 276 : this->_solver_type = BICGSTAB;
51 276 : }
52 :
53 :
54 :
55 : template <typename T>
56 523 : void EigenSparseLinearSolver<T>::clear ()
57 : {
58 799 : if (this->initialized())
59 : {
60 587 : this->_is_initialized = false;
61 :
62 587 : this->_solver_type = BICGSTAB;
63 587 : this->_preconditioner_type = ILU_PRECOND;
64 : }
65 523 : }
66 :
67 :
68 :
69 : template <typename T>
70 6900 : void EigenSparseLinearSolver<T>::init (const char * /*name*/)
71 : {
72 : // Initialize the data structures if not done so already.
73 6900 : if (!this->initialized())
74 : {
75 587 : this->_is_initialized = true;
76 : }
77 6900 : }
78 :
79 :
80 :
81 : template <typename T>
82 : std::pair<unsigned int, Real>
83 3609 : EigenSparseLinearSolver<T>::solve (SparseMatrix<T> & matrix_in,
84 : NumericVector<T> & solution_in,
85 : NumericVector<T> & rhs_in,
86 : const std::optional<double> tol,
87 : const std::optional<unsigned int> m_its)
88 : {
89 0 : LOG_SCOPE("solve()", "EigenSparseLinearSolver");
90 3609 : this->init ();
91 :
92 : // Make sure the data passed in are really Eigen types
93 0 : EigenSparseMatrix<T> & matrix = cast_ref<EigenSparseMatrix<T> &>(matrix_in);
94 0 : EigenSparseVector<T> & solution = cast_ref<EigenSparseVector<T> &>(solution_in);
95 0 : EigenSparseVector<T> & rhs = cast_ref<EigenSparseVector<T> &>(rhs_in);
96 :
97 : // Close the matrix and vectors in case this wasn't already done.
98 0 : matrix.close();
99 0 : solution.close();
100 0 : rhs.close();
101 :
102 0 : std::pair<unsigned int, Real> retval(0,0.);
103 :
104 : // Eigen doesn't give us a solver base class? We'll just use a
105 : // generic lambda, then.
106 7176 : auto do_solve = [this, &rhs, &solution, tol, m_its]
107 0 : (auto & e_solver, std::string_view msg) {
108 7134 : const int max_its = this->get_int_solver_setting("max_its", m_its);
109 3567 : const double abs_tol = this->get_real_solver_setting("abs_tol", tol);
110 :
111 3567 : e_solver.setMaxIterations(max_its);
112 0 : e_solver.setTolerance(abs_tol);
113 0 : libMesh::out << msg << std::endl;
114 :
115 3567 : solution._vec = e_solver.solveWithGuess(rhs._vec,solution._vec);
116 :
117 0 : libMesh::out << "#iterations: " << e_solver.iterations() << " / " << max_its << std::endl;
118 0 : libMesh::out << "estimated error: " << e_solver.error() << " / " << abs_tol << std::endl;
119 3567 : _comp_info = e_solver.info();
120 3567 : return std::make_pair(e_solver.iterations(), e_solver.error());
121 : };
122 :
123 : using Eigen::DiagonalPreconditioner;
124 : using Eigen::IdentityPreconditioner;
125 : using Eigen::IncompleteCholesky;
126 : using Eigen::IncompleteLUT;
127 :
128 : // Solve the linear system
129 3609 : switch (this->_solver_type)
130 : {
131 : // Conjugate-Gradient
132 30 : case CG:
133 : {
134 0 : const int UPLO = Eigen::Lower|Eigen::Upper;
135 :
136 30 : switch (this->_preconditioner_type)
137 : {
138 0 : case IDENTITY_PRECOND:
139 : {
140 0 : Eigen::ConjugateGradient<EigenSM,UPLO,IdentityPreconditioner> solver (matrix._mat);
141 0 : retval = do_solve(solver, "Eigen CG solver without preconditioning");
142 0 : break;
143 : }
144 30 : case ILU_PRECOND:
145 : {
146 : Eigen::ConjugateGradient<EigenSM,UPLO,IncompleteLUT<Number,eigen_idx_type>>
147 30 : solver (matrix._mat);
148 30 : retval = do_solve(solver, "Eigen CG solver with Incomplete Cholesky preconditioning");
149 0 : break;
150 : }
151 0 : case ICC_PRECOND:
152 : {
153 : Eigen::ConjugateGradient<EigenSM,UPLO,
154 : IncompleteCholesky<Number,Eigen::Lower,Eigen::AMDOrdering<eigen_idx_type>>>
155 0 : solver (matrix._mat);
156 0 : retval = do_solve(solver, "Eigen CG solver with Incomplete Cholesky preconditioning");
157 0 : break;
158 : }
159 0 : default: // For our default let's use their default
160 : libmesh_warning("No EigenSparseLinearSolver support for " <<
161 : Utility::enum_to_string<PreconditionerType>(this->_preconditioner_type)
162 : << " preconditioning.");
163 : libmesh_fallthrough();
164 : case JACOBI_PRECOND:
165 : {
166 0 : Eigen::ConjugateGradient<EigenSM,UPLO,DiagonalPreconditioner<Number>> solver (matrix._mat);
167 0 : retval = do_solve(solver, "Eigen CG solver with Jacobi preconditioning");
168 0 : break;
169 : }
170 : }
171 0 : break;
172 : }
173 :
174 : // Bi-Conjugate Gradient Stabilized
175 3530 : case BICGSTAB:
176 : {
177 3530 : switch (this->_preconditioner_type)
178 : {
179 129 : case IDENTITY_PRECOND:
180 : {
181 129 : Eigen::BiCGSTAB<EigenSM, IdentityPreconditioner> solver (matrix._mat);
182 129 : retval = do_solve(solver, "Eigen BiCGStab solver");
183 0 : break;
184 : }
185 3401 : case ICC_PRECOND:
186 : case ILU_PRECOND:
187 : {
188 : Eigen::BiCGSTAB<EigenSM,IncompleteLUT<Number, eigen_idx_type>>
189 3401 : solver (matrix._mat);
190 3401 : retval = do_solve(solver, "Eigen BiCGSTAB solver with ILU preconditioning");
191 0 : break;
192 : }
193 0 : default: // For our default let's use their default
194 : libmesh_warning("No EigenSparseLinearSolver support for " <<
195 : Utility::enum_to_string<PreconditionerType>(this->_preconditioner_type)
196 : << " preconditioning.");
197 : libmesh_fallthrough();
198 : case JACOBI_PRECOND:
199 : {
200 0 : Eigen::BiCGSTAB<EigenSM,DiagonalPreconditioner<Number>> solver (matrix._mat);
201 0 : retval = do_solve(solver, "Eigen BiCGSTAB solver with Jacobi preconditioning");
202 0 : break;
203 : }
204 : }
205 0 : break;
206 : }
207 :
208 : // Generalized Minimum Residual
209 7 : case GMRES:
210 : {
211 14 : auto set_restart_and_solve = [this, &do_solve]
212 0 : (auto & gm_solver, std::string_view msg)
213 : {
214 : // If there is an int parameter called "gmres_restart" in the
215 : // SolverConfiguration object, pass it to the Eigen GMRES
216 : // solver.
217 7 : if (this->_solver_configuration)
218 0 : if (const auto it = this->_solver_configuration->int_valued_data.find("gmres_restart");
219 0 : it != this->_solver_configuration->int_valued_data.end())
220 0 : gm_solver.set_restart(it->second);
221 :
222 7 : std::ostringstream full_msg;
223 7 : full_msg << msg << ", restart = " << gm_solver.get_restart();
224 14 : return do_solve(gm_solver, full_msg.str());
225 7 : };
226 :
227 7 : switch (this->_preconditioner_type)
228 : {
229 1 : case IDENTITY_PRECOND:
230 : {
231 1 : Eigen::GMRES<EigenSM,IdentityPreconditioner> solver (matrix._mat);
232 1 : retval = set_restart_and_solve(solver, "Eigen GMRES solver without preconditioning");
233 0 : break;
234 : }
235 6 : case ICC_PRECOND:
236 : case ILU_PRECOND:
237 : {
238 : Eigen::GMRES<EigenSM,IncompleteLUT<Number, eigen_idx_type>>
239 6 : solver (matrix._mat);
240 6 : retval = set_restart_and_solve(solver, "Eigen GMRES solver with ILU preconditioning");
241 0 : break;
242 : }
243 0 : default: // For our default let's use their default
244 : libmesh_warning("No EigenSparseLinearSolver support for " <<
245 : Utility::enum_to_string<PreconditionerType>(this->_preconditioner_type)
246 : << " preconditioning.");
247 : libmesh_fallthrough();
248 : case JACOBI_PRECOND:
249 : {
250 0 : Eigen::GMRES<EigenSM,DiagonalPreconditioner<Number>> solver (matrix._mat);
251 0 : retval = set_restart_and_solve(solver, "Eigen CG solver with Jacobi preconditioning");
252 0 : break;
253 : }
254 : }
255 0 : break;
256 : }
257 :
258 29 : case SPARSELU:
259 : {
260 : // SparseLU solver code adapted from:
261 : // http://eigen.tuxfamily.org/dox-devel/classEigen_1_1SparseLU.html
262 : //
263 : // From Eigen docs:
264 : // The input matrix A should be in a compressed and
265 : // column-major form. Otherwise an expensive copy will be
266 : // made. You can call the inexpensive makeCompressed() to get
267 : // a compressed matrix.
268 : //
269 : // Note: we don't have a column-major storage format here, so
270 : // I think a copy must be made in order to use SparseLU. It
271 : // appears that we also have to call makeCompressed(),
272 : // otherwise you get a segfault.
273 29 : matrix._mat.makeCompressed();
274 :
275 : // Build the SparseLU solver object. Note, there is one other
276 : // sparse direct solver available in Eigen:
277 : //
278 : // Eigen::SparseQR<EigenSM, Eigen::AMDOrdering<int>> solver;
279 : //
280 : // I've tested it, and it works, but it is much slower than
281 : // SparseLU. The main benefit of SparseQR is that it can
282 : // handle non-square matrices, but we don't allow non-square
283 : // sparse matrices to be built in libmesh...
284 29 : Eigen::SparseLU<EigenSM> solver;
285 :
286 0 : libMesh::out << "Eigen Sparse LU solver" << std::endl;
287 :
288 : // Compute the ordering permutation vector from the structural pattern of the matrix.
289 29 : solver.analyzePattern(matrix._mat);
290 :
291 : // Compute the numerical factorization
292 29 : solver.factorize(matrix._mat);
293 :
294 : // Use the factors to solve the linear system
295 29 : solution._vec = solver.solve(rhs._vec);
296 :
297 : // Set up the return value. The SparseLU solver doesn't
298 : // support asking for the number of iterations or the final
299 : // error, so we'll just report back 1 and 0, respectively.
300 0 : retval = std::make_pair(/*n. iterations=*/1, /*error=*/0);
301 :
302 : // Store the success/failure reason and break out.
303 29 : _comp_info = solver.info();
304 0 : break;
305 29 : }
306 :
307 : // Unknown solver, use BICGSTAB
308 0 : default:
309 : {
310 0 : libMesh::err << "ERROR: Unsupported Eigen Solver: "
311 13 : << Utility::enum_to_string(this->_solver_type) << std::endl
312 0 : << "Continuing with BICGSTAB" << std::endl;
313 :
314 13 : this->_solver_type = BICGSTAB;
315 :
316 13 : return this->solve (matrix,
317 : solution,
318 : rhs,
319 : tol,
320 0 : m_its);
321 : }
322 : }
323 :
324 3596 : return retval;
325 : }
326 :
327 :
328 :
329 : template <typename T>
330 : std::pair<unsigned int, Real>
331 76 : EigenSparseLinearSolver<T>::adjoint_solve (SparseMatrix<T> & matrix_in,
332 : NumericVector<T> & solution_in,
333 : NumericVector<T> & rhs_in,
334 : const std::optional<double> tol,
335 : const std::optional<unsigned int> m_its)
336 : {
337 0 : LOG_SCOPE("adjoint_solve()", "EigenSparseLinearSolver");
338 :
339 : libmesh_experimental();
340 76 : EigenSparseMatrix<T> mat_trans(this->comm());
341 76 : matrix_in.get_transpose(mat_trans);
342 :
343 76 : std::pair<unsigned int, Real> retval = this->solve (mat_trans,
344 : solution_in,
345 : rhs_in,
346 : tol,
347 : m_its);
348 :
349 76 : return retval;
350 76 : }
351 :
352 :
353 :
354 :
355 : template <typename T>
356 : std::pair<unsigned int, Real>
357 0 : EigenSparseLinearSolver<T>::solve (const ShellMatrix<T> & /*shell_matrix*/,
358 : NumericVector<T> & /*solution_in*/,
359 : NumericVector<T> & /*rhs_in*/,
360 : const std::optional<double> /*tol*/,
361 : const std::optional<unsigned int> /*m_its*/)
362 : {
363 0 : libmesh_not_implemented();
364 : return std::make_pair(0,0.0);
365 : }
366 :
367 :
368 :
369 : template <typename T>
370 : std::pair<unsigned int, Real>
371 0 : EigenSparseLinearSolver<T>::solve (const ShellMatrix<T> & /*shell_matrix*/,
372 : const SparseMatrix<T> & /*precond_matrix*/,
373 : NumericVector<T> & /*solution_in*/,
374 : NumericVector<T> & /*rhs_in*/,
375 : const std::optional<double> /*tol*/,
376 : const std::optional<unsigned int> /*m_its*/)
377 : {
378 0 : libmesh_not_implemented();
379 : return std::make_pair(0,0.0);
380 : }
381 :
382 :
383 :
384 : template <typename T>
385 0 : void EigenSparseLinearSolver<T>::set_eigen_preconditioner_type ()
386 : {
387 0 : libmesh_not_implemented();
388 :
389 : // switch (this->_preconditioner_type)
390 : // {
391 : // case IDENTITY_PRECOND:
392 : // _precond_type = nullptr; return;
393 :
394 : // case ILU_PRECOND:
395 : // _precond_type = ILUPrecond; return;
396 :
397 : // case JACOBI_PRECOND:
398 : // _precond_type = JacobiPrecond; return;
399 :
400 : // case SSOR_PRECOND:
401 : // _precond_type = SSORPrecond; return;
402 :
403 :
404 : // default:
405 : // libMesh::err << "ERROR: Unsupported LASPACK Preconditioner: "
406 : // << this->_preconditioner_type << std::endl
407 : // << "Continuing with ILU" << std::endl;
408 : // this->_preconditioner_type = ILU_PRECOND;
409 : // this->set_laspack_preconditioner_type();
410 : // }
411 : }
412 :
413 :
414 :
415 : template <typename T>
416 208 : LinearConvergenceReason EigenSparseLinearSolver<T>::get_converged_reason() const
417 : {
418 : // If later versions of Eigen start returning new enumerations,
419 : // we'll need to add them to the map...
420 208 : if (auto it = _convergence_reasons.find(_comp_info);
421 0 : it == _convergence_reasons.end())
422 : {
423 : libmesh_warning("Warning: unknown Eigen::ComputationInfo: " \
424 : << _comp_info \
425 : << " returning CONVERGED_ITS." \
426 : << std::endl);
427 0 : return CONVERGED_ITS;
428 : }
429 : else
430 208 : return it->second;
431 : }
432 :
433 :
434 :
435 : //------------------------------------------------------------------
436 : // Explicit instantiations
437 : template class LIBMESH_EXPORT EigenSparseLinearSolver<Number>;
438 :
439 : } // namespace libMesh
440 :
441 :
442 : #endif // #ifdef LIBMESH_HAVE_EIGEN
|
