laspack_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 #if defined(LIBMESH_HAVE_LASPACK)
23 
24 
25 // Local Includes
28 #include "libmesh/string_to_enum.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>
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>
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>
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>
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>
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>
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>
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>
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>
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
virtual void print_converged_reason() const override
virtual void init(const char *name=nullptr) override
Provides a uniform interface to vector storage schemes for different linear algebra libraries...
Definition: diff_context.h:40
virtual void close() override
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 clear() 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 zero() override
bool initialized()
Definition: libmesh.C:258
OStreamProxy out(std::cout)
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