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 : #if defined(LIBMESH_HAVE_LASPACK)
23 :
24 :
25 : // Local Includes
26 : #include "libmesh/laspack_linear_solver.h"
27 : #include "libmesh/libmesh_logging.h"
28 : #include "libmesh/enum_to_string.h"
29 : #include "libmesh/enum_solver_type.h"
30 : #include "libmesh/enum_preconditioner_type.h"
31 : #include "libmesh/enum_convergence_flags.h"
32 :
33 : namespace libMesh
34 : {
35 :
36 : // #ifndef LIBMESH_USE_COMPLEX_NUMBERS
37 : // extern "C"
38 : // {
39 : // #endif
40 :
41 : // void print_iter_accuracy(int Iter,
42 : // _LPReal rNorm,
43 : // _LPReal bNorm,
44 : // IterIdType IterId)
45 : // /* put out accuracy reached after each solver iteration */
46 : // {
47 :
48 : // //FILE * out = fopen("residual.hist", "a");
49 : // static int icall=0;
50 :
51 : // if (!icall)
52 : // {
53 : // printf("Iter ||r||/||f||\n");
54 : // printf("------------------\n");
55 : // icall=1;
56 : // }
57 :
58 : // if ( Iter%1==0 && (IterId == CGIterId ||
59 : // IterId == CGNIterId ||
60 : // IterId == GMRESIterId ||
61 : // IterId == BiCGIterId ||
62 : // IterId == QMRIterId ||
63 : // IterId == CGSIterId ||
64 : // IterId == BiCGSTABIterId) )
65 : // {
66 : // if (!_LPIsZeroReal(bNorm))
67 : // printf("%d \t %g\n", Iter, rNorm/bNorm);
68 : // else
69 : // printf("%d (fnorm == 0)\n", Iter);
70 : // }
71 :
72 : // //fclose(out);
73 : // }
74 :
75 : // #ifndef LIBMESH_USE_COMPLEX_NUMBERS
76 : // }
77 : // #endif
78 :
79 : /*----------------------- functions ----------------------------------*/
80 : template <typename T>
81 0 : void LaspackLinearSolver<T>::clear ()
82 : {
83 0 : if (this->initialized())
84 : {
85 0 : this->_is_initialized = false;
86 :
87 0 : this->_solver_type = GMRES;
88 0 : this->_preconditioner_type = ILU_PRECOND;
89 : }
90 0 : }
91 :
92 :
93 :
94 : template <typename T>
95 0 : void LaspackLinearSolver<T>::init (const char * /* name */)
96 : {
97 : // Initialize the data structures if not done so already.
98 0 : if (!this->initialized())
99 : {
100 0 : this->_is_initialized = true;
101 : }
102 :
103 : // SetRTCAuxProc (print_iter_accuracy);
104 0 : }
105 :
106 :
107 :
108 : template <typename T>
109 : std::pair<unsigned int, Real>
110 0 : LaspackLinearSolver<T>::solve (SparseMatrix<T> & matrix_in,
111 : NumericVector<T> & solution_in,
112 : NumericVector<T> & rhs_in,
113 : const std::optional<double> tol,
114 : const std::optional<unsigned int> m_its)
115 : {
116 0 : LOG_SCOPE("solve()", "LaspackLinearSolver");
117 0 : this->init ();
118 :
119 : // Make sure the data passed in are really in Laspack types
120 0 : LaspackMatrix<T> * matrix = cast_ptr<LaspackMatrix<T> *>(&matrix_in);
121 0 : LaspackVector<T> * solution = cast_ptr<LaspackVector<T> *>(&solution_in);
122 0 : LaspackVector<T> * rhs = cast_ptr<LaspackVector<T> *>(&rhs_in);
123 :
124 0 : const int max_its = this->get_int_solver_setting("max_its", m_its);
125 0 : const double abs_tol = this->get_real_solver_setting("abs_tol", tol);
126 :
127 : // Zero-out the solution to prevent the solver from exiting in 0
128 : // iterations (?)
129 : //TODO:[BSK] Why does Laspack do this? Comment out this and try ex13...
130 0 : solution->zero();
131 :
132 : // Close the matrix and vectors in case this wasn't already done.
133 0 : matrix->close ();
134 0 : solution->close ();
135 0 : rhs->close ();
136 :
137 : // Set the preconditioner type
138 0 : this->set_laspack_preconditioner_type ();
139 :
140 : // Set the solver tolerance
141 0 : SetRTCAccuracy (abs_tol);
142 :
143 : // Solve the linear system
144 0 : switch (this->_solver_type)
145 : {
146 : // Conjugate-Gradient
147 0 : case CG:
148 : {
149 0 : CGIter (&matrix->_QMat,
150 : &solution->_vec,
151 : &rhs->_vec,
152 : max_its,
153 : _precond_type,
154 : 1.);
155 0 : break;
156 : }
157 :
158 : // Conjugate-Gradient Normalized
159 0 : case CGN:
160 : {
161 0 : CGNIter (&matrix->_QMat,
162 : &solution->_vec,
163 : &rhs->_vec,
164 : max_its,
165 : _precond_type,
166 : 1.);
167 0 : break;
168 : }
169 :
170 : // Conjugate-Gradient Squared
171 0 : case CGS:
172 : {
173 0 : CGSIter (&matrix->_QMat,
174 : &solution->_vec,
175 : &rhs->_vec,
176 : max_its,
177 : _precond_type,
178 : 1.);
179 0 : break;
180 : }
181 :
182 : // Bi-Conjugate Gradient
183 0 : case BICG:
184 : {
185 0 : BiCGIter (&matrix->_QMat,
186 : &solution->_vec,
187 : &rhs->_vec,
188 : max_its,
189 : _precond_type,
190 : 1.);
191 0 : break;
192 : }
193 :
194 : // Bi-Conjugate Gradient Stabilized
195 0 : case BICGSTAB:
196 : {
197 0 : BiCGSTABIter (&matrix->_QMat,
198 : &solution->_vec,
199 : &rhs->_vec,
200 : max_its,
201 : _precond_type,
202 : 1.);
203 0 : break;
204 : }
205 :
206 : // Quasi-Minimum Residual
207 0 : case QMR:
208 : {
209 0 : QMRIter (&matrix->_QMat,
210 : &solution->_vec,
211 : &rhs->_vec,
212 : max_its,
213 : _precond_type,
214 : 1.);
215 0 : break;
216 : }
217 :
218 : // Symmetric over-relaxation
219 0 : case SSOR:
220 : {
221 0 : SSORIter (&matrix->_QMat,
222 : &solution->_vec,
223 : &rhs->_vec,
224 : max_its,
225 : _precond_type,
226 : 1.);
227 0 : break;
228 : }
229 :
230 : // Jacobi Relaxation
231 0 : case JACOBI:
232 : {
233 0 : JacobiIter (&matrix->_QMat,
234 : &solution->_vec,
235 : &rhs->_vec,
236 : max_its,
237 : _precond_type,
238 : 1.);
239 0 : break;
240 : }
241 :
242 : // Generalized Minimum Residual
243 0 : case GMRES:
244 : {
245 0 : SetGMRESRestart (30);
246 0 : GMRESIter (&matrix->_QMat,
247 : &solution->_vec,
248 : &rhs->_vec,
249 : max_its,
250 : _precond_type,
251 : 1.);
252 0 : break;
253 : }
254 :
255 : // Unknown solver, use GMRES
256 0 : default:
257 : {
258 0 : libMesh::err << "ERROR: Unsupported LASPACK Solver: "
259 0 : << Utility::enum_to_string(this->_solver_type) << std::endl
260 0 : << "Continuing with GMRES" << std::endl;
261 :
262 0 : this->_solver_type = GMRES;
263 :
264 : return this->solve (*matrix,
265 : *solution,
266 : *rhs,
267 : tol,
268 0 : m_its);
269 : }
270 : }
271 :
272 : // Check for an error
273 0 : if (LASResult() != LASOK)
274 : {
275 0 : WriteLASErrDescr(stdout);
276 0 : libmesh_error_msg("Exiting after LASPACK Error!");
277 : }
278 :
279 : // Get the convergence step # and residual
280 0 : return std::make_pair(GetLastNoIter(), GetLastAccuracy());
281 : }
282 :
283 :
284 :
285 : template <typename T>
286 : std::pair<unsigned int, Real>
287 0 : LaspackLinearSolver<T>::adjoint_solve (SparseMatrix<T> & matrix_in,
288 : NumericVector<T> & solution_in,
289 : NumericVector<T> & rhs_in,
290 : const std::optional<double> tol,
291 : const std::optional<unsigned int> m_its)
292 : {
293 0 : LOG_SCOPE("adjoint_solve()", "LaspackLinearSolver");
294 0 : this->init ();
295 :
296 : // Make sure the data passed in are really in Laspack types
297 0 : LaspackMatrix<T> * matrix = cast_ptr<LaspackMatrix<T> *>(&matrix_in);
298 0 : LaspackVector<T> * solution = cast_ptr<LaspackVector<T> *>(&solution_in);
299 0 : LaspackVector<T> * rhs = cast_ptr<LaspackVector<T> *>(&rhs_in);
300 :
301 0 : const int max_its = this->get_int_solver_setting("max_its", m_its);
302 0 : const double abs_tol = this->get_real_solver_setting("abs_tol", tol);
303 :
304 : // Zero-out the solution to prevent the solver from exiting in 0
305 : // iterations (?)
306 : //TODO:[BSK] Why does Laspack do this? Comment out this and try ex13...
307 0 : solution->zero();
308 :
309 : // Close the matrix and vectors in case this wasn't already done.
310 0 : matrix->close ();
311 0 : solution->close ();
312 0 : rhs->close ();
313 :
314 : // Set the preconditioner type
315 0 : this->set_laspack_preconditioner_type ();
316 :
317 : // Set the solver tolerance
318 0 : SetRTCAccuracy (abs_tol);
319 :
320 : // Solve the linear system
321 0 : switch (this->_solver_type)
322 : {
323 : // Conjugate-Gradient
324 0 : case CG:
325 : {
326 0 : CGIter (Transp_Q(&matrix->_QMat),
327 : &solution->_vec,
328 : &rhs->_vec,
329 : max_its,
330 : _precond_type,
331 : 1.);
332 0 : break;
333 : }
334 :
335 : // Conjugate-Gradient Normalized
336 0 : case CGN:
337 : {
338 0 : CGNIter (Transp_Q(&matrix->_QMat),
339 : &solution->_vec,
340 : &rhs->_vec,
341 : max_its,
342 : _precond_type,
343 : 1.);
344 0 : break;
345 : }
346 :
347 : // Conjugate-Gradient Squared
348 0 : case CGS:
349 : {
350 0 : CGSIter (Transp_Q(&matrix->_QMat),
351 : &solution->_vec,
352 : &rhs->_vec,
353 : max_its,
354 : _precond_type,
355 : 1.);
356 0 : break;
357 : }
358 :
359 : // Bi-Conjugate Gradient
360 0 : case BICG:
361 : {
362 0 : BiCGIter (Transp_Q(&matrix->_QMat),
363 : &solution->_vec,
364 : &rhs->_vec,
365 : max_its,
366 : _precond_type,
367 : 1.);
368 0 : break;
369 : }
370 :
371 : // Bi-Conjugate Gradient Stabilized
372 0 : case BICGSTAB:
373 : {
374 0 : BiCGSTABIter (Transp_Q(&matrix->_QMat),
375 : &solution->_vec,
376 : &rhs->_vec,
377 : max_its,
378 : _precond_type,
379 : 1.);
380 0 : break;
381 : }
382 :
383 : // Quasi-Minimum Residual
384 0 : case QMR:
385 : {
386 0 : QMRIter (Transp_Q(&matrix->_QMat),
387 : &solution->_vec,
388 : &rhs->_vec,
389 : max_its,
390 : _precond_type,
391 : 1.);
392 0 : break;
393 : }
394 :
395 : // Symmetric over-relaxation
396 0 : case SSOR:
397 : {
398 0 : SSORIter (Transp_Q(&matrix->_QMat),
399 : &solution->_vec,
400 : &rhs->_vec,
401 : max_its,
402 : _precond_type,
403 : 1.);
404 0 : break;
405 : }
406 :
407 : // Jacobi Relaxation
408 0 : case JACOBI:
409 : {
410 0 : JacobiIter (Transp_Q(&matrix->_QMat),
411 : &solution->_vec,
412 : &rhs->_vec,
413 : max_its,
414 : _precond_type,
415 : 1.);
416 0 : break;
417 : }
418 :
419 : // Generalized Minimum Residual
420 0 : case GMRES:
421 : {
422 0 : SetGMRESRestart (30);
423 0 : GMRESIter (Transp_Q(&matrix->_QMat),
424 : &solution->_vec,
425 : &rhs->_vec,
426 : max_its,
427 : _precond_type,
428 : 1.);
429 0 : break;
430 : }
431 :
432 : // Unknown solver, use GMRES
433 0 : default:
434 : {
435 0 : libMesh::err << "ERROR: Unsupported LASPACK Solver: "
436 0 : << Utility::enum_to_string(this->_solver_type) << std::endl
437 0 : << "Continuing with GMRES" << std::endl;
438 :
439 0 : this->_solver_type = GMRES;
440 :
441 : return this->solve (*matrix,
442 : *solution,
443 : *rhs,
444 : tol,
445 0 : m_its);
446 : }
447 : }
448 :
449 : // Check for an error
450 0 : if (LASResult() != LASOK)
451 : {
452 0 : WriteLASErrDescr(stdout);
453 0 : libmesh_error_msg("Exiting after LASPACK Error!");
454 : }
455 :
456 : // Get the convergence step # and residual
457 0 : return std::make_pair(GetLastNoIter(), GetLastAccuracy());
458 : }
459 :
460 :
461 :
462 :
463 : template <typename T>
464 : std::pair<unsigned int, Real>
465 0 : LaspackLinearSolver<T>::solve (const ShellMatrix<T> & /*shell_matrix*/,
466 : NumericVector<T> & /*solution_in*/,
467 : NumericVector<T> & /*rhs_in*/,
468 : const std::optional<double> /*tol*/,
469 : const std::optional<unsigned int> /*m_its*/)
470 : {
471 0 : libmesh_not_implemented();
472 0 : return std::make_pair(0,0.0);
473 : }
474 :
475 :
476 :
477 : template <typename T>
478 : std::pair<unsigned int, Real>
479 0 : LaspackLinearSolver<T>::solve (const ShellMatrix<T> & /*shell_matrix*/,
480 : const SparseMatrix<T> & /*precond_matrix*/,
481 : NumericVector<T> & /*solution_in*/,
482 : NumericVector<T> & /*rhs_in*/,
483 : const std::optional<double> /*tol*/,
484 : const std::optional<unsigned int> /*m_its*/)
485 : {
486 0 : libmesh_not_implemented();
487 : return std::make_pair(0,0.0);
488 : }
489 :
490 :
491 :
492 : template <typename T>
493 0 : void LaspackLinearSolver<T>::set_laspack_preconditioner_type ()
494 : {
495 0 : switch (this->_preconditioner_type)
496 : {
497 0 : case IDENTITY_PRECOND:
498 0 : _precond_type = nullptr; return;
499 :
500 0 : case ILU_PRECOND:
501 0 : _precond_type = ILUPrecond; return;
502 :
503 0 : case JACOBI_PRECOND:
504 0 : _precond_type = JacobiPrecond; return;
505 :
506 0 : case SSOR_PRECOND:
507 0 : _precond_type = SSORPrecond; return;
508 :
509 :
510 0 : default:
511 0 : libMesh::err << "ERROR: Unsupported LASPACK Preconditioner: "
512 0 : << this->_preconditioner_type << std::endl
513 0 : << "Continuing with ILU" << std::endl;
514 0 : this->_preconditioner_type = ILU_PRECOND;
515 0 : this->set_laspack_preconditioner_type();
516 : }
517 : }
518 :
519 :
520 :
521 : template <typename T>
522 0 : void LaspackLinearSolver<T>::print_converged_reason() const
523 : {
524 0 : switch (LASResult())
525 : {
526 0 : case LASOK :
527 0 : libMesh::out << "Laspack converged.\n";
528 0 : break;
529 0 : default :
530 0 : libMesh::out << "Laspack diverged.\n";
531 : }
532 0 : libMesh::out << "Detailed reporting is currently only supported"
533 0 : << "with Petsc 2.3.1 and later." << std::endl;
534 0 : }
535 :
536 :
537 :
538 : template <typename T>
539 0 : LinearConvergenceReason LaspackLinearSolver<T>::get_converged_reason() const
540 : {
541 0 : switch (LASResult())
542 : {
543 0 : case LASOK : return CONVERGED_RTOL_NORMAL;
544 0 : default : return DIVERGED_NULL;
545 : }
546 : }
547 :
548 :
549 :
550 : //------------------------------------------------------------------
551 : // Explicit instantiations
552 : template class LIBMESH_EXPORT LaspackLinearSolver<Number>;
553 :
554 : } // namespace libMesh
555 :
556 :
557 : #endif // #ifdef LIBMESH_HAVE_LASPACK
|
