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