34 #ifndef _BACKENDTRILINOS_HPP_
35 #define _BACKENDTRILINOS_HPP_
38 #include <boost/program_options/variables_map.hpp>
39 #include <feel/feelconfig.h>
41 #if defined ( FEELPP_HAS_TRILINOS_EPETRA )
42 #undef PACKAGE_BUGREPORT
45 #undef PACKAGE_TARNAME
46 #undef PACKAGE_VERSION
49 #include <feel/feelalg/preconditionerifpack.hpp>
50 #include <feel/feelalg/preconditionerml.hpp>
53 #include <feel/feelalg/datamap.hpp>
57 #include <Teuchos_ParameterList.hpp>
61 namespace po = boost::program_options;
64 template<
class T>
class Backend;
72 class BackendTrilinos :
public Backend<double>
74 typedef Backend<double> super;
78 typedef super::value_type value_type;
80 typedef super::vector_ptrtype base_vector_ptrtype;
81 typedef super::sparse_matrix_ptrtype base_sparse_matrix_ptrtype;
82 typedef super::vector_type vector_type;
83 typedef super::sparse_matrix_type sparse_matrix_type;
85 typedef sparse_matrix_type::graph_type graph_type;
86 typedef sparse_matrix_type::graph_ptrtype graph_ptrtype;
88 typedef MatrixEpetra epetra_sparse_matrix_type;
89 typedef VectorEpetra<value_type> epetra_vector_type;
90 typedef Epetra_Operator operator_type;
92 typedef OperatorMatrix op_mat_type;
93 typedef boost::shared_ptr<OperatorMatrix> op_mat_ptrtype;
95 typedef boost::shared_ptr<epetra_sparse_matrix_type> epetra_sparse_matrix_ptrtype;
96 typedef boost::shared_ptr<epetra_vector_type> epetra_vector_ptrtype;
97 typedef boost::shared_ptr<operator_type> operator_ptrtype;
99 typedef super::solve_return_type solve_return_type;
100 typedef super::nl_solve_return_type nl_solve_return_type;
102 typedef std::map< size_type, std::vector< size_type > > procdist_type;
104 typedef Teuchos::ParameterList list_type;
108 BackendTrilinos( WorldComm
const& worldComm=Environment::worldComm())
120 BackendTrilinos( po::variables_map
const& vm, std::string
const& prefix =
"", WorldComm
const& worldComm=Environment::worldComm() );
122 BackendTrilinos(
const BackendTrilinos& tc );
126 static Epetra_Map epetraMap( DataMap
const& dmap )
128 std::vector<int> e( dmap.nMyElements() );
129 std::copy( dmap.myGlobalElements().begin(),
130 dmap.myGlobalElements().end(),
132 return Epetra_Map( -1, dmap.nMyElements(), e.data(), 0, Epetra_MpiComm( dmap.comm() ) );
136 static Epetra_Map epetraMapStatic( DataMap
const& dmap )
138 return Epetra_Map( dmap.nGlobalElements(), dmap.nMyElements(), 0, Epetra_MpiComm( dmap.comm() ) );
142 template<
typename DomainSpace,
typename DualImageSpace>
143 sparse_matrix_ptrtype newMatrix( DomainSpace
const& Xh,
144 DualImageSpace
const& Yh,
148 return newMatrix( Xh->map(), Yh->map(), matrix_properties );
151 sparse_matrix_ptrtype
160 sparse_matrix_ptrtype mat(
new epetra_sparse_matrix_type( m,n,m_l,n_l,nnz,noz ) );
161 mat->setMatrixProperties( matrix_properties );
165 sparse_matrix_ptrtype
170 graph_ptrtype
const & graph,
173 sparse_matrix_ptrtype mat(
new epetra_sparse_matrix_type( m,n,m_l,n_l,30,10 ) );
174 mat->setMatrixProperties( matrix_properties );
178 sparse_matrix_ptrtype newMatrix( DataMap
const& domainmap,
179 DataMap
const& imagemap,
183 Epetra_Map erowmap = BackendTrilinos::epetraMap( imagemap );
184 Epetra_Map ecolmap = BackendTrilinos::epetraMapStatic( domainmap );
185 Epetra_Map edomainmap = BackendTrilinos::epetraMap( imagemap );
193 if ( !erowmap.SameAs( ecolmap ) )
195 Epetra_Map eimagemap = BackendTrilinos::epetraMap( domainmap );
197 auto A= sparse_matrix_ptrtype(
new epetra_sparse_matrix_type( erowmap, ecolmap, edomainmap, eimagemap ) );
198 A->setMatrixProperties( matrix_properties );
204 auto A= sparse_matrix_ptrtype(
new epetra_sparse_matrix_type( erowmap, ecolmap ) );
205 A->setMatrixProperties( matrix_properties );
210 sparse_matrix_ptrtype
216 sparse_matrix_ptrtype mat(
new epetra_sparse_matrix_type( m,n,m_l,n_l,0,0 ) );
220 sparse_matrix_ptrtype
221 newZeroMatrix( DataMap
const& domainmap, DataMap
const& imagemap )
223 Epetra_Map erowmap = BackendTrilinos::epetraMap( imagemap );
224 Epetra_Map ecolmap = BackendTrilinos::epetraMapStatic( domainmap );
226 auto A= sparse_matrix_ptrtype(
new epetra_sparse_matrix_type( erowmap, ecolmap ) );
232 template<
typename SpaceT>
233 vector_ptrtype newVector( SpaceT
const& space )
235 return newVector( space->map() );
237 vector_ptrtype newVector( DataMap
const& domainmap )
239 Epetra_Map emap = BackendTrilinos::epetraMap( domainmap );
241 return vector_ptrtype(
new epetra_vector_type( emap ) );
246 return vector_ptrtype(
new epetra_vector_type( ) );
251 static operator_ptrtype IfpackPrec( sparse_matrix_ptrtype
const& M, list_type options, std::string precType =
"Amesos" )
253 PreconditionerIfpack P( options, precType );
255 P.buildPreconditioner( M );
261 static operator_ptrtype MLPrec( sparse_matrix_ptrtype& M, list_type options )
263 PreconditionerML P( options );
265 P.buildPreconditioner( M );
271 template<
typename element_type >
272 static void Epetra2Ublas( vector_ptrtype
const& u, element_type& x )
274 epetra_vector_ptrtype
const& _v( dynamic_cast<epetra_vector_ptrtype const&>( u ) );
275 Epetra_Map v_map( _v->Map() );
285 DVLOG(2) <<
"x(" << x.firstLocalIndex() + i <<
")="
286 <<
"v[" << v_map.GID( i ) <<
"] = "
289 x( x.firstLocalIndex() + i ) = v( i );
292 DVLOG(2) <<
"Epetra2Ublas:" << x <<
"\n";
296 template<
typename element_type >
297 static void Ublas2Epetra( element_type
const& x, vector_ptrtype& v )
299 epetra_vector_type& _v( dynamic_cast<epetra_vector_type&>( *v ) );
300 Epetra_Map v_map( _v.Map() );
302 DVLOG(2) <<
"Local size of ublas vector" << x.localSize() <<
"\n";
303 DVLOG(2) <<
"Local size of epetra vector" << v->localSize() <<
"\n";
309 DVLOG(2) <<
"v[" << v_map.GID( i ) <<
"] = "
310 <<
"x(" << x.firstLocalIndex() + i <<
")="
311 << x( x.firstLocalIndex() + i ) <<
"\n";
313 v->set( v_map.GID( i ), x( x.firstLocalIndex() + i ) );
318 template<
int index,
typename spaceT >
319 static Epetra_MultiVector getComponent( spaceT
const& Xh, Epetra_MultiVector
const& sol )
321 Epetra_Map componentMap ( epetraMap( Xh->template functionSpace<index>()->map() ) );
322 Epetra_Map globalMap ( epetraMap( Xh->map() ) );
326 Epetra_MultiVector component( componentMap, 1 );
328 int Length = component.MyLength();
330 int shift = Xh->nDofStart( index );
332 for (
int i=0; i < Length; i++ )
334 int compGlobalID = componentMap.GID( i );
336 if ( compGlobalID >= 0 )
338 int compLocalID = componentMap.LID( compGlobalID );
340 int localID = globalMap.LID( compGlobalID+shift );
343 DVLOG(2) <<
"[MyBackend] Copy entry sol[" << localID <<
"]=" << sol[0][localID]
344 <<
" to component[" << compLocalID <<
"]\n";
346 component[0][compLocalID] = sol[0][localID];
348 DVLOG(2) << component[0][compLocalID] <<
"\n";
356 template<
int index,
typename spaceT >
357 static void UpdateComponent( spaceT
const& Xh, Epetra_MultiVector& sol, Epetra_MultiVector& comp )
359 Epetra_Map componentMap ( epetraMap( Xh->template functionSpace<index>()->map() ) );
360 Epetra_Map globalMap ( epetraMap( Xh->map() ) );
362 int shift = Xh->nDofStart( index );
364 int Length = comp.MyLength();
366 for (
int i=0; i < Length; i++ )
368 int compGlobalID = componentMap.GID( i );
370 if ( compGlobalID >= 0 )
372 int compLocalID = componentMap.LID( compGlobalID );
374 int localID = globalMap.LID( compGlobalID+shift );
377 DVLOG(2) <<
"Copy entry component[" << compLocalID <<
"] to sol[" << localID <<
"]="
381 sol[0][localID] = comp[0][compLocalID] ;
383 DVLOG(2) << comp[0][compLocalID] <<
"\n";
391 void set_maxiter (
int maxiter );
393 void set_tol (
double tol );
395 void set_drop (
double drop );
397 void set_overlap (
int overlap );
399 void set_restart (
int restart );
401 void set_fillin (
int fillin );
403 void set_order (
int order );
405 void set_overlap_type( std::string str );
407 void set_residual ( std::string str );
409 void set_localparts (
int localparts );
411 void set_solver( std::string str );
413 void set_prec_local_solver( std::string str );
415 void set_prec_local_solver_type( std::string str );
417 void set_prec( std::string str );
419 void set_verbose(
const std::string verb=
"none" );
421 void set_options( list_type opts );
423 list_type get_options();
425 void set_symmetric(
bool ) {}
430 applyMatrix( sparse_matrix_type
const& A,
431 const Vector<value_type>& x,
434 real_type dot(
const vector_type& f,
const vector_type& x )
const;
436 void prod( sparse_matrix_type
const& A,
437 vector_type
const& x,
438 vector_type& b )
const
440 epetra_sparse_matrix_type
const& _A =
dynamic_cast<epetra_sparse_matrix_type const&
>( A );
441 epetra_vector_type
const& _x =
dynamic_cast<epetra_vector_type const&
>( x );
442 epetra_vector_type& _b =
dynamic_cast<epetra_vector_type&
>( b );
443 _A.mat().Apply( _x.vec(), _b.vec() );
446 solve_return_type solve( base_sparse_matrix_ptrtype
const& A,
447 base_sparse_matrix_ptrtype
const& B,
448 base_vector_ptrtype& x,
449 base_vector_ptrtype
const& b );
458 mpi::communicator M_comm;
462 std::string M_prec_type, M_local_solver, M_local_solver_type;
464 operator_ptrtype M_Prec;
473 po::options_description backendtrilinos_options( std::string
const& prefix );
478 #endif // FEELPP_HAS_TRILINOS_EPETRA