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
backendpetsc.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: 2007-05-25
7 
8  Copyright (C) 2007-2011 Université Joseph Fourier (Grenoble I)
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 */
30 #ifndef _BACKENDPETSC_HPP_
31 #define _BACKENDPETSC_HPP_
32 
33 #include <boost/program_options/variables_map.hpp>
34 
39 #include <feel/feelalg/backend.hpp>
40 
41 
42 
43 namespace Feel
44 {
45 namespace po = boost::program_options;
46 
47 
48 
49 
50 
51 #if defined( FEELPP_HAS_PETSC_H )
52 
57 template<typename T>
58 class BackendPetsc : public Backend<T>
59 {
60  typedef Backend<T> super;
61 public:
62 
63  // -- TYPEDEFS --
64  typedef typename super::value_type value_type;
65 
66  /* matrix */
67  typedef typename super::sparse_matrix_type sparse_matrix_type;
68  typedef typename super::sparse_matrix_ptrtype sparse_matrix_ptrtype;
69  typedef MatrixPetsc<value_type> petsc_sparse_matrix_type;
70  typedef boost::shared_ptr<sparse_matrix_type> petsc_sparse_matrix_ptrtype;
71  typedef MatrixPetscMPI<value_type> petscMPI_sparse_matrix_type;
72  //typedef boost::shared_ptr<sparse_matrix_type> petscMPI_sparse_matrix_ptrtype;
73 
74  typedef typename sparse_matrix_type::graph_type graph_type;
75  typedef typename sparse_matrix_type::graph_ptrtype graph_ptrtype;
76 
77  /* vector */
78  typedef typename super::vector_type vector_type;
79  typedef typename super::vector_ptrtype vector_ptrtype;
80  typedef VectorPetsc<value_type> petsc_vector_type;
81  typedef boost::shared_ptr<vector_type> petsc_vector_ptrtype;
82  typedef VectorPetscMPI<value_type> petscMPI_vector_type;
83 
84  typedef typename super::solve_return_type solve_return_type;
85  typedef typename super::nl_solve_return_type nl_solve_return_type;
86 
87  typedef typename super::datamap_type datamap_type;
88  typedef typename super::datamap_ptrtype datamap_ptrtype;
89 
90 
91  // -- CONSTRUCTOR --
92  BackendPetsc( WorldComm const& worldComm=Environment::worldComm() )
93  :
94  super( worldComm ),
95  M_solver_petsc( worldComm ),
96  M_nl_solver_petsc( worldComm )
97  {}
98 
99  BackendPetsc( po::variables_map const& vm, std::string const& prefix = "",
100  WorldComm const& worldComm=Environment::worldComm() )
101  :
102  super( vm, prefix, worldComm ),
103  M_solver_petsc( vm, worldComm ),
104  M_nl_solver_petsc( worldComm )
105  {
106  std::string _prefix = prefix;
107 
108  if ( !_prefix.empty() )
109  _prefix += "-";
110  }
111  ~BackendPetsc();
112  void clear();
113  // -- FACTORY METHODS --
114  template<typename DomainSpace, typename DualImageSpace>
115  static sparse_matrix_ptrtype newMatrix( DomainSpace const& Xh,
116  DualImageSpace const& Yh,
117  size_type matrix_properties = NON_HERMITIAN )
118  {
119  auto s = stencil( _test=Yh,_trial=Xh );
120 
121  sparse_matrix_ptrtype mat;
122  if ( Yh->worldComm().globalSize()>1 )
123  mat = sparse_matrix_ptrtype( new petscMPI_sparse_matrix_type( Yh->dof(),Xh->dof() ) );
124  else // seq
125  mat = sparse_matrix_ptrtype( new petsc_sparse_matrix_type( Yh->dof(),Xh->dof() ) );
126 
127  mat->setMatrixProperties( matrix_properties );
128  mat->init( Yh->nDof(), Xh->nDof(),
129  Yh->nLocalDofWithoutGhost(), Xh->nLocalDofWithoutGhost(),
130  s->graph() );
131  //Yh->nLocalDof(), Xh->nLocalDof() );
132 #if 0
133  auto nSpace = DomainSpace::nSpaces;
134 
135  std::vector < std::vector<int> > is( nSpace );
136  uint cptSpaces=0;
137 
138  //boost::tuple< typename DomainSpace::functionspace_vector_type, uint, std::vector < std::vector<int> > > hola;
139  // auto result = boost::make_tuple(Xh->functionSpaces(),cptSpaces,is);
140  auto result = boost::make_tuple( cptSpaces,is );
141  boost::fusion::fold( Xh->functionSpaces(), result, computeNDofForEachSpace() );
142 
143  for ( uint i = 0; i<nSpace; i++ )
144  {
145  //is[i].resize()
146  }
147 
148 #endif
149 
150  return mat;
151  }
152 
153  sparse_matrix_ptrtype
154  newMatrix( const size_type m,
155  const size_type n,
156  const size_type m_l,
157  const size_type n_l,
158  const size_type nnz=30,
159  const size_type noz=10,
160  size_type matrix_properties = NON_HERMITIAN )
161  {
162  sparse_matrix_ptrtype mat;
163 
164  if ( this->comm().globalSize()>1 )
165  mat = sparse_matrix_ptrtype( new petscMPI_sparse_matrix_type );
166  else
167  mat = sparse_matrix_ptrtype( new petsc_sparse_matrix_type );
168 
169  mat->setMatrixProperties( matrix_properties );
170  mat->init( m,n,m_l,n_l,nnz,noz );
171  return mat;
172  }
173 
174  sparse_matrix_ptrtype
175  newMatrix( datamap_ptrtype const& domainmap,
176  datamap_ptrtype const& imagemap,
177  size_type matrix_properties = NON_HERMITIAN,
178  bool init = true )
179  {
180  sparse_matrix_ptrtype mat;
181 
182  if ( imagemap->worldComm().globalSize()>1 )
183  mat = sparse_matrix_ptrtype( new petscMPI_sparse_matrix_type( imagemap,domainmap,imagemap->worldComm() ) );
184  else
185  mat = sparse_matrix_ptrtype( new petsc_sparse_matrix_type( imagemap,domainmap,imagemap->worldComm() ) );
186 
187  mat->setMatrixProperties( matrix_properties );
188 
189  if ( init )
190  {
191  mat->init( imagemap->nDof(), domainmap->nDof(),
192  imagemap->nLocalDofWithoutGhost(), domainmap->nLocalDofWithoutGhost() );
193  }
194 
195  return mat;
196  }
197 
198  sparse_matrix_ptrtype
199  newMatrix( const size_type m,
200  const size_type n,
201  const size_type m_l,
202  const size_type n_l,
203  graph_ptrtype const & graph,
204  size_type matrix_properties = NON_HERMITIAN )
205  {
206  sparse_matrix_ptrtype mat;
207  auto const& mapGraphRow = graph->mapRowPtr();
208  auto const& mapGraphCol = graph->mapColPtr();
209  if ( this->comm().globalSize()>1 )
210  mat = sparse_matrix_ptrtype( new petscMPI_sparse_matrix_type( mapGraphRow,mapGraphCol,mapGraphRow->worldComm() ) );
211  else
212  mat = sparse_matrix_ptrtype( new petsc_sparse_matrix_type( mapGraphRow,mapGraphCol,mapGraphRow->worldComm() ) ) ;
213 
214  mat->setMatrixProperties( matrix_properties );
215  //mat->init( m,n,m_l,n_l,graph );
216  mat->init( mapGraphRow->nDof(), mapGraphCol->nDof(),
217  mapGraphRow->nLocalDofWithoutGhost(), mapGraphCol->nLocalDofWithoutGhost(),
218  graph );
219 
220 
221  return mat;
222  }
223 
224  sparse_matrix_ptrtype
225  newZeroMatrix( datamap_ptrtype const& domainmap,
226  datamap_ptrtype const& imagemap )
227  {
228  graph_ptrtype sparsity_graph( new graph_type( imagemap, domainmap ) );
229  sparsity_graph->zero();
230  sparsity_graph->close();
231 
232  sparse_matrix_ptrtype mat;
233  if ( imagemap->worldComm().globalSize()>1 )
234  mat = sparse_matrix_ptrtype( new petscMPI_sparse_matrix_type( imagemap,domainmap,imagemap->worldComm() ) );
235  else
236  mat = sparse_matrix_ptrtype( new petsc_sparse_matrix_type( imagemap,domainmap,imagemap->worldComm() ) );
237 
238  mat->init( imagemap->nDof(), domainmap->nDof(),
239  imagemap->nLocalDofWithoutGhost(), domainmap->nLocalDofWithoutGhost(),
240  sparsity_graph );
241 
242  return mat;
243  }
244 
245  sparse_matrix_ptrtype
246  newZeroMatrix( const size_type m,
247  const size_type n,
248  const size_type m_l,
249  const size_type n_l )
250  {
251  graph_ptrtype sparsity_graph( new graph_type( 0,0,m_l-1,0,n_l-1 ) );
252  sparsity_graph->zero();
253  sparsity_graph->close();
254 
255  sparse_matrix_ptrtype mat;
256 
257  if ( this->comm().globalSize()>1 )
258  mat = sparse_matrix_ptrtype( new petscMPI_sparse_matrix_type );
259  else
260  mat = sparse_matrix_ptrtype( new petsc_sparse_matrix_type );
261 
262  //mat->setMatrixProperties( matrix_properties );
263  mat->init( m, n, m_l, n_l, sparsity_graph );
264 
265  return mat;
266  }
267 
268  template<typename SpaceT>
269  static vector_ptrtype newVector( SpaceT const& space )
270  {
271  if ( space->worldComm().globalSize()>1 )
272  return vector_ptrtype( new petscMPI_vector_type( space->dof() ) );
273  else
274  return vector_ptrtype( new petsc_vector_type( space->dof() ) );
275  }
276 
277  vector_ptrtype newVector( datamap_ptrtype const& dm )
278  {
279  if ( this->comm().globalSize()>1 ) return vector_ptrtype( new petscMPI_vector_type( dm ) );
280 
281  else return vector_ptrtype( new petsc_vector_type( dm ) );
282  }
283 
284  vector_ptrtype newVector( const size_type n, const size_type n_local )
285  {
286  return vector_ptrtype( new petsc_vector_type( n, n_local ) );
287  }
288 
289 
290  void set_symmetric( bool /*is_sym*/ ) {}
291 
292  // -- LINEAR ALGEBRA INTERFACE --
293  template <class Vector>
294  static void applyMatrix( sparse_matrix_type const& A,
295  const Vector& x,
296  vector_type& b )
297  {
298  // petsc mat/vec here
299  }
300 
301  void prod( sparse_matrix_type const& A,
302  vector_type const& x,
303  vector_type& b ) const
304  {
305  int ierr = 0;
306  petsc_sparse_matrix_type const& _A = dynamic_cast<petsc_sparse_matrix_type const&>( A );
307  petsc_vector_type const& _x = dynamic_cast<petsc_vector_type const&>( x );
308  petsc_vector_type const& _b = dynamic_cast<petsc_vector_type const&>( b );
309  if ( _A.mapCol().worldComm().globalSize() == x.map().worldComm().globalSize() )
310  {
311  //std::cout << "BackendPetsc::prod STANDART"<< std::endl;
312  ierr = MatMult( _A.mat(), _x.vec(), _b.vec() );
313  CHKERRABORT( _A.comm().globalComm(),ierr );
314  }
315  else
316  {
317  //std::cout << "BackendPetsc::prod with convert"<< std::endl;
318  auto x_convert = petscMPI_vector_type(_A.mapColPtr());
319  x_convert.duplicateFromOtherPartition(x);
320  x_convert.close();
321  ierr = MatMult( _A.mat(), x_convert.vec(), _b.vec() );
322  CHKERRABORT( _A.comm().globalComm(),ierr );
323  }
324  b.close();
325  }
326 
327 
328  solve_return_type solve( sparse_matrix_type const& A,
329  vector_type& x,
330  vector_type const& b );
331 
332  solve_return_type solve( sparse_matrix_ptrtype const& A,
333  sparse_matrix_ptrtype const& B,
334  vector_ptrtype& x,
335  vector_ptrtype const& b );
336 
337  template <class Vector>
338  static value_type dot( const vector_type& f,
339  const Vector& x )
340  {
341  value_type result( 0 );
342 
343  // petsc dot here
344 
345  return result;
346  }
347 
348 private:
349 
350  SolverLinearPetsc<double> M_solver_petsc;
351  SolverNonLinearPetsc<double> M_nl_solver_petsc;
352 
353 }; // class BackendPetsc
354 
355 
356 template<typename T>
357 BackendPetsc<T>::~BackendPetsc()
358 {
359  this->clear();
360 }
361 template<typename T>
362 void
363 BackendPetsc<T>::clear()
364 {
365  LOG(INFO) << "Deleting linear solver petsc";
366  M_solver_petsc.clear();
367  LOG(INFO) << "Deleting non linear solver petsc";
368  M_nl_solver_petsc.clear();
369  LOG(INFO) << "Deleting backend petsc";
370 
371  super::clear();
372 
373 }
374 
375 template<typename T>
376 typename BackendPetsc<T>::solve_return_type
377 BackendPetsc<T>::solve( sparse_matrix_ptrtype const& A,
378  sparse_matrix_ptrtype const& B,
379  vector_ptrtype& x,
380  vector_ptrtype const& b )
381 {
382  M_solver_petsc.setPrefix( this->prefix() );
383  M_solver_petsc.setPreconditionerType( this->pcEnumType() );
384  M_solver_petsc.setSolverType( this->kspEnumType() );
385  if (!M_solver_petsc.initialized())
386  M_solver_petsc.attachPreconditioner( this->M_preconditioner );
387  M_solver_petsc.setConstantNullSpace( this->hasConstantNullSpace() );
388  M_solver_petsc.setFieldSplitType( this->fieldSplitEnumType() );
389  M_solver_petsc.setTolerances( _rtolerance=this->rTolerance(),
390  _atolerance=this->aTolerance(),
391  _dtolerance=this->dTolerance(),
392  _maxit = this->maxIterations() );
393  M_solver_petsc.setPrecMatrixStructure( this->precMatrixStructure() );
394  M_solver_petsc.setMatSolverPackageType( this->matSolverPackageEnumType() );
395  M_solver_petsc.setShowKSPMonitor( this->showKSPMonitor() );
396  M_solver_petsc.setShowKSPConvergedReason( this->showKSPConvergedReason() );
397 
398  auto res = M_solver_petsc.solve( *A, *B, *x, *b, this->rTolerance(), this->maxIterations(), this->transpose() );
399  DVLOG(2) << "[BackendPetsc::solve] number of iterations : " << res.template get<1>() << "\n";
400  DVLOG(2) << "[BackendPetsc::solve] residual : " << res.template get<2>() << "\n";
401 
402  if ( !res.template get<0>() )
403  LOG(ERROR) << "Backend " << this->prefix() << " : linear solver failed to converge" << std::endl;
404 
405  return res;
406 } // BackendPetsc::solve
407 
408 
409 template<typename T>
410 typename BackendPetsc<T>::solve_return_type
411 BackendPetsc<T>::solve( sparse_matrix_type const& A,
412  vector_type& x,
413  vector_type const& b )
414 {
415  M_solver_petsc.setPrefix( this->prefix() );
416  M_solver_petsc.setPreconditionerType( this->pcEnumType() );
417  M_solver_petsc.setSolverType( this->kspEnumType() );
418  if (!M_solver_petsc.initialized())
419  M_solver_petsc.attachPreconditioner( this->M_preconditioner );
420  M_solver_petsc.setConstantNullSpace( this->hasConstantNullSpace() );
421  M_solver_petsc.setFieldSplitType( this->fieldSplitEnumType() );
422  M_solver_petsc.setTolerances( _rtolerance=this->rTolerance(),
423  _atolerance=this->aTolerance(),
424  _dtolerance=this->dTolerance(),
425  _maxit = this->maxIterations() );
426  M_solver_petsc.setPrecMatrixStructure( this->precMatrixStructure() );
427  M_solver_petsc.setMatSolverPackageType( this->matSolverPackageEnumType() );
428  M_solver_petsc.setShowKSPMonitor( this->showKSPMonitor() );
429  M_solver_petsc.setShowKSPConvergedReason( this->showKSPConvergedReason() );
430 
431  auto res = M_solver_petsc.solve( A, x, b, this->rTolerance(), this->maxIterations() );
432  DVLOG(2) << "[BackendPetsc::solve] number of iterations : " << res.template get<1>() << "\n";
433  DVLOG(2) << "[BackendPetsc::solve] residual : " << res.template get<2>() << "\n";
434 
435  if ( !res.template get<0>() )
436  LOG(ERROR) << "Backend " << this->prefix() << " : linear solver failed to converge" << std::endl;
437 
438  return res;
439 } // BackendPetsc::solve
440 
441 
442 
443 #endif // FEELPP_HAS_PETSC_H
444 } // Feel
445 
446 #endif /* _BACKENDPETSC_HPP_ */

Generated on Sun Oct 20 2013 08:24:54 for Feel++ by doxygen 1.8.4