Logo  0.95.0-final
Finite Element Embedded Library and Language in C++
Feel++ Feel++ on Github Feel++ community
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
solverlinearpetsc.hpp
Go to the documentation of this file.
1 /* -*- mode: c++; coding: utf-8; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; show-trailing-whitespace: t -*- vim:fenc=utf-8:ft=tcl:et:sw=4:ts=4:sts=4
2 
3  This file is part of the Feel library
4 
5  Author(s): Christophe Prud'homme <christophe.prudhomme@feelpp.org>
6  Date: 2005-11-27
7 
8  Copyright (C) 2005,2006 EPFL
9 
10  This library is free software; you can redistribute it and/or
11  modify it under the terms of the GNU Lesser General Public
12  License as published by the Free Software Foundation; either
13  version 3.0 of the License, or (at your option) any later version.
14 
15  This library is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18  Lesser General Public License for more details.
19 
20  You should have received a copy of the GNU Lesser General Public
21  License along with this library; if not, write to the Free Software
22  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
23 */
24 // $Id: petsc_linear_solver.h,v 1.2 2005/02/22 22:17:34 jwpeterson Exp $
25 
26 // The libMesh Finite Element Library.
27 // Copyright (C) 2002-2005 Benjamin S. Kirk, John W. Peterson
28 
29 // This library is free software; you can redistribute it and/or
30 // modify it under the terms of the GNU Lesser General Public
31 // License as published by the Free Software Foundation; either
32 // version 3.0 of the License, or (at your option) any later version.
33 
34 // This library is distributed in the hope that it will be useful,
35 // but WITHOUT ANY WARRANTY; without even the implied warranty of
36 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
37 // Lesser General Public License for more details.
38 
39 // You should have received a copy of the GNU Lesser General Public
40 // License along with this library; if not, write to the Free Software
41 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
48 #ifndef __petsc_linear_solver_h__
49 #define __petsc_linear_solver_h__
50 
51 
52 #include <feel/feelcore/feel.hpp>
53 
55 #include <feel/feelalg/solverlinear.hpp>
56 
60 
61 
65 #ifdef FEELPP_HAS_PETSC_H
67 
68 #ifndef USE_COMPLEX_NUMBERS
69 extern "C" {
70 # include <petscversion.h>
71 # if (PETSC_VERSION_MAJOR == 2) && (PETSC_VERSION_MINOR <= 1)
72 # include <petscsles.h>
73 # else
74 # include <petscksp.h>
75 # endif
76 }
77 #else
78 # include <petscversion.h>
79 # if (PETSC_VERSION_MAJOR == 2) && (PETSC_VERSION_MINOR <= 1)
80 # include <petscsles.h>
81 # else
82 # include <petscksp.h>
83 # endif
84 #endif
85 
86 //--------------------------------------------------------------------
87 // Functions with C linkage to pass to PETSc. PETSc will call these
88 // methods as needed for preconditioning
89 //
90 // Since they must have C linkage they have no knowledge of a namespace.
91 // Give them an obscure name to avoid namespace pollution.
92 extern "C"
93 {
94 #if PETSC_VERSION_LESS_THAN(3,0,1) && PETSC_VERSION_RELEASE
95 
99  PetscErrorCode __feel_petsc_preconditioner_setup ( void * ctx );
100 
105  PetscErrorCode __feel_petsc_preconditioner_apply( void *ctx, Vec x, Vec y );
106 #else
107  PetscErrorCode __feel_petsc_preconditioner_setup ( PC );
108  PetscErrorCode __feel_petsc_preconditioner_apply( PC, Vec x, Vec y );
109 #endif
110 } // end extern "C"
111 
112 
113 
114 namespace Feel
115 {
116 
117 
127 template <typename T>
128 class SolverLinearPetsc : public SolverLinear<T>
129 {
130  typedef SolverLinear<T> super;
131 
132 public:
133 
134  typedef typename super::value_type value_type;
135  typedef typename super::real_type real_type;
136 
140  SolverLinearPetsc ( WorldComm const& worldComm=Environment::worldComm() );
141 
145  SolverLinearPetsc ( po::variables_map const& vm, WorldComm const& worldComm=Environment::worldComm() );
146 
150  ~SolverLinearPetsc ();
151 
155  void clear ();
156 
160  void init ();
161 
165  void setConstantNullSpace( bool cns )
166  {
167  M_constant_null_space = cns;
168  }
169 
182  //std::pair<unsigned int, real_type>
183  boost::tuple<bool,unsigned int, real_type>
184  solve ( MatrixSparse<T> const &mat,
185  Vector<T> & x,
186  Vector<T> const& b,
187  const double tolerance,
188  const unsigned int maxit,
189  bool transpose )
190  {
191  return this->solve( mat, mat, x, b, tolerance, maxit, transpose );
192  }
193 
194  boost::tuple<bool,unsigned int, real_type>
195  solve ( MatrixShell<T> const &mat,
196  Vector<T> & x,
197  Vector<T> const& b,
198  const double tolerance,
199  const unsigned int maxit,
200  bool transpose );
201 
227  //std::pair<unsigned int, real_type>
228  boost::tuple<bool,unsigned int, real_type>
229  solve ( MatrixSparse<T> const& mat,
230  MatrixSparse<T> const& prec,
231  Vector<T> & x,
232  Vector<T> const& b,
233  const double tolerance,
234  const unsigned int maxit,
235  bool transpose );
236 
242  PC pc()
243  {
244  this->init();
245  return M_pc;
246  }
247 
252  void getResidualHistory( std::vector<double>& hist );
253 
260  real_type getInitialResidual();
261 
262 private:
263 
268  void setPetscSolverType ();
269 
274  void setPetscPreconditionerType ();
275 
276  void setPetscConstantNullSpace ();
277 
278  // SLES removed from >= PETSc 2.2.0
279 #if (PETSC_VERSION_MAJOR == 2) && (PETSC_VERSION_MINOR <= 1)
280 
284  SLES M_sles;
285 
286 #endif
287 
291  PC M_pc;
292 
296  KSP M_ksp;
297 
298  bool M_constant_null_space;
299 };
300 
301 
302 /*----------------------- functions ----------------------------------*/
303 template <typename T>
304 inline
305 SolverLinearPetsc<T>::SolverLinearPetsc ( WorldComm const& worldComm )
306  :
307  super( worldComm ),
308  M_constant_null_space( false )
309 {
310  if ( this->worldComm().globalSize() == 1 )
311  this->setPreconditionerType( LU_PRECOND );
312 
313  else
314  this->setPreconditionerType( BLOCK_JACOBI_PRECOND );
315 }
316 
317 template <typename T>
318 inline
319 SolverLinearPetsc<T>::SolverLinearPetsc ( po::variables_map const& vm, WorldComm const& worldComm )
320  :
321  super( vm, worldComm ),
322  M_constant_null_space( false )
323 {
324 }
325 
326 
327 
328 template <typename T>
329 inline
330 SolverLinearPetsc<T>::~SolverLinearPetsc ()
331 {
332  this->clear ();
333 }
334 } // Feel
335 
336 
337 #endif // #ifdef FEELPP_HAS_PETSC_H
338 #endif // #ifdef __petsc_linear_solver_h__

Generated on Sun Oct 20 2013 08:25:05 for Feel++ by doxygen 1.8.4