laspack_linear_solver.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 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
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>
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  START_LOG("solve()", "LaspackLinearSolver");
116  this->init ();
117 
118  // Make sure the data passed in are really in Laspack types
119  LaspackMatrix<T>* matrix = libmesh_cast_ptr<LaspackMatrix<T>*>(&matrix_in);
120  LaspackVector<T>* solution = libmesh_cast_ptr<LaspackVector<T>*>(&solution_in);
121  LaspackVector<T>* rhs = libmesh_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  libMesh::err << "ERROR: LASPACK Error: " << std::endl;
272  WriteLASErrDescr(stdout);
273  libmesh_error();
274  }
275 
276  STOP_LOG("solve()", "LaspackLinearSolver");
277  // Get the convergence step # and residual
278  return std::make_pair(GetLastNoIter(), GetLastAccuracy());
279 }
280 
281 
282 
283 template <typename T>
284 std::pair<unsigned int, Real>
286  NumericVector<T> &solution_in,
287  NumericVector<T> &rhs_in,
288  const double tol,
289  const unsigned int m_its)
290 {
291  START_LOG("adjoint_solve()", "LaspackLinearSolver");
292  this->init ();
293 
294  // Make sure the data passed in are really in Laspack types
295  LaspackMatrix<T>* matrix = libmesh_cast_ptr<LaspackMatrix<T>*>(&matrix_in);
296  LaspackVector<T>* solution = libmesh_cast_ptr<LaspackVector<T>*>(&solution_in);
297  LaspackVector<T>* rhs = libmesh_cast_ptr<LaspackVector<T>*>(&rhs_in);
298 
299  // Zero-out the solution to prevent the solver from exiting in 0
300  // iterations (?)
301  //TODO:[BSK] Why does Laspack do this? Comment out this and try ex13...
302  solution->zero();
303 
304  // Close the matrix and vectors in case this wasn't already done.
305  matrix->close ();
306  solution->close ();
307  rhs->close ();
308 
309  // Set the preconditioner type
310  this->set_laspack_preconditioner_type ();
311 
312  // Set the solver tolerance
313  SetRTCAccuracy (tol);
314 
315  // Solve the linear system
316  switch (this->_solver_type)
317  {
318  // Conjugate-Gradient
319  case CG:
320  {
321  CGIter (Transp_Q(&matrix->_QMat),
322  &solution->_vec,
323  &rhs->_vec,
324  m_its,
325  _precond_type,
326  1.);
327  break;
328  }
329 
330  // Conjugate-Gradient Normalized
331  case CGN:
332  {
333  CGNIter (Transp_Q(&matrix->_QMat),
334  &solution->_vec,
335  &rhs->_vec,
336  m_its,
337  _precond_type,
338  1.);
339  break;
340  }
341 
342  // Conjugate-Gradient Squared
343  case CGS:
344  {
345  CGSIter (Transp_Q(&matrix->_QMat),
346  &solution->_vec,
347  &rhs->_vec,
348  m_its,
349  _precond_type,
350  1.);
351  break;
352  }
353 
354  // Bi-Conjugate Gradient
355  case BICG:
356  {
357  BiCGIter (Transp_Q(&matrix->_QMat),
358  &solution->_vec,
359  &rhs->_vec,
360  m_its,
361  _precond_type,
362  1.);
363  break;
364  }
365 
366  // Bi-Conjugate Gradient Stabilized
367  case BICGSTAB:
368  {
369  BiCGSTABIter (Transp_Q(&matrix->_QMat),
370  &solution->_vec,
371  &rhs->_vec,
372  m_its,
373  _precond_type,
374  1.);
375  break;
376  }
377 
378  // Quasi-Minimum Residual
379  case QMR:
380  {
381  QMRIter (Transp_Q(&matrix->_QMat),
382  &solution->_vec,
383  &rhs->_vec,
384  m_its,
385  _precond_type,
386  1.);
387  break;
388  }
389 
390  // Symmetric over-relaxation
391  case SSOR:
392  {
393  SSORIter (Transp_Q(&matrix->_QMat),
394  &solution->_vec,
395  &rhs->_vec,
396  m_its,
397  _precond_type,
398  1.);
399  break;
400  }
401 
402  // Jacobi Relaxation
403  case JACOBI:
404  {
405  JacobiIter (Transp_Q(&matrix->_QMat),
406  &solution->_vec,
407  &rhs->_vec,
408  m_its,
409  _precond_type,
410  1.);
411  break;
412  }
413 
414  // Generalized Minimum Residual
415  case GMRES:
416  {
417  SetGMRESRestart (30);
418  GMRESIter (Transp_Q(&matrix->_QMat),
419  &solution->_vec,
420  &rhs->_vec,
421  m_its,
422  _precond_type,
423  1.);
424  break;
425  }
426 
427  // Unknown solver, use GMRES
428  default:
429  {
430  libMesh::err << "ERROR: Unsupported LASPACK Solver: "
431  << Utility::enum_to_string(this->_solver_type) << std::endl
432  << "Continuing with GMRES" << std::endl;
433 
434  this->_solver_type = GMRES;
435 
436  return this->solve (*matrix,
437  *solution,
438  *rhs,
439  tol,
440  m_its);
441  }
442  }
443 
444  // Check for an error
445  if (LASResult() != LASOK)
446  {
447  libMesh::err << "ERROR: LASPACK Error: " << std::endl;
448  WriteLASErrDescr(stdout);
449  libmesh_error();
450  }
451 
452  STOP_LOG("adjoint_solve()", "LaspackLinearSolver");
453  // Get the convergence step # and residual
454  return std::make_pair(GetLastNoIter(), GetLastAccuracy());
455 }
456 
457 
458 
459 
460 template <typename T>
461 std::pair<unsigned int, Real>
463  NumericVector<T>& /*solution_in*/,
464  NumericVector<T>& /*rhs_in*/,
465  const double /*tol*/,
466  const unsigned int /*m_its*/)
467 {
468  libmesh_not_implemented();
469  return std::make_pair(0,0.0);
470 }
471 
472 
473 
474 template <typename T>
475 std::pair<unsigned int, Real>
477  const SparseMatrix<T>& /*precond_matrix*/,
478  NumericVector<T>& /*solution_in*/,
479  NumericVector<T>& /*rhs_in*/,
480  const double /*tol*/,
481  const unsigned int /*m_its*/)
482 {
483  libmesh_not_implemented();
484  return std::make_pair(0,0.0);
485 }
486 
487 
488 
489 template <typename T>
491 {
492  switch (this->_preconditioner_type)
493  {
494  case IDENTITY_PRECOND:
495  _precond_type = NULL; return;
496 
497  case ILU_PRECOND:
498  _precond_type = ILUPrecond; return;
499 
500  case JACOBI_PRECOND:
501  _precond_type = JacobiPrecond; return;
502 
503  case SSOR_PRECOND:
504  _precond_type = SSORPrecond; return;
505 
506 
507  default:
508  libMesh::err << "ERROR: Unsupported LASPACK Preconditioner: "
509  << this->_preconditioner_type << std::endl
510  << "Continuing with ILU" << std::endl;
511  this->_preconditioner_type = ILU_PRECOND;
512  this->set_laspack_preconditioner_type();
513  }
514 }
515 
516 
517 
518 template <typename T>
520 {
521  libMesh::out << "print_converged_reason() is currently only supported"
522  << "with Petsc 2.3.1 and later." << std::endl;
523 }
524 
525 
526 
527 //------------------------------------------------------------------
528 // Explicit instantiations
529 template class LaspackLinearSolver<Number>;
530 
531 } // namespace libMesh
532 
533 
534 #endif // #ifdef LIBMESH_HAVE_LASPACK

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:05 UTC

Hosted By:
SourceForge.net Logo