libMesh
trilinos_aztec_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 #ifdef LIBMESH_TRILINOS_HAVE_AZTECOO
23 
24 
25 // C++ includes
26 
27 // Local Includes
28 #include "libmesh/libmesh_logging.h"
29 #include "libmesh/string_to_enum.h"
30 #include "libmesh/trilinos_aztec_linear_solver.h"
31 #include "libmesh/trilinos_epetra_matrix.h"
32 #include "libmesh/trilinos_epetra_vector.h"
33 
34 namespace libMesh
35 {
36 
37 
38 /*----------------------- functions ----------------------------------*/
39 template <typename T>
41 {
42  if (this->initialized())
43  {
44  this->_is_initialized = false;
45 
46  // Mimic PETSc default solver and preconditioner
47  this->_solver_type = GMRES;
48 
49  if (this->n_processors() == 1)
50  this->_preconditioner_type = ILU_PRECOND;
51  else
52  this->_preconditioner_type = BLOCK_JACOBI_PRECOND;
53  }
54 }
55 
56 
57 
58 template <typename T>
59 void AztecLinearSolver<T>::init (const char * /*name*/)
60 {
61  // Initialize the data structures if not done so already.
62  if (!this->initialized())
63  {
64  this->_is_initialized = true;
65 
66  _linear_solver = new AztecOO();
67 
68  set_solver_type();
69 
70  switch(this->_preconditioner_type)
71  {
72  case ILU_PRECOND:
73  _linear_solver->SetAztecOption(AZ_precond,AZ_dom_decomp);
74  _linear_solver->SetAztecOption(AZ_subdomain_solve,AZ_ilu);
75  break;
76 
78  _linear_solver->SetAztecOption(AZ_precond,AZ_Jacobi);
79  break;
80 
81  case ICC_PRECOND:
82  _linear_solver->SetAztecOption(AZ_precond,AZ_dom_decomp);
83  _linear_solver->SetAztecOption(AZ_subdomain_solve,AZ_icc);
84  break;
85 
86  case LU_PRECOND:
87  _linear_solver->SetAztecOption(AZ_precond,AZ_dom_decomp);
88  _linear_solver->SetAztecOption(AZ_subdomain_solve,AZ_lu);
89  break;
90 
91  default:
92  _linear_solver->SetAztecOption(AZ_precond,AZ_dom_decomp);
93  _linear_solver->SetAztecOption(AZ_subdomain_solve,AZ_ilu);
94  }
95  }
96 }
97 
98 
99 
100 
101 template <typename T>
102 std::pair<unsigned int, Real>
104  SparseMatrix<T> & precond_in,
105  NumericVector<T> & solution_in,
106  NumericVector<T> & rhs_in,
107  const double tol,
108  const unsigned int m_its)
109 {
110  LOG_SCOPE("solve()", "AztecLinearSolver");
111 
112  // Make sure the data passed in are really of Epetra types
113  EpetraMatrix<T> * matrix = cast_ptr<EpetraMatrix<T> *>(&matrix_in);
114  EpetraMatrix<T> * precond = cast_ptr<EpetraMatrix<T> *>(&precond_in);
115  EpetraVector<T> * solution = cast_ptr<EpetraVector<T> *>(&solution_in);
116  EpetraVector<T> * rhs = cast_ptr<EpetraVector<T> *>(&rhs_in);
117 
118  this->init();
119 
120  // Close the matrices and vectors in case this wasn't already done.
121  matrix->close ();
122  precond->close ();
123  solution->close ();
124  rhs->close ();
125 
126  _linear_solver->SetAztecOption(AZ_max_iter,m_its);
127  _linear_solver->SetAztecParam(AZ_tol,tol);
128 
129  Epetra_FECrsMatrix * emat = matrix->mat();
130  Epetra_Vector * esol = solution->vec();
131  Epetra_Vector * erhs = rhs->vec();
132 
133  _linear_solver->Iterate(emat, esol, erhs, m_its, tol);
134 
135  // return the # of its. and the final residual norm.
136  return std::make_pair(_linear_solver->NumIters(), _linear_solver->TrueResidual());
137 }
138 
139 
140 
141 template <typename T>
142 std::pair<unsigned int, Real>
146  const double,
147  const unsigned int)
148 {
149  libmesh_not_implemented();
150 }
151 
152 
153 
154 template <typename T>
155 std::pair<unsigned int, Real>
157  const SparseMatrix<T> &,
160  const double,
161  const unsigned int)
162 {
163  libmesh_not_implemented();
164 }
165 
166 
167 
168 template <typename T>
169 void AztecLinearSolver<T>::get_residual_history(std::vector<double> & /* hist */)
170 {
171  libmesh_not_implemented();
172 }
173 
174 
175 
176 
177 template <typename T>
179 {
180  return _linear_solver->TrueResidual();
181 }
182 
183 
184 
185 template <typename T>
187 {
188  const double *status = _linear_solver->GetAztecStatus();
189 
190  switch (static_cast<int>(status[AZ_why]))
191  {
192  case AZ_normal :
193  libMesh::out << "AztecOO converged.\n";
194  break;
195  case AZ_maxits :
196  libMesh::out << "AztecOO failed to converge within maximum iterations.\n";
197  break;
198  case AZ_param :
199  libMesh::out << "AztecOO failed to support a user-requested parameter.\n";
200  break;
201  case AZ_breakdown :
202  libMesh::out << "AztecOO encountered numerical breakdown.\n";
203  break;
204  case AZ_loss :
205  libMesh::out << "AztecOO encountered numerical loss of precision.\n";
206  break;
207  case AZ_ill_cond :
208  libMesh::out << "AztecOO encountered an ill-conditioned GMRES Hessian.\n";
209  break;
210  default:
211  libMesh::out << "AztecOO reported an unrecognized condition.\n";
212  break;
213  }
214 }
215 
216 
217 
218 template <typename T>
220 {
221  const double *status = _linear_solver->GetAztecStatus();
222 
223  switch (static_cast<int>(status[AZ_why]))
224  {
225  case AZ_normal :
226  return CONVERGED_RTOL_NORMAL;
227  case AZ_maxits :
228  return DIVERGED_ITS;
229  }
230  return DIVERGED_NULL;
231 }
232 
233 
234 
235 template <typename T>
237 {
238  switch (this->_solver_type)
239  {
240  case CG:
241  _linear_solver->SetAztecOption(AZ_solver, AZ_cg); return;
242 
243  case CGS:
244  _linear_solver->SetAztecOption(AZ_solver, AZ_cgs); return;
245 
246  case TFQMR:
247  _linear_solver->SetAztecOption(AZ_solver, AZ_tfqmr); return;
248 
249  case BICGSTAB:
250  _linear_solver->SetAztecOption(AZ_solver, AZ_bicgstab); return;
251 
252  case GMRES:
253  _linear_solver->SetAztecOption(AZ_solver, AZ_gmres); return;
254 
255  default:
256  libMesh::err << "ERROR: Unsupported AztecOO Solver: "
257  << Utility::enum_to_string(this->_solver_type) << std::endl
258  << "Continuing with AztecOO defaults" << std::endl;
259  }
260 }
261 
262 //------------------------------------------------------------------
263 // Explicit instantiations
264 template class AztecLinearSolver<Number>;
265 
266 } // namespace libMesh
267 
268 
269 
270 #endif // #ifdef LIBMESH_TRILINOS_HAVE_AZTECOO
OStreamProxy err
virtual void init(const char *name=libmesh_nullptr) libmesh_override
Initialize data structures if not done so already.
virtual void print_converged_reason() const libmesh_override
Prints a useful message about why the latest linear solve con(di)verged.
virtual void clear() libmesh_override
Release all memory and clear data structures.
void get_residual_history(std::vector< double > &hist)
Fills the input vector with the sequence of residual norms from the latest iterative solve...
Numeric vector.
Definition: dof_map.h:66
The libMesh namespace provides an interface to certain functionality in the library.
Generic sparse matrix.
Definition: dof_map.h:65
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.
MPI_Status status
Status object for querying messages.
Definition: parallel.h:176
virtual LinearConvergenceReason get_converged_reason() const libmesh_override
std::string enum_to_string(const T e)
This class provides an interface to AztecOO iterative solvers that is compatible with the libMesh Lin...
void set_solver_type()
Tells AztecOO to use the user-specified solver stored in _solver_type.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual 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) libmesh_override
Call the Aztec solver.
OStreamProxy out
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.
processor_id_type n_processors()
Definition: libmesh_base.h:88