30 #ifndef __trilinos_linear_solver_H
31 #define __trilinos_linear_solver_H 1
35 #include <feel/feelalg/solverlinear.hpp>
39 #if defined( FEELPP_HAS_TRILINOS_AZTECOO )
40 #undef PACKAGE_BUGREPORT
43 #undef PACKAGE_TARNAME
44 #undef PACKAGE_VERSION
45 #include <AztecOO_config.h>
47 #undef PACKAGE_BUGREPORT
50 #undef PACKAGE_TARNAME
51 #undef PACKAGE_VERSION
56 class BackendTrilinos;
59 template<
typename T =
double >
60 class SolverLinearTrilinos :
public SolverLinear<T>
62 typedef SolverLinear<T> super;
66 typedef typename super::value_type value_type;
67 typedef typename super::real_type real_type;
71 typedef MatrixEpetra sparse_matrix_type;
72 typedef VectorEpetra<value_type> vector_type;
73 typedef boost::shared_ptr<sparse_matrix_type> sparse_matrix_ptrtype;
74 typedef boost::shared_ptr<vector_type> vector_ptrtype;
76 typedef Teuchos::ParameterList list_type;
80 SolverLinearTrilinos()
84 SolverLinearTrilinos( po::variables_map
const& vm )
90 ~SolverLinearTrilinos()
97 if ( this->initialized() )
99 this->setInitialized(
false );
105 if ( !this->initialized() )
107 this->setInitialized(
true );
109 M_Solver.SetParameters( M_List,
true );
113 double status[AZ_STATUS_SIZE];
114 M_Solver.GetAllAztecStatus( status );
118 void setOptions( list_type _list )
123 list_type& getOptions()
129 boost::tuple<bool,unsigned int, real_type>
130 solve ( MatrixSparse<T>
const& matrix,
131 Vector<T> & solution,
132 Vector<T>
const& rhs,
134 const unsigned int m_its,
135 bool transpose =
false )
137 DVLOG(2) <<
"Matrix solver...\n";
141 setUserOperator( matrix );
143 M_Solver.SetParameters( M_List,
true );
144 M_Solver.Iterate( m_its, tol );
148 return boost::make_tuple(
true, M_Solver.NumIters(), M_Solver.TrueResidual() );
152 boost::tuple<bool,unsigned int, real_type>
153 solve ( MatrixSparse<T>
const& matrix,
154 MatrixSparse<T>
const& preconditioner,
156 Vector<T>
const& rhs,
158 const unsigned int m_its,
159 bool transpose =
false )
161 DVLOG(2) <<
"Matrix solver with preconditioner...\n";
165 setUserOperator( matrix );
167 M_Solver.SetParameters( M_List,
true );
168 M_Solver.Iterate( m_its, tol );
172 return boost::make_tuple(
true, M_Solver.NumIters(), M_Solver.TrueResidual() );
176 std::pair<unsigned int, real_type>
177 solve ( boost::shared_ptr<OperatorMatrix>
const& op,
178 Vector<T> & solution,
179 Vector<T>
const& rhs,
181 const unsigned int m_its,
182 bool transpose =
false )
184 DVLOG(2) <<
"Operator solver...\n";
188 setUserOperator( op );
190 M_Solver.SetParameters( M_List,
true );
191 M_Solver.Iterate( m_its, tol );
193 return std::make_pair( M_Solver.NumIters(), M_Solver.TrueResidual() );
196 template<
typename operator_type >
197 std::pair<unsigned int, real_type>
198 solve ( operator_type
const& op,
199 Epetra_MultiVector & solution,
200 Epetra_MultiVector
const& rhs,
202 const unsigned int m_its )
204 DVLOG(2) <<
"Operator solver...\n";
208 setUserOperator( op );
210 M_Solver.SetParameters( M_List,
true );
211 M_Solver.Iterate( m_its, tol );
213 return std::make_pair( M_Solver.NumIters(), M_Solver.TrueResidual() );
217 template<
typename op1_type,
typename op2_type >
218 std::pair<unsigned int, real_type>
219 solve ( op1_type
const& op1,
221 Epetra_MultiVector & solution,
222 Epetra_MultiVector
const& rhs,
224 const unsigned int m_its )
226 DVLOG(2) <<
"Operator solver with preconditioner...\n";
230 setUserOperator( op1 );
232 if ( op2.get() != 0 )
233 M_Solver.SetPrecOperator( &( *op2 ) );
235 M_Solver.SetParameters( M_List,
true );
236 M_Solver.Iterate( m_its, tol );
238 return std::make_pair( M_Solver.NumIters(), M_Solver.TrueResidual() );
248 template<
typename operator_type >
249 void setUserOperator( operator_type
const& D )
251 setUserOperator( mpl::bool_<
252 mpl::or_< boost::is_same< operator_type, MatrixSparse<T> >,
253 boost::is_same< operator_type, MatrixEpetra > >::value >(),
257 void setUserOperator( mpl::bool_<true>, MatrixSparse<T>
const& A )
259 DVLOG(2) <<
"Set matrix operator...\n";
260 sparse_matrix_type* A_ptr =
const_cast<sparse_matrix_type *
>(
dynamic_cast< sparse_matrix_type const*
>( &A ) );
262 M_Solver.SetUserMatrix( &A_ptr->mat() );
267 template<
typename operator_type >
268 void setUserOperator( mpl::bool_<false>, operator_type
const& A )
270 DVLOG(2) <<
"Set epetra operator...\n";
272 M_Solver.SetUserOperator( &( *A ) );
276 void setRHS( Vector<T>
const& x )
278 VectorEpetra<T>* aux =
const_cast< VectorEpetra<T>*
>(
dynamic_cast<VectorEpetra<T> const*
>( &x ) );
280 M_Solver.SetRHS( &aux->vec() );
283 void setLHS( Vector<T>& x )
285 VectorEpetra<T>* aux =
dynamic_cast< VectorEpetra<T>*
>( &x );
287 M_Solver.SetLHS( &aux->vec() );
290 void setRHS( Epetra_MultiVector
const& x )
292 Epetra_MultiVector* aux =
const_cast< Epetra_MultiVector*
>(
dynamic_cast<Epetra_MultiVector const*
>( &x ) );
294 M_Solver.SetRHS( &( *aux ) );
297 void setLHS( Epetra_MultiVector& x )
299 M_Solver.SetLHS( &x );