petsc_linear_solver.h
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 #ifndef LIBMESH_PETSC_LINEAR_SOLVER_H
21 #define LIBMESH_PETSC_LINEAR_SOLVER_H
22 
23 #include "libmesh/libmesh_config.h"
24 
25 #ifdef LIBMESH_HAVE_PETSC
26 
27 #include "libmesh/petsc_macro.h"
28 
33 EXTERN_C_FOR_PETSC_BEGIN
34 #if PETSC_VERSION_LESS_THAN(2,2,0)
35 # include <petscsles.h>
36 #else
37 # include <petscksp.h>
38 #endif
39 EXTERN_C_FOR_PETSC_END
40 
41 // Local includes
42 #include "libmesh/linear_solver.h"
43 
44 // C++ includes
45 #include <cstddef>
46 #include <vector>
47 
48 //--------------------------------------------------------------------
49 // Functions with C linkage to pass to PETSc. PETSc will call these
50 // methods as needed for preconditioning
51 //
52 // Since they must have C linkage they have no knowledge of a namespace.
53 // Give them an obscure name to avoid namespace pollution.
54 extern "C"
55 {
56  // Older versions of PETSc do not have the different int typedefs.
57  // On 64-bit machines, PetscInt may actually be a long long int.
58  // This change occurred in Petsc-2.2.1.
59 #if PETSC_VERSION_LESS_THAN(2,2,1)
60  typedef int PetscErrorCode;
61  typedef int PetscInt;
62 #endif
63 
64 #if PETSC_RELEASE_LESS_THAN(3,0,1)
65 
69  PetscErrorCode __libmesh_petsc_preconditioner_setup (void * ctx);
70 
75  PetscErrorCode __libmesh_petsc_preconditioner_apply(void *ctx, Vec x, Vec y);
76 #else
77  PetscErrorCode __libmesh_petsc_preconditioner_setup (PC);
78  PetscErrorCode __libmesh_petsc_preconditioner_apply(PC, Vec x, Vec y);
79 #endif
80 } // end extern "C"
81 
82 
83 namespace libMesh
84 {
85 
86 // forward declarations
87 template <typename T> class PetscMatrix;
88 
97 template <typename T>
99 {
100 public:
105  LIBMESH_CAN_DEFAULT_TO_COMMWORLD);
106 
111 
115  void clear ();
116 
120  void init ();
121 
125  void init (PetscMatrix<T>* matrix);
126 
134  virtual void restrict_solve_to (const std::vector<unsigned int>* const dofs,
135  const SubsetSolveMode subset_solve_mode=SUBSET_ZERO);
136 
141  std::pair<unsigned int, Real>
142  solve (SparseMatrix<T> &matrix_in,
143  NumericVector<T> &solution_in,
144  NumericVector<T> &rhs_in,
145  const double tol,
146  const unsigned int m_its)
147  {
148  return this->solve(matrix_in, matrix_in, solution_in, rhs_in, tol, m_its);
149  }
150 
151 
156  std::pair<unsigned int, Real>
157  adjoint_solve (SparseMatrix<T> &matrix_in,
158  NumericVector<T> &solution_in,
159  NumericVector<T> &rhs_in,
160  const double tol,
161  const unsigned int m_its);
162 
163 
179  std::pair<unsigned int, Real>
180  solve (SparseMatrix<T> &matrix,
181  SparseMatrix<T> &preconditioner,
182  NumericVector<T> &solution,
183  NumericVector<T> &rhs,
184  const double tol,
185  const unsigned int m_its);
186 
190  std::pair<unsigned int, Real>
191  solve (const ShellMatrix<T>& shell_matrix,
192  NumericVector<T>& solution_in,
193  NumericVector<T>& rhs_in,
194  const double tol,
195  const unsigned int m_its);
196 
202  virtual std::pair<unsigned int, Real>
203  solve (const ShellMatrix<T>& shell_matrix,
204  const SparseMatrix<T>& precond_matrix,
205  NumericVector<T>& solution_in,
206  NumericVector<T>& rhs_in,
207  const double tol,
208  const unsigned int m_its);
209 
215  PC pc() { this->init(); return _pc; }
216 
222  KSP ksp() { this->init(); return _ksp; }
223 
228  void get_residual_history(std::vector<double>& hist);
229 
237 
242  virtual void print_converged_reason();
243 
244 private:
245 
250  void set_petsc_solver_type ();
251 
255  static PetscErrorCode _petsc_shell_matrix_mult(Mat mat, Vec arg, Vec dest);
256 
260  static PetscErrorCode _petsc_shell_matrix_mult_add(Mat mat, Vec arg, Vec add, Vec dest);
261 
265  static PetscErrorCode _petsc_shell_matrix_get_diagonal(Mat mat, Vec dest);
266 
267  // SLES removed from >= PETSc 2.2.0
268 #if PETSC_VERSION_LESS_THAN(2,2,0)
269 
273  SLES _sles;
274 
275 #endif
276 
280  PC _pc;
281 
285  KSP _ksp;
286 
292 
299 
304  size_t _restrict_solve_to_is_local_size(void)const;
305 
311  void _create_complement_is (const NumericVector<T> &vec_in);
312 
318 
319 };
320 
321 
322 /*----------------------- functions ----------------------------------*/
323 template <typename T>
324 inline
326  LinearSolver<T>(comm),
327  _restrict_solve_to_is(NULL),
328  _restrict_solve_to_is_complement(NULL),
329  _subset_solve_mode(SUBSET_ZERO)
330 {
331  if (this->n_processors() == 1)
333  else
335 }
336 
337 
338 
339 template <typename T>
340 inline
342 {
343  this->clear ();
344 }
345 
346 
347 
348 template <typename T>
349 inline size_t
352 {
353  libmesh_assert(_restrict_solve_to_is);
354 
355  PetscInt s;
356  int ierr = ISGetLocalSize(_restrict_solve_to_is,&s);
357  LIBMESH_CHKERRABORT(ierr);
358 
359  return static_cast<size_t>(s);
360 }
361 
362 
363 
364 template <typename T>
365 void
367 #if PETSC_VERSION_LESS_THAN(3,0,0)
368  // unnamed to avoid compiler "unused parameter" warning
369 #else
370  vec_in
371 #endif
372  )
373 {
374  libmesh_assert(_restrict_solve_to_is);
375 #if PETSC_VERSION_LESS_THAN(3,0,0)
376  // No ISComplement in PETSc 2.3.3
377  libmesh_not_implemented();
378 #else
379  if(_restrict_solve_to_is_complement==NULL)
380  {
381  int ierr = ISComplement(_restrict_solve_to_is,
382  vec_in.first_local_index(),
383  vec_in.last_local_index(),
384  &_restrict_solve_to_is_complement);
385  LIBMESH_CHKERRABORT(ierr);
386  }
387 #endif
388 }
389 
390 
391 
392 } // namespace libMesh
393 
394 
395 #endif // #ifdef LIBMESH_HAVE_PETSC
396 #endif // LIBMESH_PETSC_LINEAR_SOLVER_H

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

Hosted By:
SourceForge.net Logo