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
solverlineartrilinos.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@epfl.ch>
6  Date: 2006-03-06
7 
8  Copyright (C) 2006 EPFL
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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
23 */
30 #ifndef __trilinos_linear_solver_H
31 #define __trilinos_linear_solver_H 1
32 
33 
34 
35 #include <feel/feelalg/solverlinear.hpp>
36 
38 
39 #if defined( FEELPP_HAS_TRILINOS_AZTECOO )
40 #undef PACKAGE_BUGREPORT
41 #undef PACKAGE_NAME
42 #undef PACKAGE_STRING
43 #undef PACKAGE_TARNAME
44 #undef PACKAGE_VERSION
45 #include <AztecOO_config.h>
46 #include <AztecOO.h>
47 #undef PACKAGE_BUGREPORT
48 #undef PACKAGE_NAME
49 #undef PACKAGE_STRING
50 #undef PACKAGE_TARNAME
51 #undef PACKAGE_VERSION
52 
53 namespace Feel
54 {
55 
56 class BackendTrilinos;
57 class OperatorMatrix;
58 
59 template< typename T = double >
60 class SolverLinearTrilinos : public SolverLinear<T>
61 {
62  typedef SolverLinear<T> super;
63 
64 public:
65 
66  typedef typename super::value_type value_type;
67  typedef typename super::real_type real_type;
68 
69  //typedef BackendTrilinos backend_type;
70 
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;
75 
76  typedef Teuchos::ParameterList list_type;
77  //typedef typename backend_type::list_type list_type;
78 
79 
80  SolverLinearTrilinos()
81  {
82  }
83 
84  SolverLinearTrilinos( po::variables_map const& vm )
85  :
86  super( vm )
87  {
88  }
89 
90  ~SolverLinearTrilinos()
91  {
92  this->clear();
93  }
94 
95  void clear()
96  {
97  if ( this->initialized() )
98  {
99  this->setInitialized( false );
100  }
101  }
102 
103  void init()
104  {
105  if ( !this->initialized() )
106  {
107  this->setInitialized( true );
108 
109  M_Solver.SetParameters( M_List, true );
110 
111  // AztecOO defined a certain number of output parameters, and store them
112  // in a double vector called status.
113  double status[AZ_STATUS_SIZE];
114  M_Solver.GetAllAztecStatus( status );
115  }
116  }
117 
118  void setOptions( list_type _list )
119  {
120  M_List = _list;
121  }
122 
123  list_type& getOptions()
124  {
125  return M_List;
126  }
127 
128  //std::pair<unsigned int, real_type>
129  boost::tuple<bool,unsigned int, real_type>
130  solve ( MatrixSparse<T> const& matrix,
131  Vector<T> & solution,
132  Vector<T> const& rhs,
133  const double tol,
134  const unsigned int m_its,
135  bool transpose = false )
136  {
137  DVLOG(2) << "Matrix solver...\n";
138 
139  setRHS( rhs );
140  setLHS( solution );
141  setUserOperator( matrix );
142 
143  M_Solver.SetParameters( M_List, true );
144  M_Solver.Iterate( m_its, tol );
145 
146  //return std::make_pair( M_Solver.NumIters(), M_Solver.TrueResidual() );
147 #warning todo!
148  return boost::make_tuple( true, M_Solver.NumIters(), M_Solver.TrueResidual() );
149  }
150 
151  //std::pair<unsigned int, real_type>
152  boost::tuple<bool,unsigned int, real_type>
153  solve ( MatrixSparse<T> const& matrix,
154  MatrixSparse<T> const& preconditioner,
155  Vector<T>& solution,
156  Vector<T> const& rhs,
157  const double tol,
158  const unsigned int m_its,
159  bool transpose = false )
160  {
161  DVLOG(2) << "Matrix solver with preconditioner...\n";
162 
163  setRHS( rhs );
164  setLHS( solution );
165  setUserOperator( matrix );
166 
167  M_Solver.SetParameters( M_List, true );
168  M_Solver.Iterate( m_its, tol );
169 
170  //return std::make_pair( M_Solver.NumIters(), M_Solver.TrueResidual() );
171 #warning todo!
172  return boost::make_tuple( true, M_Solver.NumIters(), M_Solver.TrueResidual() );
173 
174  }
175 
176  std::pair<unsigned int, real_type>
177  solve ( boost::shared_ptr<OperatorMatrix> const& op,
178  Vector<T> & solution,
179  Vector<T> const& rhs,
180  const double tol,
181  const unsigned int m_its,
182  bool transpose = false )
183  {
184  DVLOG(2) << "Operator solver...\n";
185 
186  setRHS( rhs );
187  setLHS( solution );
188  setUserOperator( op );
189 
190  M_Solver.SetParameters( M_List, true );
191  M_Solver.Iterate( m_its, tol );
192 
193  return std::make_pair( M_Solver.NumIters(), M_Solver.TrueResidual() );
194  }
195 
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,
201  const double tol,
202  const unsigned int m_its )
203  {
204  DVLOG(2) << "Operator solver...\n";
205 
206  setRHS( rhs );
207  setLHS( solution );
208  setUserOperator( op );
209 
210  M_Solver.SetParameters( M_List, true );
211  M_Solver.Iterate( m_its, tol );
212 
213  return std::make_pair( M_Solver.NumIters(), M_Solver.TrueResidual() );
214  }
215 
216 
217  template< typename op1_type, typename op2_type >
218  std::pair<unsigned int, real_type>
219  solve ( op1_type const& op1,
220  op2_type const& op2,
221  Epetra_MultiVector & solution,
222  Epetra_MultiVector const& rhs,
223  const double tol,
224  const unsigned int m_its )
225  {
226  DVLOG(2) << "Operator solver with preconditioner...\n";
227 
228  setRHS( rhs );
229  setLHS( solution );
230  setUserOperator( op1 );
231 
232  if ( op2.get() != 0 )
233  M_Solver.SetPrecOperator( &( *op2 ) );
234 
235  M_Solver.SetParameters( M_List, true );
236  M_Solver.Iterate( m_its, tol );
237 
238  return std::make_pair( M_Solver.NumIters(), M_Solver.TrueResidual() );
239  }
240 
241 private:
242 
243  AztecOO M_Solver;
244 
245  list_type M_List;
246 
247  // sets the user specified operator
248  template< typename operator_type >
249  void setUserOperator( operator_type const& D )
250  {
251  setUserOperator( mpl::bool_<
252  mpl::or_< boost::is_same< operator_type, MatrixSparse<T> >,
253  boost::is_same< operator_type, MatrixEpetra > >::value >(),
254  D );
255  }
256 
257  void setUserOperator( mpl::bool_<true>, MatrixSparse<T> const& A )
258  {
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 ) );
261 
262  M_Solver.SetUserMatrix( &A_ptr->mat() );
263  }
264 
265 
266  // case the operator is passed as an Epetra_Operator
267  template< typename operator_type >
268  void setUserOperator( mpl::bool_<false>, operator_type const& A )
269  {
270  DVLOG(2) << "Set epetra operator...\n";
271 
272  M_Solver.SetUserOperator( &( *A ) );
273  }
274 
275 
276  void setRHS( Vector<T> const& x )
277  {
278  VectorEpetra<T>* aux = const_cast< VectorEpetra<T>* >( dynamic_cast<VectorEpetra<T> const*>( &x ) );
279 
280  M_Solver.SetRHS( &aux->vec() );
281  }
282 
283  void setLHS( Vector<T>& x )
284  {
285  VectorEpetra<T>* aux = dynamic_cast< VectorEpetra<T>* >( &x );
286 
287  M_Solver.SetLHS( &aux->vec() );
288  }
289 
290  void setRHS( Epetra_MultiVector const& x )
291  {
292  Epetra_MultiVector* aux = const_cast< Epetra_MultiVector* >( dynamic_cast<Epetra_MultiVector const*>( &x ) );
293 
294  M_Solver.SetRHS( &( *aux ) );
295  }
296 
297  void setLHS( Epetra_MultiVector& x )
298  {
299  M_Solver.SetLHS( &x );
300  }
301 };
302 
303 
304 } // Feel
305 
306 #endif /* FEELPP_HAS_TRILINOS_AZTECOO */
307 #endif /* __trilinos_linear_solver_H */
308 
309 

Generated on Sun Oct 20 2013 08:25:05 for Feel++ by doxygen 1.8.4