30 #ifndef __operator_trilinos_matrix_H
31 #define __operator_trilinos_matrix_H 1
40 #if defined( FEELPP_HAS_TRILINOS_EPETRA )
41 #include "Epetra_ConfigDefs.h"
42 #include "Epetra_Operator.h"
45 #if defined( FEELPP_HAS_MPI )
47 #include "Epetra_MpiComm.h"
49 #include "Epetra_SerialComm.h"
51 #include "Epetra_Map.h"
52 #include "Epetra_Vector.h"
53 #include "Epetra_MultiVector.h"
54 #include "Epetra_Operator.h"
61 class OperatorMatrix :
public Epetra_Operator
65 typedef MatrixSparse<double> sparse_matrix_type;
66 typedef boost::shared_ptr<sparse_matrix_type> sparse_matrix_ptrtype;
68 typedef MatrixEpetra epetra_sparse_matrix_type;
69 typedef boost::shared_ptr<epetra_sparse_matrix_type> epetra_sparse_matrix_ptrtype;
71 typedef Vector<double> vector_type;
72 typedef boost::shared_ptr<vector_type> vector_ptrtype;
74 typedef VectorEpetra<double> epetra_vector_type;
75 typedef boost::shared_ptr<epetra_vector_type> epetra_vector_ptrtype;
77 typedef Epetra_Operator prec_type;
78 typedef boost::shared_ptr<prec_type> prec_ptrtype;
80 typedef SolverLinearTrilinos<double> solver_type;
81 typedef boost::shared_ptr<solver_type> solver_ptrtype;
83 typedef Teuchos::ParameterList list_type;
85 typedef solver_type::real_type real_type;
88 OperatorMatrix( sparse_matrix_ptrtype
const& F, std::string _label,
bool transpose = 0 )
91 M_Matrix( dynamic_cast<epetra_sparse_matrix_type const*>( F.get() ) ),
93 M_Solver( new solver_type() ),
96 M_useTranspose( transpose ),
97 M_domainMap( M_Matrix->getDomainMap() ),
98 M_rangeMap( M_Matrix->getRangeMap() ),
103 DVLOG(2) <<
"Create operator " << Label() <<
" ...\n";
106 OperatorMatrix( sparse_matrix_ptrtype
const& F,
107 list_type
const& options,
109 prec_ptrtype Prec,
bool transpose = 0 )
112 M_Matrix( dynamic_cast<epetra_sparse_matrix_type const*>( F.get() ) ),
114 M_Solver( new solver_type() ),
117 M_useTranspose( transpose ),
118 M_domainMap( M_Matrix->getDomainMap() ),
119 M_rangeMap( M_Matrix->getRangeMap() ),
124 DVLOG(2) <<
"Create operator " << Label() <<
" ...\n";
126 M_Solver->setOptions( options );
128 M_tol = M_Solver->getOptions().get(
"tol", 1e-10 );
129 M_maxiter = M_Solver->getOptions().get(
"max_iter", 100 );
134 OperatorMatrix(
const OperatorMatrix& tc )
137 M_Matrix( tc.M_Matrix ),
139 M_Solver( tc.M_Solver ),
140 M_hasInverse( tc.M_hasInverse ),
141 M_hasApply( tc.M_hasApply ),
142 M_useTranspose( tc.M_useTranspose ),
143 M_domainMap( tc.M_domainMap ),
144 M_rangeMap( tc.M_rangeMap ),
145 M_label( tc.M_label ),
146 M_maxiter( tc.M_maxiter ),
149 DVLOG(2) <<
"Copy operator " << Label() <<
" ...\n";
152 bool hasInverse()
const
157 bool hasApply()
const
163 int Apply(
const vector_ptrtype& X, vector_ptrtype& Y )
const
165 Apply( dynamic_cast<epetra_vector_type const&>( *X ).vec(),
166 dynamic_cast<epetra_vector_type&
>( *Y ).vec() );
171 int Apply(
const Epetra_MultiVector & X, Epetra_MultiVector & Y )
const
173 DVLOG(2) <<
"Apply Operator " << Label() <<
"\n";
175 M_Matrix->multiply(
false,X,Y );
177 DVLOG(2) <<
"Apply Operator " << Label() <<
" successfully\n";
181 int ApplyInverse (
const vector_ptrtype& X, vector_ptrtype& Y )
const
183 FEELPP_ASSERT( hasInverse() ).error(
"This operator cannot be inverted." );
185 int result = ApplyInverse( dynamic_cast<epetra_vector_type const&>( *X ).vec(),
186 dynamic_cast<epetra_vector_type&
>( *Y ).vec() );
191 int ApplyInverse (
const Epetra_MultiVector& X, Epetra_MultiVector& Y )
const
193 FEELPP_ASSERT( hasInverse() ).error(
"This operator cannot be inverted." );
195 std::pair<unsigned int, real_type> result;
197 DVLOG(2) <<
"Apply Inverse Operator " << Label() <<
"\n";
199 if ( M_Prec.get() != 0 )
200 result = M_Solver->solve( *M_Matrix, M_Prec, Y, X, M_tol, M_maxiter );
203 result = M_Solver->solve( *M_Matrix, Y, X, M_tol, M_maxiter );
205 DVLOG(2) <<
"Apply Inverse Operator " << Label() <<
" successfully\n";
207 int result_integer = result.first;
209 return result_integer;
213 int SetUseTranspose(
bool UseTranspose = 0 )
215 M_useTranspose = UseTranspose;
220 double NormInf()
const
222 return M_Matrix->linftyNorm();
225 const char * Label ()
const
227 return( M_label.c_str() );
230 bool UseTranspose()
const
232 return M_useTranspose;
235 bool HasNormInf ()
const
240 const Epetra_Comm& Comm()
const
242 return( M_Matrix->getDomainMap().Comm() );
245 const Epetra_Map& OperatorDomainMap()
const
247 return( M_domainMap );
250 const Epetra_Map& OperatorRangeMap()
const
252 return( M_rangeMap );
257 DVLOG(2) <<
"Destroyed matrix operator: " << this->Label() <<
" ...\n";
262 sparse_matrix_ptrtype M_F;
263 epetra_sparse_matrix_type
const* M_Matrix;
267 solver_ptrtype M_Solver;
269 bool M_hasInverse, M_hasApply, M_useTranspose;
271 Epetra_Map M_domainMap, M_rangeMap;
288 template<
typename operator_type >
289 class OperatorInverse :
public Epetra_Operator
293 typedef boost::shared_ptr<operator_type> operator_ptrtype;
296 OperatorInverse( operator_ptrtype& F )
301 DVLOG(2) <<
"Create inverse operator " << this->Label() <<
"...\n";
304 OperatorInverse(
const OperatorInverse& tc )
307 M_Label( tc.M_Label )
309 DVLOG(2) <<
"Copy inverse operator " << this->Label() <<
"...\n";
312 bool hasInverse()
const
314 return M_F->hasApply();
317 bool hasApply()
const
319 return M_F->hasInverse();
322 int Apply(
const Epetra_MultiVector & X, Epetra_MultiVector & Y )
const
324 DVLOG(2) <<
"Apply matrix " << Label() <<
"\n";
326 FEELPP_ASSERT( hasApply() ).error(
"This operator cannot be applied." );
327 M_F->ApplyInverse( X,Y );
332 int ApplyInverse (
const Epetra_MultiVector& X, Epetra_MultiVector& Y )
const
334 DVLOG(2) <<
"ApplyInverse matrix " << Label() <<
"\n";
336 FEELPP_ASSERT( hasInverse() ).error(
"This operator cannot be inverted." );
340 return !hasInverse();
344 int SetUseTranspose(
bool UseTranspose )
349 double NormInf()
const
354 const char * Label ()
const
356 return M_Label.c_str();
359 bool UseTranspose()
const
364 bool HasNormInf ()
const
370 const Epetra_Comm & Comm()
const
372 return( M_F->Comm() );
375 const Epetra_Map & OperatorDomainMap()
const
377 return( M_F->OperatorRangeMap() );
380 const Epetra_Map & OperatorRangeMap()
const
382 return( M_F->OperatorDomainMap() );
387 DVLOG(2) <<
"Destroyed inverse operator: " << this->Label() <<
" ...\n";
392 operator_ptrtype M_F;
398 std::string L = M_F->Label();
400 std::string temp(
"inv(" );
412 template<
typename op1_type,
typename op2_type >
413 class OperatorCompose :
public Epetra_Operator
417 typedef boost::shared_ptr<op1_type> op1_ptrtype;
418 typedef boost::shared_ptr<op2_type> op2_ptrtype;
427 std::string t( M_F->Label() );
428 std::string u( M_G->Label() );
435 OperatorCompose( op1_ptrtype& F, op2_ptrtype& G )
440 std::string t( F->Label() );
441 std::string u( G->Label() );
446 DVLOG(2) <<
"Create operator " << Label() <<
" ...\n";
449 OperatorCompose(
const OperatorCompose& tc )
453 M_Label( tc.M_Label )
455 DVLOG(2) <<
"Copy operator " << Label() <<
" ...\n";
458 bool hasInverse()
const
460 return M_F->hasInverse() * M_G->hasInverse();
463 bool hasApply()
const
465 return M_F->hasApply() * M_G->hasApply();
468 int Apply(
const Epetra_MultiVector & X, Epetra_MultiVector & Y )
const
470 FEELPP_ASSERT( hasApply() ).error(
"This operator cannot be applied." );
472 DVLOG(2) <<
"Apply operator " << Label() <<
" ...\n";
474 Epetra_MultiVector Z( M_G->OperatorRangeMap(), 1 );
482 int ApplyInverse (
const Epetra_MultiVector& X, Epetra_MultiVector& Y )
const
484 FEELPP_ASSERT( hasInverse() ).error(
"This operator cannot be inverted." );
486 DVLOG(2) <<
"Apply Inverse operator " << Label() <<
" ...\n";
488 Epetra_MultiVector Z( X );
490 M_F->ApplyInverse( X,Z );
491 M_G->ApplyInverse( Z,Y );
497 int SetUseTranspose(
bool UseTranspose )
502 double NormInf()
const
507 const char * Label ()
const
510 return M_Label.c_str();
513 bool UseTranspose()
const
518 bool HasNormInf ()
const
524 const Epetra_Comm & Comm()
const
526 return( M_F->Comm() );
529 const Epetra_Map & OperatorDomainMap()
const
531 return( M_G->OperatorDomainMap() );
534 const Epetra_Map & OperatorRangeMap()
const
536 return( M_F->OperatorRangeMap() );
541 DVLOG(2) <<
"Destroyed compose operator: " << this->Label() <<
" ...\n";
564 template<
typename op1_type>
565 class OperatorScale :
public Epetra_Operator
569 typedef boost::shared_ptr<op1_type> op1_ptrtype;
579 OperatorScale( op1_ptrtype& F )
584 std::string temp(
"alpha" );
585 std::string t = M_F->Label();
592 DVLOG(2) <<
"Create scale operator " << Label() <<
" ...\n";
595 OperatorScale( op1_ptrtype& F,
double alpha )
600 std::string temp(
"alpha" );
601 std::string t = M_F->Label();
608 DVLOG(2) <<
"Create scale operator " << Label() <<
" ...\n";
611 OperatorScale(
const OperatorScale& tc )
614 M_alpha( tc.M_alpha ),
615 M_Label( tc.M_Label )
617 DVLOG(2) <<
"Copy scale operator " << Label() <<
" ...\n";
620 bool hasInverse()
const
622 return M_F->hasApply()*( M_alpha != 0 );
625 bool hasApply()
const
627 return M_F->hasApply();
632 int Apply(
const Epetra_MultiVector & X, Epetra_MultiVector & Y )
const
634 DVLOG(2) <<
"Apply scale operator " << Label() <<
"\n";
643 int ApplyInverse (
const Epetra_MultiVector& X, Epetra_MultiVector& Y )
const
645 DVLOG(2) <<
"ApplyInverse scale operator " << Label() <<
"\n";
647 FEELPP_ASSERT( hasInverse() && ( M_alpha != 0 ) ).error(
"This operator cannot be inverted." );
649 Epetra_MultiVector Z( X );
650 Z.Scale( 1./M_alpha );
652 M_F->ApplyInverse( Z,Y );
654 return !hasInverse();
658 int SetUseTranspose(
bool UseTranspose )
663 double NormInf()
const
668 const char * Label ()
const
670 return M_Label.c_str();
673 bool UseTranspose()
const
678 bool HasNormInf ()
const
684 const Epetra_Comm & Comm()
const
686 return( M_F->Comm() );
689 const Epetra_Map & OperatorDomainMap()
const
691 return( M_F->OperatorDomainMap() );
694 const Epetra_Map & OperatorRangeMap()
const
696 return( M_F->OperatorRangeMap() );
718 template<
typename operator_type >
719 class OperatorFree :
public Epetra_Operator
723 typedef boost::shared_ptr<operator_type> operator_ptrtype;
725 typedef Epetra_Operator prec_type;
726 typedef boost::shared_ptr<prec_type> prec_ptrtype;
728 typedef SolverLinearTrilinos<double> solver_type;
729 typedef boost::shared_ptr<solver_type> solver_ptrtype;
731 typedef solver_type::real_type real_type;
733 typedef Teuchos::ParameterList list_type;
735 OperatorFree( operator_ptrtype F )
738 M_Solver( new solver_type() ),
740 M_hasApply( F->hasApply() ),
741 M_useTranspose( F->UseTranspose() ),
742 M_label( F->Label() )
744 DVLOG(2) <<
"Create operator " << Label() <<
" ...\n";
747 OperatorFree( operator_ptrtype F, prec_ptrtype Prec, list_type options )
750 M_Solver( new solver_type() ),
752 M_hasApply( F->hasApply() ),
753 M_useTranspose( F->UseTranspose() ),
754 M_label( F->Label() ),
757 DVLOG(2) <<
"Create operator " << Label() <<
" ...\n";
758 M_tol = options.get(
"tol", 1e-7 );
759 M_maxiter = options.get(
"max_iter", 100 );
761 M_Solver->setOptions( options );
765 OperatorFree(
const OperatorFree& tc )
768 M_hasInverse( tc.M_hasInverse ),
769 M_hasApply( tc.M_hasApply ),
770 M_useTranspose( tc.M_useTranspose ),
771 M_Solver( tc.M_Solver ),
772 M_label( tc.M_label ),
775 DVLOG(2) <<
"Copy operator " << Label() <<
" ...\n";
779 bool hasInverse()
const
784 bool hasApply()
const
790 int Apply(
const Epetra_MultiVector & X, Epetra_MultiVector & Y )
const
792 DVLOG(2) <<
"Apply operator " << Label() <<
"\n";
794 DVLOG(2) <<
"Finished Apply operator " << Label() <<
"\n";
799 int ApplyInverse (
const Epetra_MultiVector& X, Epetra_MultiVector& Y )
const
801 DVLOG(2) <<
"ApplyInverse operator " << Label() <<
"\n";
803 FEELPP_ASSERT( hasInverse() ).error(
"This operator cannot be inverted." );
805 std::pair<unsigned int, real_type> result = M_Solver->solve( M_op, M_Prec, Y, X, M_tol, M_maxiter );
807 DVLOG(2) <<
"Finished ApplyInverse operator " << Label() <<
"\n";
808 return !hasInverse();
812 int SetUseTranspose(
bool UseTranspose = 0 )
814 M_useTranspose = UseTranspose;
819 double NormInf()
const
824 const char * Label ()
const
826 return( M_label.c_str() );
829 bool UseTranspose()
const
831 return M_useTranspose;
834 bool HasNormInf ()
const
839 const Epetra_Comm & Comm()
const
841 return( M_op->OperatorDomainMap().Comm() );
844 const Epetra_Map & OperatorDomainMap()
const
846 return( M_op->OperatorDomainMap() );
849 const Epetra_Map & OperatorRangeMap()
const
851 return( M_op->OperatorRangeMap() );
856 DVLOG(2) <<
"Destroyed operator: " << Label() <<
" ...\n";
861 operator_ptrtype M_op;
863 solver_ptrtype M_Solver;
865 bool M_hasInverse, M_hasApply, M_useTranspose;
884 #endif // FEELPP_HAS_TRILINOS_EPETRA