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
backend.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-12-23
7 
8  Copyright (C) 2007-2012 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 */
29 #ifndef __Backend_H
30 #define __Backend_H 1
31 
32 #include <boost/timer.hpp>
33 #include <boost/tuple/tuple.hpp>
34 #include <boost/fusion/include/fold.hpp>
35 #include <boost/smart_ptr/enable_shared_from_this.hpp>
36 
37 #include <feel/feelcore/feel.hpp>
40 #include <feel/feelcore/parameter.hpp>
41 #include <feel/feelalg/enums.hpp>
42 #include <feel/feelalg/vector.hpp>
45 #include <feel/feelalg/vectorblock.hpp>
46 #include <feel/feelalg/datamap.hpp>
47 
51 
54 //#include <feel/feelvf/vf.hpp>
55 //#include <boost/fusion/support/pair.hpp>
56 //#include <boost/fusion/container.hpp>
57 //#include <boost/fusion/sequence.hpp>
58 //#include <boost/fusion/algorithm.hpp>
59 
60 //namespace fusion = boost::fusion;
61 
62 //#include <feel/feelvf/bilinearform.hpp>
63 #include <feel/feelvf/pattern.hpp>
64 #include <feel/feelvf/block.hpp>
65 
66 namespace Feel
67 {
68 /*enum Pattern
69 {
70 DEFAULT = 1 << 0,
71 EXTENDED = 1 << 1,
72 COUPLED = 1 << 2,
73 SYMMETRIC = 1 << 3
74 };
75 */
77 namespace detail
78 {
79 template<typename T>
80 boost::shared_ptr<DataMap> datamap( T const& t, mpl::true_ )
81 {
82  return t->mapPtr();
83 }
84 template<typename T>
85 boost::shared_ptr<DataMap> datamap( T const& t, mpl::false_ )
86 {
87  return t.mapPtr();
88 }
89 template<typename T>
90 boost::shared_ptr<DataMap> datamap( T const& t )
91 {
92  return datamap( t, detail::is_shared_ptr<T>() );
93 }
94 
95 template<typename T>
96 #if BOOST_VERSION >= 105300
97 typename boost::detail::sp_dereference< typename T::element_type >::type
98 #else
99 typename T::reference
100 #endif
101 ref( T t, mpl::true_ )
102 {
103  return *t;
104 }
105 template<typename T>
106 T& ref( T& t, mpl::false_ )
107 {
108  return t;
109 }
110 template<typename T>
111 auto ref( T& t ) -> decltype( ref( t, detail::is_shared_ptr<T>() ) )
112 {
113  return ref( t, detail::is_shared_ptr<T>() );
114 }
115 
116 
117 }
119 
120 template<typename T> class MatrixBlockBase;
121 template<int NR, int NC, typename T> class MatrixBlock;
122 template<typename T> class VectorBlockBase;
123 template<int NR, typename T> class VectorBlock;
124 
132 template<typename T>
133 class Backend:
134  public boost::enable_shared_from_this<Backend<T> >
135 {
136 public:
137 
138 
142  typedef T value_type;
143  typedef typename type_traits<T>::real_type real_type;
144 
146  typedef boost::shared_ptr<vector_type> vector_ptrtype;
148  typedef boost::shared_ptr<sparse_matrix_type> sparse_matrix_ptrtype;
149 
151  typedef boost::shared_ptr<shell_matrix_type> shell_matrix_ptrtype;
152 
154  typedef typename sparse_matrix_type::graph_ptrtype graph_ptrtype;
155 
157  typedef boost::shared_ptr<backend_type> backend_ptrtype;
158 
160  typedef boost::shared_ptr<solvernonlinear_type> solvernonlinear_ptrtype;
161 
162  typedef boost::tuple<bool, size_type, value_type> solve_return_type;
163  typedef boost::tuple<bool, size_type, value_type> nl_solve_return_type;
164 
165  typedef DataMap datamap_type;
166  typedef boost::shared_ptr<datamap_type> datamap_ptrtype;
167 
169 
173 
174  Backend( WorldComm const& worldComm=Environment::worldComm() );
175  Backend( po::variables_map const& vm, std::string const& prefix = "", WorldComm const& worldComm=Environment::worldComm() );
176  Backend( Backend const & );
177  virtual ~Backend();
178 
179 
184  static backend_ptrtype build(
185 #if defined( FEELPP_HAS_PETSC_H )
186  BackendType = BACKEND_PETSC
187 #else
188  BackendType = BACKEND_GMM
189 #endif
190  , WorldComm const& worldComm=Environment::worldComm()
191  );
192 
196  static backend_ptrtype build( po::variables_map const& vm, std::string const& prefix = "", WorldComm const& worldComm=Environment::worldComm() );
197 
201  virtual sparse_matrix_ptrtype newMatrix( const size_type m,
202  const size_type n,
203  const size_type m_l,
204  const size_type n_l,
205  const size_type nnz=30,
206  const size_type noz=10,
207  size_type prop = NON_HERMITIAN ) = 0;
208 
212  virtual sparse_matrix_ptrtype newMatrix( const size_type m,
213  const size_type n,
214  const size_type m_l,
215  const size_type n_l,
216  graph_ptrtype const & graph,
217  size_type matrix_properties = NON_HERMITIAN ) = 0;
218 
222  sparse_matrix_ptrtype newMatrix( const size_type m,
223  const size_type n,
224  const size_type m_l,
225  const size_type n_l,
226  graph_ptrtype const & graph,
227  std::vector < std::vector<size_type> > indexSplit,
228  size_type matrix_properties = NON_HERMITIAN )
229  {
230  auto mat = this->newMatrix( m,n,m_l,n_l,graph,matrix_properties );
231  mat->setIndexSplit( indexSplit );
232  return mat;
233  }
234 
235 
239  virtual sparse_matrix_ptrtype newMatrix( datamap_ptrtype const& dm1,
240  datamap_ptrtype const& dm2,
241  size_type prop = NON_HERMITIAN,
242  bool init = true ) = 0;
243 
247  sparse_matrix_ptrtype newMatrix( datamap_ptrtype const& domainmap,
248  datamap_ptrtype const& imagemap,
249  graph_ptrtype const & graph,
250  size_type matrix_properties = NON_HERMITIAN,
251  bool init = true )
252  {
253  auto mat = this->newMatrix( domainmap,imagemap, matrix_properties, false );
254 
255  if ( init ) mat->init( imagemap->nDof(), domainmap->nDof(),
256  imagemap->nLocalDofWithoutGhost(), domainmap->nLocalDofWithoutGhost(),
257  graph );
258 
259  mat->zero();
260  // todo!
261  //mat->setIndexSplit( trial->dofIndexSplit() );
262  return mat;
263  }
264 
265 
269  virtual
270  sparse_matrix_ptrtype
271  newZeroMatrix( const size_type m,
272  const size_type n,
273  const size_type m_l,
274  const size_type n_l ) =0;
275 
276  virtual sparse_matrix_ptrtype newZeroMatrix( datamap_ptrtype const& dm1, datamap_ptrtype const& dm2 ) = 0;
277 
281  BOOST_PARAMETER_MEMBER_FUNCTION( ( sparse_matrix_ptrtype ),
282  newMatrix,
283  tag,
284  ( required
285  ( trial,*( boost::is_convertible<mpl::_,boost::shared_ptr<FunctionSpaceBase> > ) )
286  ( test,*( boost::is_convertible<mpl::_,boost::shared_ptr<FunctionSpaceBase> > ) ) )
287  ( optional
288  ( pattern,( size_type ),Pattern::COUPLED )
289  ( properties,( size_type ),NON_HERMITIAN )
290  ( buildGraphWithTranspose, ( bool ),false )
291  ( pattern_block, *, ( BlocksStencilPattern(1,1,size_type( Pattern::HAS_NO_BLOCK_PATTERN ) ) ) )
292  ( diag_is_nonzero, *( boost::is_integral<mpl::_> ), true )
293  ( verbose,( int ),0 )
294  ( collect_garbage, *( boost::is_integral<mpl::_> ), true )
295  ) )
296  {
297 
298  //auto mat = this->newMatrix( trial->map(), test->map(), properties, false );
299  auto mat = this->newMatrix( trial->dofOnOff(), test->dofOn(), properties, false );
300 
301  if ( !buildGraphWithTranspose )
302  {
303  auto s = stencil( _test=test,
304  _trial=trial,
305  _pattern=pattern,
306  _pattern_block=pattern_block,
307  _diag_is_nonzero=diag_is_nonzero,
308  _collect_garbage=collect_garbage);
309 
310  mat->init( test->nDof(), trial->nDof(),
311  test->nLocalDofWithoutGhost(), trial->nLocalDofWithoutGhost(),
312  s->graph() );
313  }
314  else
315  {
316  auto s = stencil( _test=trial,
317  _trial=test,
318  _pattern=pattern,
319  _pattern_block=pattern_block.transpose(),
320  _diag_is_nonzero=false,// because transpose(do just after)
321  _close=false,
322  _collect_garbage=collect_garbage );
323  // get the good graph
324  auto graph = s->graph()->transpose(false);
325  if ( diag_is_nonzero )
326  graph->addMissingZeroEntriesDiagonal();
327  graph->close();
328 
329  //maybe do that
330  //stencilManagerGarbage(boost::make_tuple( trial, test, pattern, pattern_block.transpose().getSetOfBlocks(), false/*diag_is_nonzero*/));
331  //now save the good graph in the StencilManager(only if entry is empty)
332  stencilManagerAdd(boost::make_tuple( test, trial, pattern, pattern_block.getSetOfBlocks(), diag_is_nonzero), graph);
333 
334  mat->init( test->nDof(), trial->nDof(),
335  test->nLocalDofWithoutGhost(), trial->nLocalDofWithoutGhost(),
336  graph );
337  }
338 
339  mat->zero();
340  mat->setIndexSplit( trial->dofIndexSplit() );
341  return mat;
342  }
343 
344  template<typename DomainSpace, typename ImageSpace>
345  sparse_matrix_ptrtype newMatrix( DomainSpace const& dm, ImageSpace const& im, sparse_matrix_ptrtype const& M, size_type prop = NON_HERMITIAN )
346  {
347  sparse_matrix_ptrtype m = newMatrix( dm, im, prop );
348  m->init( im->nDof(), dm->nDof(), im->nLocalDof(), dm->nLocalDof(), M->graph() );
349  return m;
350  }
351 
355  sparse_matrix_ptrtype newBlockMatrixImpl( vf::BlocksBase<sparse_matrix_ptrtype> const & b,
356  bool copy_values=true,
357  bool diag_is_nonzero=true )
358  {
360  boost::shared_ptr<matrix_block_type> mb( new matrix_block_type( b, *this, copy_values, diag_is_nonzero ) );
361  return mb->getSparseMatrix();
362  }
363 
364  sparse_matrix_ptrtype newBlockMatrixImpl( vf::BlocksBase<boost::shared_ptr<GraphCSR> > const & b,
365  bool copy_values=true,
366  bool diag_is_nonzero=true )
367  {
368  typedef MatrixBlockBase<value_type> matrix_block_type;
369  boost::shared_ptr<matrix_block_type> mb( new matrix_block_type( b, *this, diag_is_nonzero ) );
370  return mb->getSparseMatrix();
371  }
372 
376  BOOST_PARAMETER_MEMBER_FUNCTION( ( sparse_matrix_ptrtype ),
377  newBlockMatrix,
378  tag,
379  ( required
380  ( block,* )
381  )
382  ( optional
383  ( copy_values,*( boost::is_integral<mpl::_> ),true )
384  ( diag_is_nonzero, *( boost::is_integral<mpl::_> ), true )
385  )
386  )
387  {
388  return newBlockMatrixImpl( block,copy_values,diag_is_nonzero );
389  }
390 
394  template < typename BlockType=vector_ptrtype >
395  vector_ptrtype newBlockVectorImpl( vf::BlocksBase<BlockType> const & b,
396  bool copy_values=true )
397  {
398  typedef VectorBlockBase<typename BlockType::element_type::value_type> vector_block_type;
399  boost::shared_ptr<vector_block_type> mb( new vector_block_type( b, *this, copy_values ) );
400  return mb->getVector();
401  }
402 
406  BOOST_PARAMETER_MEMBER_FUNCTION( ( vector_ptrtype ),
407  newBlockVector,
408  tag,
409  ( required
410  ( block,* )
411  )
412  ( optional
413  ( copy_values,*( boost::is_integral<mpl::_> ),true )
414  )
415  )
416  {
417  return newBlockVectorImpl( block,copy_values );
418  }
419 
420 
424  BOOST_PARAMETER_MEMBER_FUNCTION( ( sparse_matrix_ptrtype ),
426  tag,
427  ( required
428  ( test,*( boost::is_convertible<mpl::_,boost::shared_ptr<FunctionSpaceBase> >) )
429  ( trial,*( boost::is_convertible<mpl::_,boost::shared_ptr<FunctionSpaceBase> >) )
430  )
431  )
432  {
433  return this->newZeroMatrix( trial->dofOnOff(), test->dofOn() );
434  }
435 
439  virtual vector_ptrtype newVector( datamap_ptrtype const& dm ) = 0;
440 
444  virtual vector_ptrtype newVector( const size_type n, const size_type n_local ) = 0;
445 
449  BOOST_PARAMETER_MEMBER_FUNCTION( ( vector_ptrtype ),
450  newVector,
451  tag,
452  ( required
453  ( test,*( boost::is_convertible<mpl::_,boost::shared_ptr<FunctionSpaceBase> >) )
454  )
455  )
456  {
457  return this->newVector( test->dof() );
458  }
459 
461 
465 
466 
468 
472 
476  std::string prefix() const
477  {
478  return M_prefix;
479  }
480 
484  std::string kspType() const
485  {
486  return M_ksp;
487  }
488 
492  std::string pcType() const
493  {
494  return M_pc;
495  }
496 
500  bool hasConstantNullSpace() const
501  {
502  return M_constant_null_space;
503  }
504 
508  std::string fieldsplitType() const
509  {
510  return M_fieldSplit;
511  }
512 
517 
521  SolverType kspEnumType() const;
522 
527 
531  std::string pcFactorMatSolverPackageType() const
532  {
533  return M_pcFactorMatSolverPackage;
534  }
535 
539  MatSolverPackageType matSolverPackageEnumType() const;
540 
545  {
546  return M_prec_matrix_structure;
547  }
548 
552  value_type rTolerance() const
553  {
554  return M_rtolerance;
555  }
556 
560  value_type rToleranceSNES() const
561  {
562  return M_rtoleranceSNES;
563  }
564 
568  value_type dTolerance() const
569  {
570  return M_dtolerance;
571  }
572 
576  value_type sToleranceSNES() const
577  {
578  return M_stoleranceSNES;
579  }
580 
584  value_type aTolerance() const
585  {
586  return M_atolerance;
587  }
588 
592  value_type aToleranceSNES() const
593  {
594  return M_atoleranceSNES;
595  }
596 
601  {
602  return M_maxit;
603  }
604 
609  {
610  return M_maxitSNES;
611  }
612 
616  value_type rtoleranceKSPinSNES() const
617  {
618  return M_rtoleranceKSPinSNES;
619  }
620 
621 
622  bool converged() const
623  {
624  return M_converged;
625  }
626 
627  size_type nIterations() const
628  {
629  return M_iteration;
630  }
631 
632  bool transpose() const
633  {
634  return M_transpose;
635  }
636 
640  WorldComm const& comm() const
641  {
642  return M_worldComm;
643  }
644 
645  bool showKSPMonitor() const { return M_showKSPMonitor; }
646  bool showSNESMonitor() const { return M_showSNESMonitor; }
647  bool showKSPConvergedReason() const { return M_showKSPConvergedReason; }
648  bool showSNESConvergedReason() const { return M_showSNESConvergedReason; }
649 
650  bool reusePrec() const { return M_reuse_prec; }
651  bool reuseJac() const { return M_reuse_jac; }
652 
653  bool reusePrecRebuildAtFirstNewtonStep() const { return M_reusePrecRebuildAtFirstNewtonStep; }
654  bool reuseJacRebuildAtFirstNewtonStep() const { return M_reuseJacRebuildAtFirstNewtonStep; }
655 
656  BackendType type() const { return M_backend; }
657 
659 
663 
668  BOOST_PARAMETER_MEMBER_FUNCTION( ( void ),
669  setTolerances,
670  tag,
671  ( required
672  ( rtolerance, ( double ) )
673  )
674  ( optional
675  ( maxit, ( size_type ), 1000 )
676  ( atolerance, ( double ), 1e-50 )
677  ( dtolerance, ( double ), 1e5 )
678  ) )
679  {
680  M_rtolerance = rtolerance;
681  M_dtolerance = dtolerance;
682  M_atolerance = atolerance;
683  M_maxit = maxit;
684  }
685 
686  BOOST_PARAMETER_MEMBER_FUNCTION( ( void ),
687  setTolerancesSNES,
688  tag,
689  ( required
690  ( rtolerance, ( double ) )
691  )
692  ( optional
693  ( maxit, ( size_type ), 50 )
694  ( atolerance, ( double ), 1e-50 )
695  ( stolerance, ( double ), 1e-8 )
696  ) )
697  {
698  M_rtoleranceSNES = rtolerance;
699  M_stoleranceSNES = stolerance;
700  M_atoleranceSNES = atolerance;
701  M_maxitSNES = maxit;
702  }
703 
707  BOOST_PARAMETER_MEMBER_FUNCTION( ( void ),
708  setSolverType,
709  tag,
710  ( required
711  ( ksp, ( std::string ) )
712  )
713  ( optional
714  ( pc, ( std::string ), "lu" )
715  ( constant_null_space, ( bool ), false )
716  ( pcfactormatsolverpackage, ( std::string ), "petsc" )
717  ) )
718  {
719  M_ksp = ksp;
720  M_pc = pc;
721  M_pcFactorMatSolverPackage = pcfactormatsolverpackage;
722  M_constant_null_space = constant_null_space;
723  }
724 
729  {
730  M_prec_matrix_structure = mstruct;
731  }
732 
736  solvernonlinear_ptrtype nlSolver()
737  {
738  return M_nlsolver;
739  }
740 
741  void setTranspose( bool transpose )
742  {
743  M_transpose = transpose;
744  }
745 
746  void setShowKSPMonitor( bool b ) { M_showKSPMonitor=b; }
747  void setShowSNESMonitor( bool b ) { M_showSNESMonitor=b; }
748  void setShowKSPConvergedReason( bool b ) { M_showKSPConvergedReason=b; }
749  void setShowSNESConvergedReason( bool b ) { M_showSNESConvergedReason=b; }
750 
751  void setReusePrec( bool b ) { M_reuse_prec=b; }
752  void setReuseJac( bool b) { M_reuse_jac=b; }
753 
754  void setReusePrecRebuildAtFirstNewtonStep(bool b) { M_reusePrecRebuildAtFirstNewtonStep=b; }
755  void setReuseJacRebuildAtFirstNewtonStep(bool b) { M_reuseJacRebuildAtFirstNewtonStep=b; }
756 
757 
759 
763 
767  virtual void clear();
768 
772  virtual real_type dot( vector_type const& x, vector_type const& y ) const;
773 
774 
778  real_type dot( vector_ptrtype const& x, vector_ptrtype const& y ) const
779  {
780  return this->dot( *x, *y );
781  }
785  virtual void prod( sparse_matrix_type const& A, vector_type const& x, vector_type& y ) const = 0;
786 
790  void prod( sparse_matrix_ptrtype const& A, vector_ptrtype const& x, vector_ptrtype& y ) const
791  {
792  this->prod( *A, *x, *y );
793  }
794 
811  BOOST_PARAMETER_MEMBER_FUNCTION( ( solve_return_type ),
812  solve,
813  tag,
814  ( required
815  //( matrix,*(mpl::or_<sparse_matrix_ptrtype,shell_matrix_ptrtype>) )
816  ( matrix,(sparse_matrix_ptrtype) )
817  //( in_out( solution ),*( mpl::or_<mpl::or_<boost::is_convertible<mpl::_,vector_type&>,boost::is_convertible<mpl::_,vector_type> >,boost::is_convertible<mpl::_,vector_ptrtype> > ) )
818  ( in_out( solution ),* )
819  ( rhs,( vector_ptrtype ) ) )
820  ( optional
821  //(prec,(sparse_matrix_ptrtype), matrix )
822  ( prec,( preconditioner_ptrtype ), preconditioner( _prefix=this->prefix(),_matrix=matrix,_pc=this->pcEnumType()/*LU_PRECOND*/,
823  _pcfactormatsolverpackage=this->matSolverPackageEnumType(), _backend=this->shared_from_this() ) )
824  ( maxit,( size_type ), M_maxit/*1000*/ )
825  ( rtolerance,( double ), M_rtolerance/*1e-13*/ )
826  ( atolerance,( double ), M_atolerance/*1e-50*/ )
827  ( dtolerance,( double ), M_dtolerance/*1e5*/ )
828  ( reuse_prec,( bool ), M_reuse_prec )
829  ( transpose,( bool ), false )
830  ( constant_null_space,( bool ), false )
831  ( pc,( std::string ),M_pc/*"lu"*/ )
832  ( ksp,( std::string ),M_ksp/*"gmres"*/ )
833  ( pcfactormatsolverpackage,( std::string ), M_pcFactorMatSolverPackage )
834  )
835  )
836  {
837  this->setTolerances( _dtolerance=dtolerance,
838  _rtolerance=rtolerance,
839  _atolerance=atolerance,
840  _maxit=maxit );
841 
842  this->setSolverType( _pc=pc, _ksp=ksp,
843  _constant_null_space=constant_null_space,
844  _pcfactormatsolverpackage = pcfactormatsolverpackage );
845  this->attachPreconditioner( prec );
846  // make sure matrix and rhs are closed
847  matrix->close();
848  rhs->close();
849 
850  // print them in matlab format
851  if ( !M_export.empty() )
852  {
853  matrix->printMatlab( M_export+"_A.m" );
854  rhs->printMatlab( M_export+"_b.m" );
855  }
856 
857  vector_ptrtype _sol( this->newVector( Feel::detail::datamap( solution ) ) );
858  // initialize
859  *_sol = Feel::detail::ref( solution );
860  this->setTranspose( transpose );
861  solve_return_type ret;
862 
863  if ( reuse_prec == false )
864  {
865  ret = solve( matrix, matrix, _sol, rhs );
866  }
867 
868  else
869  ret = solve( matrix, matrix, _sol, rhs, reuse_prec );
870 
871  //new
872  _sol->close();
873  detail::ref( solution ) = *_sol;
874  return ret;
875  }
876 
891  virtual solve_return_type solve( sparse_matrix_ptrtype const& A,
892  sparse_matrix_ptrtype const& P,
893  vector_ptrtype& x,
894  vector_ptrtype const& b
895  ) = 0;
896 
912  solve_return_type solve( sparse_matrix_ptrtype const& A,
913  sparse_matrix_ptrtype const& P,
914  vector_ptrtype& x,
915  vector_ptrtype const& b,
916  bool reuse_prec
917  );
918 
919 
923  BOOST_PARAMETER_MEMBER_FUNCTION( ( nl_solve_return_type ),
924  nlSolve,
925  tag,
926  ( required
927  //( in_out( solution ),*( mpl::or_<boost::is_convertible<mpl::_,vector_type&>,boost::is_convertible<mpl::_,vector_ptrtype> > ) ) )
928  ( in_out( solution ),*))
929  ( optional
930  ( jacobian,( sparse_matrix_ptrtype ), sparse_matrix_ptrtype() )
931  ( residual,( vector_ptrtype ), vector_ptrtype() )
932  //(prec,(sparse_matrix_ptrtype), jacobian )
933  ( prec,( preconditioner_ptrtype ), preconditioner( _prefix=this->prefix(),_pc=this->pcEnumType()/*LU_PRECOND*/,_backend=this->shared_from_this(),
934  _pcfactormatsolverpackage=this->matSolverPackageEnumType() ) )
935  ( maxit,( size_type ), M_maxitSNES/*50*/ )
936  ( rtolerance,( double ), M_rtoleranceSNES/*1e-8*/ )
937  ( atolerance,( double ), M_atoleranceSNES/*1e-50*/ )
938  ( stolerance,( double ), M_stoleranceSNES/*1e-8*/ )
939  ( reuse_prec,( bool ), M_reuse_prec )
940  ( reuse_jac,( bool ), M_reuse_jac )
941  ( transpose,( bool ), false )
942  ( pc,( std::string ),M_pc/*"lu"*/ )
943  ( ksp,( std::string ),M_ksp/*"gmres"*/ )
944  ( pcfactormatsolverpackage,( std::string ), M_pcFactorMatSolverPackage )
945  )
946  )
947  {
948  this->setTolerancesSNES( _stolerance=stolerance,
949  _rtolerance=rtolerance,
950  _atolerance=atolerance,
951  _maxit=maxit );
952  this->setSolverType( _pc=pc, _ksp=ksp,
953  _pcfactormatsolverpackage = pcfactormatsolverpackage );
954  vector_ptrtype _sol( this->newVector( detail::datamap( solution ) ) );
955  // initialize
956  *_sol = detail::ref( solution );
957  this->setTranspose( transpose );
958  solve_return_type ret;
959 
960  // this is done with nonlinerarsolver
961  if ( !residual )
962  {
963  residual = this->newVector( ( detail::datamap( solution ) ) );
964  //this->nlSolver()->residual( _sol, residual );
965  }
966 
967  this->nlSolver()->setPrefix( this->prefix() );
968  if ( !jacobian )
969  this->nlSolver()->jacobian( _sol, jacobian );
970 
971  if ( prec && !this->nlSolver()->initialized() )
972  this->nlSolver()->attachPreconditioner( prec );
973 
974  if ( reuse_prec == false && reuse_jac == false )
975  ret = nlSolve( jacobian, _sol, residual, rtolerance, maxit );
976 
977  else
978  ret = nlSolve( jacobian, _sol, residual, rtolerance, maxit, reuse_prec, reuse_jac );
979 
980  //new
981  _sol->close();
982  detail::ref( solution ) = *_sol;
983  detail::ref( solution ).close();
984  return ret;
985  }
986 
990  virtual nl_solve_return_type nlSolve( sparse_matrix_ptrtype& A,
991  vector_ptrtype& x,
992  vector_ptrtype& b,
993  const double, const int );
994 
999  virtual nl_solve_return_type nlSolve( sparse_matrix_ptrtype& A,
1000  vector_ptrtype& x,
1001  vector_ptrtype& b,
1002  const double, const int,
1003  bool reusePC, bool reuseJAC );
1004 
1008  void attachPreconditioner( preconditioner_ptrtype preconditioner )
1009  {
1010  if ( M_preconditioner && M_preconditioner != preconditioner )
1011  M_preconditioner->clear();
1012  M_preconditioner = preconditioner;
1013  }
1014 
1018  template<typename Observer>
1019  void
1020  addDeleteObserver( Observer const& obs )
1021  {
1022  M_deleteObservers.connect( obs );
1023  }
1028  template<typename Observer>
1029  void
1030  addDeleteObserver( boost::shared_ptr<Observer> const& obs )
1031  {
1032  M_deleteObservers.connect(boost::bind(&Observer::operator(), obs));
1033  }
1037  void
1039  {
1040  M_deleteObservers();
1041  }
1043 
1044 
1045 
1046 protected:
1047  preconditioner_ptrtype M_preconditioner;
1048 private:
1049 
1050  void start();
1051 
1052  void stop();
1053 
1054  void reset();
1055 
1056 private:
1057  WorldComm M_worldComm;
1058 
1059  po::variables_map M_vm;
1060 
1061  BackendType M_backend;
1062 
1063  std::string M_prefix;
1064 
1065  solvernonlinear_ptrtype M_nlsolver;
1066 
1067  MatrixStructure M_prec_matrix_structure;
1068 
1069  double M_totalSolveIter;
1070  double M_lastSolveIter;
1071  double M_firstSolveTime;
1072  double M_residual;
1073  double M_rtolerance;
1074  double M_dtolerance;
1075  double M_atolerance;
1076  double M_rtoleranceSNES, M_stoleranceSNES, M_atoleranceSNES;
1077  double M_rtoleranceKSPinSNES;
1078 
1079  bool M_reuse_prec;
1080  bool M_reuse_jac;
1081  bool M_reusePrecIsBuild,M_reusePrecRebuildAtFirstNewtonStep;
1082  bool M_reuseJacIsBuild,M_reuseJacRebuildAtFirstNewtonStep;
1083  size_t M_nUsePC;
1084  bool M_converged;
1085  bool M_reusePC;
1086  bool M_reusedPC;
1087  bool M_reuseFailed;
1088  boost::timer M_timer;
1089  bool M_transpose;
1090  size_type M_maxit, M_maxitSNES;
1091  size_type M_iteration;
1092  std::string M_export;
1093  std::string M_ksp;
1094  std::string M_pc;
1095  std::string M_fieldSplit;
1096  std::string M_pcFactorMatSolverPackage;
1097  bool M_constant_null_space;
1098  bool M_showKSPMonitor, M_showSNESMonitor;
1099  bool M_showKSPConvergedReason, M_showSNESConvergedReason;
1100  //std::map<std::string,boost::tuple<std::string,std::string> > M_sub;
1101 
1102  boost::signals2::signal<void()> M_deleteObservers;
1103 };
1104 
1105 
1106 typedef Backend<double> backend_type;
1107 typedef boost::shared_ptr<backend_type> backend_ptrtype;
1108 
1109 
1110 namespace detail
1111 {
1112 class BackendManagerImpl:
1113  public std::map<std::pair<BackendType,std::string>, backend_ptrtype >,
1114  public boost::noncopyable
1115 {
1116 public:
1117  typedef backend_ptrtype value_type;
1118  typedef std::pair<BackendType,std::string> key_type;
1119  typedef std::map<key_type, value_type> backend_manager_type;
1120 
1121 };
1123 
1124 struct BackendManagerDeleterImpl
1125 {
1126  void operator()() const
1127  {
1128  VLOG(2) << "[BackendManagerDeleter] clear BackendManager Singleton: " << detail::BackendManager::instance().size() << "\n";
1130  VLOG(2) << "[BackendManagerDeleter] clear BackendManager done\n";
1131  }
1132 };
1134 } // detail
1135 
1136 
1137 BOOST_PARAMETER_FUNCTION(
1138  ( backend_ptrtype ), // return type
1139  backend, // 2. function name
1140  tag, // 3. namespace of tag types
1141  ( optional
1142  ( vm, ( po::variables_map ), Environment::vm() )
1143  ( name, ( std::string ), "" )
1144  ( kind, ( BackendType ), BACKEND_PETSC )
1145  ( rebuild, ( bool ), false )
1146  ) )
1147 {
1148  // register the BackendManager into Feel::Environment so that it gets the
1149  // BackendManager is cleared up when the Environment is deleted
1150  static bool observed=false;
1151  if ( !observed )
1152  {
1153  Environment::addDeleteObserver( detail::BackendManagerDeleter::instance() );
1154  observed = true;
1155  }
1156 
1157 
1158  Feel::detail::ignore_unused_variable_warning( args );
1159 
1160  auto git = detail::BackendManager::instance().find( std::make_pair( kind, name ) );
1161 
1162  if ( git != detail::BackendManager::instance().end() && ( rebuild == false ) )
1163  {
1164  VLOG(2) << "[backend] found backend name=" << name << " kind=" << kind << " rebuild=" << rebuild << "\n";
1165  return git->second;
1166  }
1167 
1168  else
1169  {
1170  if ( git != detail::BackendManager::instance().end() && ( rebuild == true ) )
1171  git->second->clear();
1172 
1173  VLOG(2) << "[backend] building backend name=" << name << " kind=" << kind << " rebuild=" << rebuild << "\n";
1174 
1175  backend_ptrtype b;
1176  if ( vm.empty() )
1177  {
1178  b = Feel::backend_type::build( kind );
1179  }
1180  else
1181  b = Feel::backend_type::build( vm, name );
1182  VLOG(2) << "storing backend in singleton" << "\n";
1183  detail::BackendManager::instance().operator[]( std::make_pair( kind, name ) ) = b;
1184  return b;
1185  }
1186 
1187 }
1188 
1189 }
1190 #endif /* __Backend_H */
WorldComm const & comm() const
Definition: backend.hpp:640
void stencilManagerAdd(StencilManagerImpl::key_type const &key, StencilManagerImpl::graph_ptrtype graph)
function that add an entry in the StencilManager
Definition: stencil.cpp:87
sparse_matrix_ptrtype newBlockMatrixImpl(vf::BlocksBase< sparse_matrix_ptrtype > const &b, bool copy_values=true, bool diag_is_nonzero=true)
Definition: backend.hpp:355
!
Definition: backend.hpp:120
size_type maxIterations() const
Definition: backend.hpp:600
virtual nl_solve_return_type nlSolve(sparse_matrix_ptrtype &A, vector_ptrtype &x, vector_ptrtype &b, const double, const int)
Definition: backend.cpp:442
value_type dTolerance() const
Definition: backend.hpp:568
Definition: solverlinear.hpp:33
PreconditionerType
Definition: feelalg/enums.hpp:114
void attachPreconditioner(preconditioner_ptrtype preconditioner)
Definition: backend.hpp:1008
MatrixStructure precMatrixStructure() const
Definition: backend.hpp:544
std::string pcType() const
Definition: backend.hpp:492
virtual sparse_matrix_ptrtype newMatrix(const size_type m, const size_type n, const size_type m_l, const size_type n_l, const size_type nnz=30, const size_type noz=10, size_type prop=NON_HERMITIAN)=0
Non linear solver base interface.
Definition: solvernonlinear.hpp:62
PreconditionerType pcEnumType() const
Definition: backend.cpp:594
virtual void clear()
Definition: backend.cpp:157
bool hasConstantNullSpace() const
Definition: backend.hpp:500
static singleton_type & instance()
Definition: singleton.hpp:103
SolverType kspEnumType() const
Definition: backend.cpp:562
void addDeleteObserver(boost::shared_ptr< Observer > const &obs)
Definition: backend.hpp:1030
sparse_matrix_ptrtype newMatrix(const size_type m, const size_type n, const size_type m_l, const size_type n_l, graph_ptrtype const &graph, std::vector< std::vector< size_type > > indexSplit, size_type matrix_properties=NON_HERMITIAN)
Definition: backend.hpp:222
MatrixStructure
Definition: feelalg/enums.hpp:146
value_type rTolerance() const
Definition: backend.hpp:552
void prod(sparse_matrix_ptrtype const &A, vector_ptrtype const &x, vector_ptrtype &y) const
Definition: backend.hpp:790
base class for all linear algebra backends
Definition: backend.hpp:133
static backend_ptrtype build(BackendType=BACKEND_GMM, WorldComm const &worldComm=Environment::worldComm())
Definition: backend.cpp:167
value_type rtoleranceKSPinSNES() const
Definition: backend.hpp:616
SolverType
Definition: feelalg/enums.hpp:89
virtual real_type dot(vector_type const &x, vector_type const &y) const
Definition: backend.cpp:480
void sendDeleteSignal()
Definition: backend.hpp:1038
value_type rToleranceSNES() const
Definition: backend.hpp:560
sparse_matrix_ptrtype newMatrix(datamap_ptrtype const &domainmap, datamap_ptrtype const &imagemap, graph_ptrtype const &graph, size_type matrix_properties=NON_HERMITIAN, bool init=true)
Definition: backend.hpp:247
FieldSplitType
Definition: feelalg/enums.hpp:138
Definition: feelalg/enums.hpp:67
data layout in a multi-processor environnement
Definition: datamap.hpp:46
virtual vector_ptrtype newVector(datamap_ptrtype const &dm)=0
virtual solve_return_type solve(sparse_matrix_ptrtype const &A, sparse_matrix_ptrtype const &P, vector_ptrtype &x, vector_ptrtype const &b)=0
std::string prefix() const
Definition: backend.hpp:476
std::string pcFactorMatSolverPackageType() const
Definition: backend.hpp:531
real_type dot(vector_ptrtype const &x, vector_ptrtype const &y) const
Definition: backend.hpp:778
virtual void prod(sparse_matrix_type const &A, vector_type const &x, vector_type &y) const =0
std::string kspType() const
Definition: backend.hpp:484
size_t size_type
Indices (starting from 0)
Definition: feelcore/feel.hpp:319
value_type aTolerance() const
Definition: backend.hpp:584
Definition: matrixsparse.hpp:76
void setPrecMatrixStructure(MatrixStructure mstruct)
Definition: backend.hpp:728
matrices that define its action against a vector
Definition: matrixshell.hpp:49
value_type aToleranceSNES() const
Definition: backend.hpp:592
virtual sparse_matrix_ptrtype newZeroMatrix(const size_type m, const size_type n, const size_type m_l, const size_type n_l)=0
implement the Singleton pattern
Definition: singleton.hpp:57
value_type sToleranceSNES() const
Definition: backend.hpp:576
void addDeleteObserver(Observer const &obs)
Definition: backend.hpp:1020
BackendType
Definition: feelalg/enums.hpp:77
size_type maxIterationsSNES() const
Definition: backend.hpp:608
block of vector
Definition: backend.hpp:123
block of matrices
Definition: backend.hpp:121
Graph representation of the Compressed Sparse Row format.
Definition: graphcsr.hpp:57
std::string fieldsplitType() const
Definition: backend.hpp:508
vector_ptrtype newBlockVectorImpl(vf::BlocksBase< BlockType > const &b, bool copy_values=true)
Definition: backend.hpp:395
FieldSplitType fieldSplitEnumType() const
Definition: backend.cpp:633
solvernonlinear_ptrtype nlSolver()
Definition: backend.hpp:736
MatSolverPackageType matSolverPackageEnumType() const
Definition: backend.cpp:647

Generated on Sun Dec 22 2013 13:10:54 for Feel++ by doxygen 1.8.5