eigen_sparse_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 #ifdef LIBMESH_HAVE_EIGEN
23 
24 
25 // C++ includes
26 
27 // Local Includes
30 #include "libmesh/string_to_enum.h"
31 
32 namespace libMesh
33 {
34 
35 /*----------------------- functions ----------------------------------*/
36 template <typename T>
38 {
39  if (this->initialized())
40  {
41  this->_is_initialized = false;
42 
43  this->_solver_type = GMRES;
44  this->_preconditioner_type = ILU_PRECOND;
45  }
46 }
47 
48 
49 
50 template <typename T>
52 {
53  // Initialize the data structures if not done so already.
54  if (!this->initialized())
55  {
56  this->_is_initialized = true;
57  }
58 }
59 
60 
61 
62 template <typename T>
63 std::pair<unsigned int, Real>
65  NumericVector<T> &solution_in,
66  NumericVector<T> &rhs_in,
67  const double tol,
68  const unsigned int m_its)
69 {
70  START_LOG("solve()", "EigenSparseLinearSolver");
71  this->init ();
72 
73  // Make sure the data passed in are really Eigen types
74  EigenSparseMatrix<T>& matrix = libmesh_cast_ref<EigenSparseMatrix<T>&>(matrix_in);
75  EigenSparseVector<T>& solution = libmesh_cast_ref<EigenSparseVector<T>&>(solution_in);
76  EigenSparseVector<T>& rhs = libmesh_cast_ref<EigenSparseVector<T>&>(rhs_in);
77 
78  // Close the matrix and vectors in case this wasn't already done.
79  matrix.close();
80  solution.close();
81  rhs.close();
82 
83  std::pair<unsigned int, Real> retval(0,0.);
84 
85  // Solve the linear system
86  switch (this->_solver_type)
87  {
88  // Conjugate-Gradient
89  case CG:
90  {
91  Eigen::ConjugateGradient<EigenSM> solver (matrix._mat);
92  solver.setMaxIterations(m_its);
93  solver.setTolerance(tol);
94  solution._vec = solver.solveWithGuess(rhs._vec,solution._vec);
95  libMesh::out << "#iterations: " << solver.iterations() << std::endl;
96  libMesh::out << "estimated error: " << solver.error() << std::endl;
97  retval = std::make_pair(solver.iterations(), solver.error());
98  break;
99  }
100 
101  // Bi-Conjugate Gradient Stabilized
102  case BICGSTAB:
103  {
104  Eigen::BiCGSTAB<EigenSM> solver (matrix._mat);
105  solver.setMaxIterations(m_its);
106  solver.setTolerance(tol);
107  solution._vec = solver.solveWithGuess(rhs._vec,solution._vec);
108  libMesh::out << "#iterations: " << solver.iterations() << std::endl;
109  libMesh::out << "estimated error: " << solver.error() << std::endl;
110  retval = std::make_pair(solver.iterations(), solver.error());
111  break;
112  }
113 
114  // // Generalized Minimum Residual
115  // case GMRES:
116  // {
117  // libmesh_not_implemented();
118  // break;
119  // }
120 
121  // Unknown solver, use BICGSTAB
122  default:
123  {
124  libMesh::err << "ERROR: Unsupported Eigen Solver: "
125  << Utility::enum_to_string(this->_solver_type) << std::endl
126  << "Continuing with BICGSTAB" << std::endl;
127 
128  this->_solver_type = BICGSTAB;
129 
130  STOP_LOG("solve()", "EigenSparseLinearSolver");
131 
132  return this->solve (matrix,
133  solution,
134  rhs,
135  tol,
136  m_its);
137  }
138  }
139 
140  STOP_LOG("solve()", "EigenSparseLinearSolver");
141  return retval;
142 }
143 
144 
145 
146 template <typename T>
147 std::pair<unsigned int, Real>
149  NumericVector<T> &solution_in,
150  NumericVector<T> &rhs_in,
151  const double tol,
152  const unsigned int m_its)
153 {
154 
155  START_LOG("adjoint_solve()", "EigenSparseLinearSolver");
156 
157  libmesh_experimental();
158  EigenSparseMatrix<T> mat_trans(this->comm());
159  matrix_in.get_transpose(mat_trans);
160 
161  std::pair<unsigned int, Real> retval = this->solve (mat_trans,
162  solution_in,
163  rhs_in,
164  tol,
165  m_its);
166 
167  STOP_LOG("adjoint_solve()", "EigenSparseLinearSolver");
168 
169  return retval;
170 }
171 
172 
173 
174 
175 template <typename T>
176 std::pair<unsigned int, Real>
178  NumericVector<T>& /*solution_in*/,
179  NumericVector<T>& /*rhs_in*/,
180  const double /*tol*/,
181  const unsigned int /*m_its*/)
182 {
183  libmesh_not_implemented();
184  return std::make_pair(0,0.0);
185 }
186 
187 
188 
189 template <typename T>
190 std::pair<unsigned int, Real>
192  const SparseMatrix<T>& /*precond_matrix*/,
193  NumericVector<T>& /*solution_in*/,
194  NumericVector<T>& /*rhs_in*/,
195  const double /*tol*/,
196  const unsigned int /*m_its*/)
197 {
198  libmesh_not_implemented();
199  return std::make_pair(0,0.0);
200 }
201 
202 
203 
204 template <typename T>
206 {
207  libmesh_not_implemented();
208 
209  // switch (this->_preconditioner_type)
210  // {
211  // case IDENTITY_PRECOND:
212  // _precond_type = NULL; return;
213 
214  // case ILU_PRECOND:
215  // _precond_type = ILUPrecond; return;
216 
217  // case JACOBI_PRECOND:
218  // _precond_type = JacobiPrecond; return;
219 
220  // case SSOR_PRECOND:
221  // _precond_type = SSORPrecond; return;
222 
223 
224  // default:
225  // libMesh::err << "ERROR: Unsupported LASPACK Preconditioner: "
226  // << this->_preconditioner_type << std::endl
227  // << "Continuing with ILU" << std::endl;
228  // this->_preconditioner_type = ILU_PRECOND;
229  // this->set_laspack_preconditioner_type();
230  // }
231 }
232 
233 
234 
235 template <typename T>
237 {
238  libMesh::out << "print_converged_reason() is currently only supported"
239  << "with Petsc 2.3.1 and later." << std::endl;
240 }
241 
242 
243 
244 //------------------------------------------------------------------
245 // Explicit instantiations
246 template class EigenSparseLinearSolver<Number>;
247 
248 } // namespace libMesh
249 
250 
251 #endif // #ifdef LIBMESH_HAVE_EIGEN

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

Hosted By:
SourceForge.net Logo