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