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>
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>
63 #include <feel/feelvf/pattern.hpp>
80 boost::shared_ptr<DataMap> datamap( T
const& t, mpl::true_ )
85 boost::shared_ptr<DataMap> datamap( T
const& t, mpl::false_ )
90 boost::shared_ptr<DataMap> datamap( T
const& t )
92 return datamap( t, detail::is_shared_ptr<T>() );
96 #if BOOST_VERSION >= 105300
97 typename boost::detail::sp_dereference< typename T::element_type >::type
101 ref( T t, mpl::true_ )
106 T& ref( T& t, mpl::false_ )
111 auto ref( T& t ) -> decltype( ref( t, detail::is_shared_ptr<T>() ) )
113 return ref( t, detail::is_shared_ptr<T>() );
122 template<
typename T>
class VectorBlockBase;
134 public boost::enable_shared_from_this<Backend<T> >
142 typedef T value_type;
143 typedef typename type_traits<T>::real_type real_type;
146 typedef boost::shared_ptr<vector_type> vector_ptrtype;
148 typedef boost::shared_ptr<sparse_matrix_type> sparse_matrix_ptrtype;
151 typedef boost::shared_ptr<shell_matrix_type> shell_matrix_ptrtype;
154 typedef typename sparse_matrix_type::graph_ptrtype graph_ptrtype;
157 typedef boost::shared_ptr<backend_type> backend_ptrtype;
160 typedef boost::shared_ptr<solvernonlinear_type> solvernonlinear_ptrtype;
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;
166 typedef boost::shared_ptr<datamap_type> datamap_ptrtype;
174 Backend( WorldComm
const& worldComm=Environment::worldComm() );
175 Backend( po::variables_map
const& vm, std::string
const&
prefix =
"", WorldComm
const& worldComm=Environment::worldComm() );
184 static backend_ptrtype
build(
185 #
if defined( FEELPP_HAS_PETSC_H )
190 , WorldComm
const& worldComm=Environment::worldComm()
196 static backend_ptrtype
build( po::variables_map
const& vm, std::string
const&
prefix =
"", WorldComm
const& worldComm=Environment::worldComm() );
216 graph_ptrtype
const & graph,
226 graph_ptrtype
const & graph,
227 std::vector < std::vector<size_type> > indexSplit,
230 auto mat = this->
newMatrix( m,n,m_l,n_l,graph,matrix_properties );
231 mat->setIndexSplit( indexSplit );
239 virtual sparse_matrix_ptrtype
newMatrix( datamap_ptrtype
const& dm1,
240 datamap_ptrtype
const& dm2,
242 bool init =
true ) = 0;
247 sparse_matrix_ptrtype
newMatrix( datamap_ptrtype
const& domainmap,
248 datamap_ptrtype
const& imagemap,
249 graph_ptrtype
const & graph,
253 auto mat = this->
newMatrix( domainmap,imagemap, matrix_properties,
false );
255 if ( init ) mat->init( imagemap->nDof(), domainmap->nDof(),
256 imagemap->nLocalDofWithoutGhost(), domainmap->nLocalDofWithoutGhost(),
270 sparse_matrix_ptrtype
276 virtual sparse_matrix_ptrtype
newZeroMatrix( datamap_ptrtype
const& dm1, datamap_ptrtype
const& dm2 ) = 0;
281 BOOST_PARAMETER_MEMBER_FUNCTION( ( sparse_matrix_ptrtype ),
285 ( trial,*( boost::is_convertible<mpl::_,boost::shared_ptr<FunctionSpaceBase> > ) )
286 ( test,*( boost::is_convertible<mpl::_,boost::shared_ptr<FunctionSpaceBase> > ) ) )
288 ( pattern,(
size_type ),Pattern::COUPLED )
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 )
299 auto mat = this->
newMatrix( trial->dofOnOff(), test->dofOn(), properties, false );
301 if ( !buildGraphWithTranspose )
303 auto s = stencil( _test=test,
306 _pattern_block=pattern_block,
307 _diag_is_nonzero=diag_is_nonzero,
308 _collect_garbage=collect_garbage);
310 mat->init( test->nDof(), trial->nDof(),
311 test->nLocalDofWithoutGhost(), trial->nLocalDofWithoutGhost(),
316 auto s = stencil( _test=trial,
319 _pattern_block=pattern_block.transpose(),
320 _diag_is_nonzero=
false,
322 _collect_garbage=collect_garbage );
324 auto graph = s->graph()->transpose(
false);
325 if ( diag_is_nonzero )
326 graph->addMissingZeroEntriesDiagonal();
332 stencilManagerAdd(boost::make_tuple( test, trial, pattern, pattern_block.getSetOfBlocks(), diag_is_nonzero), graph);
334 mat->init( test->nDof(), trial->nDof(),
335 test->nLocalDofWithoutGhost(), trial->nLocalDofWithoutGhost(),
340 mat->setIndexSplit( trial->dofIndexSplit() );
344 template<
typename DomainSpace,
typename ImageSpace>
347 sparse_matrix_ptrtype m =
newMatrix( dm, im, prop );
348 m->init( im->nDof(), dm->nDof(), im->nLocalDof(), dm->nLocalDof(), M->graph() );
356 bool copy_values=
true,
357 bool diag_is_nonzero=
true )
360 boost::shared_ptr<matrix_block_type> mb(
new matrix_block_type( b, *
this, copy_values, diag_is_nonzero ) );
361 return mb->getSparseMatrix();
364 sparse_matrix_ptrtype
newBlockMatrixImpl( vf::BlocksBase<boost::shared_ptr<GraphCSR> >
const & b,
365 bool copy_values=
true,
366 bool diag_is_nonzero=
true )
369 boost::shared_ptr<matrix_block_type> mb(
new matrix_block_type( b, *
this, diag_is_nonzero ) );
370 return mb->getSparseMatrix();
376 BOOST_PARAMETER_MEMBER_FUNCTION( ( sparse_matrix_ptrtype ),
383 ( copy_values,*( boost::is_integral<mpl::_> ),
true )
384 ( diag_is_nonzero, *( boost::is_integral<mpl::_> ),
true )
394 template <
typename BlockType=vector_ptrtype >
396 bool copy_values=
true )
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();
406 BOOST_PARAMETER_MEMBER_FUNCTION( ( vector_ptrtype ),
413 ( copy_values,*( boost::is_integral<mpl::_> ),
true )
424 BOOST_PARAMETER_MEMBER_FUNCTION( ( sparse_matrix_ptrtype ),
428 ( test,*( boost::is_convertible<mpl::_,boost::shared_ptr<FunctionSpaceBase> >) )
429 ( trial,*( boost::is_convertible<mpl::_,boost::shared_ptr<FunctionSpaceBase> >) )
433 return this->
newZeroMatrix( trial->dofOnOff(), test->dofOn() );
439 virtual vector_ptrtype
newVector( datamap_ptrtype
const& dm ) = 0;
449 BOOST_PARAMETER_MEMBER_FUNCTION( ( vector_ptrtype ),
453 ( test,*( boost::is_convertible<mpl::_,boost::shared_ptr<FunctionSpaceBase> >) )
502 return M_constant_null_space;
533 return M_pcFactorMatSolverPackage;
546 return M_prec_matrix_structure;
562 return M_rtoleranceSNES;
578 return M_stoleranceSNES;
594 return M_atoleranceSNES;
618 return M_rtoleranceKSPinSNES;
622 bool converged()
const
632 bool transpose()
const
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; }
650 bool reusePrec()
const {
return M_reuse_prec; }
651 bool reuseJac()
const {
return M_reuse_jac; }
653 bool reusePrecRebuildAtFirstNewtonStep()
const {
return M_reusePrecRebuildAtFirstNewtonStep; }
654 bool reuseJacRebuildAtFirstNewtonStep()
const {
return M_reuseJacRebuildAtFirstNewtonStep; }
668 BOOST_PARAMETER_MEMBER_FUNCTION( (
void ),
672 ( rtolerance, (
double ) )
676 ( atolerance, (
double ), 1e-50 )
677 ( dtolerance, (
double ), 1e5 )
680 M_rtolerance = rtolerance;
681 M_dtolerance = dtolerance;
682 M_atolerance = atolerance;
686 BOOST_PARAMETER_MEMBER_FUNCTION( (
void ),
690 ( rtolerance, (
double ) )
694 ( atolerance, (
double ), 1e-50 )
695 ( stolerance, (
double ), 1e-8 )
698 M_rtoleranceSNES = rtolerance;
699 M_stoleranceSNES = stolerance;
700 M_atoleranceSNES = atolerance;
707 BOOST_PARAMETER_MEMBER_FUNCTION( (
void ),
711 ( ksp, ( std::string ) )
714 ( pc, ( std::string ),
"lu" )
715 ( constant_null_space, (
bool ),
false )
716 ( pcfactormatsolverpackage, ( std::string ),
"petsc" )
721 M_pcFactorMatSolverPackage = pcfactormatsolverpackage;
722 M_constant_null_space = constant_null_space;
730 M_prec_matrix_structure = mstruct;
741 void setTranspose(
bool transpose )
743 M_transpose = transpose;
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; }
751 void setReusePrec(
bool b ) { M_reuse_prec=b; }
752 void setReuseJac(
bool b) { M_reuse_jac=b; }
754 void setReusePrecRebuildAtFirstNewtonStep(
bool b) { M_reusePrecRebuildAtFirstNewtonStep=b; }
755 void setReuseJacRebuildAtFirstNewtonStep(
bool b) { M_reuseJacRebuildAtFirstNewtonStep=b; }
767 virtual void clear();
772 virtual real_type
dot( vector_type
const& x, vector_type
const& y )
const;
778 real_type
dot( vector_ptrtype
const& x, vector_ptrtype
const& y )
const
780 return this->
dot( *x, *y );
790 void prod( sparse_matrix_ptrtype
const& A, vector_ptrtype
const& x, vector_ptrtype& y )
const
792 this->
prod( *A, *x, *y );
811 BOOST_PARAMETER_MEMBER_FUNCTION( ( solve_return_type ),
816 ( matrix,(sparse_matrix_ptrtype) )
818 ( in_out( solution ),* )
819 ( rhs,( vector_ptrtype ) ) )
822 ( prec,( preconditioner_ptrtype ), preconditioner( _prefix=this->
prefix(),_matrix=matrix,_pc=this->
pcEnumType(),
825 ( rtolerance,(
double ), M_rtolerance )
826 ( atolerance,(
double ), M_atolerance )
827 ( dtolerance,(
double ), M_dtolerance )
828 ( reuse_prec,(
bool ), M_reuse_prec )
829 ( transpose,(
bool ),
false )
830 ( constant_null_space,(
bool ),
false )
831 ( pc,( std::string ),M_pc )
832 ( ksp,( std::string ),M_ksp )
833 ( pcfactormatsolverpackage,( std::string ), M_pcFactorMatSolverPackage )
837 this->setTolerances( _dtolerance=dtolerance,
838 _rtolerance=rtolerance,
839 _atolerance=atolerance,
842 this->setSolverType( _pc=pc, _ksp=ksp,
843 _constant_null_space=constant_null_space,
844 _pcfactormatsolverpackage = pcfactormatsolverpackage );
851 if ( !M_export.empty() )
853 matrix->printMatlab( M_export+
"_A.m" );
854 rhs->printMatlab( M_export+
"_b.m" );
857 vector_ptrtype _sol( this->
newVector( Feel::detail::datamap( solution ) ) );
859 *_sol = Feel::detail::ref( solution );
860 this->setTranspose( transpose );
861 solve_return_type ret;
863 if ( reuse_prec ==
false )
865 ret =
solve( matrix, matrix, _sol, rhs );
869 ret =
solve( matrix, matrix, _sol, rhs, reuse_prec );
873 detail::ref( solution ) = *_sol;
891 virtual solve_return_type
solve( sparse_matrix_ptrtype
const& A,
892 sparse_matrix_ptrtype
const& P,
894 vector_ptrtype
const& b
912 solve_return_type
solve( sparse_matrix_ptrtype
const& A,
913 sparse_matrix_ptrtype
const& P,
915 vector_ptrtype
const& b,
923 BOOST_PARAMETER_MEMBER_FUNCTION( ( nl_solve_return_type ),
928 ( in_out( solution ),*))
930 ( jacobian,( sparse_matrix_ptrtype ), sparse_matrix_ptrtype() )
931 ( residual,( vector_ptrtype ), vector_ptrtype() )
933 ( prec,( preconditioner_ptrtype ), preconditioner( _prefix=this->
prefix(),_pc=this->
pcEnumType(),_backend=this->shared_from_this(),
936 ( rtolerance,(
double ), M_rtoleranceSNES )
937 ( atolerance,(
double ), M_atoleranceSNES )
938 ( stolerance,(
double ), M_stoleranceSNES )
939 ( reuse_prec,(
bool ), M_reuse_prec )
940 ( reuse_jac,(
bool ), M_reuse_jac )
941 ( transpose,(
bool ),
false )
942 ( pc,( std::string ),M_pc )
943 ( ksp,( std::string ),M_ksp )
944 ( pcfactormatsolverpackage,( std::string ), M_pcFactorMatSolverPackage )
948 this->setTolerancesSNES( _stolerance=stolerance,
949 _rtolerance=rtolerance,
950 _atolerance=atolerance,
952 this->setSolverType( _pc=pc, _ksp=ksp,
953 _pcfactormatsolverpackage = pcfactormatsolverpackage );
954 vector_ptrtype _sol( this->
newVector( detail::datamap( solution ) ) );
956 *_sol = detail::ref( solution );
957 this->setTranspose( transpose );
958 solve_return_type ret;
963 residual = this->
newVector( ( detail::datamap( solution ) ) );
969 this->
nlSolver()->jacobian( _sol, jacobian );
971 if ( prec && !this->
nlSolver()->initialized() )
972 this->
nlSolver()->attachPreconditioner( prec );
974 if ( reuse_prec ==
false && reuse_jac ==
false )
975 ret =
nlSolve( jacobian, _sol, residual, rtolerance, maxit );
978 ret =
nlSolve( jacobian, _sol, residual, rtolerance, maxit, reuse_prec, reuse_jac );
982 detail::ref( solution ) = *_sol;
983 detail::ref( solution ).close();
990 virtual nl_solve_return_type
nlSolve( sparse_matrix_ptrtype& A,
993 const double,
const int );
999 virtual nl_solve_return_type
nlSolve( sparse_matrix_ptrtype& A,
1002 const double,
const int,
1003 bool reusePC,
bool reuseJAC );
1010 if ( M_preconditioner && M_preconditioner != preconditioner )
1011 M_preconditioner->clear();
1012 M_preconditioner = preconditioner;
1018 template<
typename Observer>
1022 M_deleteObservers.connect( obs );
1028 template<
typename Observer>
1032 M_deleteObservers.connect(boost::bind(&Observer::operator(), obs));
1040 M_deleteObservers();
1047 preconditioner_ptrtype M_preconditioner;
1057 WorldComm M_worldComm;
1059 po::variables_map M_vm;
1063 std::string M_prefix;
1065 solvernonlinear_ptrtype M_nlsolver;
1069 double M_totalSolveIter;
1070 double M_lastSolveIter;
1071 double M_firstSolveTime;
1073 double M_rtolerance;
1074 double M_dtolerance;
1075 double M_atolerance;
1076 double M_rtoleranceSNES, M_stoleranceSNES, M_atoleranceSNES;
1077 double M_rtoleranceKSPinSNES;
1081 bool M_reusePrecIsBuild,M_reusePrecRebuildAtFirstNewtonStep;
1082 bool M_reuseJacIsBuild,M_reuseJacRebuildAtFirstNewtonStep;
1088 boost::timer M_timer;
1092 std::string M_export;
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;
1102 boost::signals2::signal<void()> M_deleteObservers;
1106 typedef Backend<double> backend_type;
1107 typedef boost::shared_ptr<backend_type> backend_ptrtype;
1112 class BackendManagerImpl:
1113 public std::map<std::pair<BackendType,std::string>, backend_ptrtype >,
1114 public boost::noncopyable
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;
1124 struct BackendManagerDeleterImpl
1126 void operator()()
const
1130 VLOG(2) <<
"[BackendManagerDeleter] clear BackendManager done\n";
1137 BOOST_PARAMETER_FUNCTION(
1138 ( backend_ptrtype ),
1142 ( vm, ( po::variables_map ), Environment::vm() )
1143 ( name, ( std::string ),
"" )
1145 ( rebuild, (
bool ),
false )
1150 static bool observed=
false;
1158 Feel::detail::ignore_unused_variable_warning( args );
1164 VLOG(2) <<
"[backend] found backend name=" << name <<
" kind=" << kind <<
" rebuild=" << rebuild <<
"\n";
1171 git->second->clear();
1173 VLOG(2) <<
"[backend] building backend name=" << name <<
" kind=" << kind <<
" rebuild=" << rebuild <<
"\n";
1182 VLOG(2) <<
"storing backend in singleton" <<
"\n";