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
backendtrilinos.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): Christoph Winkelmann <christoph.winkelmann@epfl.ch>
6  Christophe Prud'homme <christophe.prudhomme@feelpp.org>
7  Date: 2007-01-18
8 
9  Copyright (C) 2007 EPFL
10  Copyright (C) 2008-2011 Université Joseph Fourier (Grenoble 1)
11 
12 
13  This library is free software; you can redistribute it and/or
14  modify it under the terms of the GNU Lesser General Public
15  License as published by the Free Software Foundation; either
16  version 3.0 of the License, or (at your option) any later version.
17 
18  This library is distributed in the hope that it will be useful,
19  but WITHOUT ANY WARRANTY; without even the implied warranty of
20  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
21  Lesser General Public License for more details.
22 
23  You should have received a copy of the GNU Lesser General Public
24  License along with this library; if not, write to the Free Software
25  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
26 */
34 #ifndef _BACKENDTRILINOS_HPP_
35 #define _BACKENDTRILINOS_HPP_
36 
37 
38 #include <boost/program_options/variables_map.hpp>
39 #include <feel/feelconfig.h>
40 
41 #if defined ( FEELPP_HAS_TRILINOS_EPETRA )
42 #undef PACKAGE_BUGREPORT
43 #undef PACKAGE_NAME
44 #undef PACKAGE_STRING
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>
54 #include <feel/feelalg/backend.hpp>
55 
56 
57 #include <Teuchos_ParameterList.hpp>
58 
59 namespace Feel
60 {
61 namespace po = boost::program_options;
62 
63 // forward declararion
64 template<class T> class Backend;
65 class OperatorMatrix;
66 
72 class BackendTrilinos : public Backend<double>
73 {
74  typedef Backend<double> super;
75 public:
76 
77  // -- TYPEDEFS --
78  typedef super::value_type value_type;
79 
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;
84 
85  typedef sparse_matrix_type::graph_type graph_type;
86  typedef sparse_matrix_type::graph_ptrtype graph_ptrtype;
87 
88  typedef MatrixEpetra epetra_sparse_matrix_type;
89  typedef VectorEpetra<value_type> epetra_vector_type;
90  typedef Epetra_Operator operator_type;
91 
92  typedef OperatorMatrix op_mat_type;
93  typedef boost::shared_ptr<OperatorMatrix> op_mat_ptrtype;
94 
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;
98 
99  typedef super::solve_return_type solve_return_type;
100  typedef super::nl_solve_return_type nl_solve_return_type;
101 
102  typedef std::map< size_type, std::vector< size_type > > procdist_type;
103 
104  typedef Teuchos::ParameterList list_type;
105 
106 
107  // -- CONSTRUCTOR --
108  BackendTrilinos( WorldComm const& worldComm=Environment::worldComm())
109  :
110  super( worldComm ),
111  M_options(),
112  M_prec_type( "" ),
113  M_Prec()
114  {
115  set_maxiter( 1000 );
116  set_tol( 1e-10 );
117  }
118 
119 
120  BackendTrilinos( po::variables_map const& vm, std::string const& prefix = "", WorldComm const& worldComm=Environment::worldComm() );
121 
122  BackendTrilinos( const BackendTrilinos& tc );
123 
124 
125  // -- FACTORY METHODS --
126  static Epetra_Map epetraMap( DataMap const& dmap )
127  {
128  std::vector<int> e( dmap.nMyElements() );
129  std::copy( dmap.myGlobalElements().begin(),
130  dmap.myGlobalElements().end(),
131  e.begin() );
132  return Epetra_Map( -1, dmap.nMyElements(), e.data(), 0, Epetra_MpiComm( dmap.comm() ) );
133  }
134 
135 
136  static Epetra_Map epetraMapStatic( DataMap const& dmap )
137  {
138  return Epetra_Map( dmap.nGlobalElements(), dmap.nMyElements(), 0, Epetra_MpiComm( dmap.comm() ) );
139  }
140 
141 
142  template<typename DomainSpace, typename DualImageSpace>
143  sparse_matrix_ptrtype newMatrix( DomainSpace const& Xh,
144  DualImageSpace const& Yh,
145  size_type matrix_properties = NON_HERMITIAN,
146  bool init = true )
147  {
148  return newMatrix( Xh->map(), Yh->map(), matrix_properties );
149  }
150 
151  sparse_matrix_ptrtype
152  newMatrix( const size_type m,
153  const size_type n,
154  const size_type m_l,
155  const size_type n_l,
156  const size_type nnz=30,
157  const size_type noz=10,
158  size_type matrix_properties = NON_HERMITIAN )
159  {
160  sparse_matrix_ptrtype mat( new epetra_sparse_matrix_type( m,n,m_l,n_l,nnz,noz ) );
161  mat->setMatrixProperties( matrix_properties );
162  return mat;
163  }
164 
165  sparse_matrix_ptrtype
166  newMatrix( const size_type m,
167  const size_type n,
168  const size_type m_l,
169  const size_type n_l,
170  graph_ptrtype const & graph,
171  size_type matrix_properties = NON_HERMITIAN )
172  {
173  sparse_matrix_ptrtype mat( new epetra_sparse_matrix_type( m,n,m_l,n_l,30,10 ) );
174  mat->setMatrixProperties( matrix_properties );
175  return mat;
176  }
177 
178  sparse_matrix_ptrtype newMatrix( DataMap const& domainmap,
179  DataMap const& imagemap,
180  size_type matrix_properties = NON_HERMITIAN,
181  bool init = true )
182  {
183  Epetra_Map erowmap = BackendTrilinos::epetraMap( imagemap );
184  Epetra_Map ecolmap = BackendTrilinos::epetraMapStatic( domainmap );
185  Epetra_Map edomainmap = BackendTrilinos::epetraMap( imagemap );
186 
187  //std::cout << "Rowmap: " << erowmap << "\n";
188  //std::cout << "Colmap: " << ecolmap << "\n";
189  //std::cout << "Domainmap: " << edomainmap << "\n";
190 
191  //std::cout << "Is matrix rectangular? " << !erowmap.SameAs( ecolmap ) << "\n";
192 
193  if ( !erowmap.SameAs( ecolmap ) )
194  {
195  Epetra_Map eimagemap = BackendTrilinos::epetraMap( domainmap );
196  //std::cout << "Imagemap: " << eimagemap << "\n";
197  auto A= sparse_matrix_ptrtype( new epetra_sparse_matrix_type( erowmap, ecolmap, edomainmap, eimagemap ) );
198  A->setMatrixProperties( matrix_properties );
199  return A;
200  }
201 
202  else
203  {
204  auto A= sparse_matrix_ptrtype( new epetra_sparse_matrix_type( erowmap, ecolmap ) );
205  A->setMatrixProperties( matrix_properties );
206  return A;
207  }
208  }
209 
210  sparse_matrix_ptrtype
211  newZeroMatrix( const size_type m,
212  const size_type n,
213  const size_type m_l,
214  const size_type n_l )
215  {
216  sparse_matrix_ptrtype mat( new epetra_sparse_matrix_type( m,n,m_l,n_l,0,0 ) );
217  return mat;
218 
219  }
220  sparse_matrix_ptrtype
221  newZeroMatrix( DataMap const& domainmap, DataMap const& imagemap )
222  {
223  Epetra_Map erowmap = BackendTrilinos::epetraMap( imagemap );
224  Epetra_Map ecolmap = BackendTrilinos::epetraMapStatic( domainmap );
225 
226  auto A= sparse_matrix_ptrtype( new epetra_sparse_matrix_type( erowmap, ecolmap ) );
227  //A->setMatrixProperties( matrix_properties );
228  return A;
229  }
230 
231 
232  template<typename SpaceT>
233  vector_ptrtype newVector( SpaceT const& space )
234  {
235  return newVector( space->map() );
236  }
237  vector_ptrtype newVector( DataMap const& domainmap )
238  {
239  Epetra_Map emap = BackendTrilinos::epetraMap( domainmap );
240 
241  return vector_ptrtype( new epetra_vector_type( emap ) );
242  }
243 
244  vector_ptrtype newVector( const size_type n, const size_type n_local )
245  {
246  return vector_ptrtype( new epetra_vector_type( /*n, n_local*/ ) );
247 #warning to fix!
248  }
249 
250 
251  static operator_ptrtype IfpackPrec( sparse_matrix_ptrtype const& M, list_type options, std::string precType = "Amesos" )
252  {
253  PreconditionerIfpack P( options, precType );
254 
255  P.buildPreconditioner( M );
256 
257  return P.getPrec();
258  }
259 
260 
261  static operator_ptrtype MLPrec( sparse_matrix_ptrtype& M, list_type options )
262  {
263  PreconditionerML P( options );
264 
265  P.buildPreconditioner( M );
266 
267  return P.getPrec();
268  }
269 
270 #if 0
271  template< typename element_type >
272  static void Epetra2Ublas( vector_ptrtype const& u, element_type& x )
273  {
274  epetra_vector_ptrtype const& _v( dynamic_cast<epetra_vector_ptrtype const&>( u ) );
275  Epetra_Map v_map( _v->Map() );
276 
277  vector_type v = *u;
278 
279  //DVLOG(2) << "Initial EpetraVector " << v << "\n";
280 
281  const size_type L = v.localSize();
282 
283  for ( size_type i=0; i<L; i++ )
284  {
285  DVLOG(2) << "x(" << x.firstLocalIndex() + i << ")="
286  << "v[" << v_map.GID( i ) << "] = "
287  << v( i ) << "\n";
288 
289  x( x.firstLocalIndex() + i ) = v( i );
290  }
291 
292  DVLOG(2) << "Epetra2Ublas:" << x << "\n";
293  }
294 
295 
296  template< typename element_type >
297  static void Ublas2Epetra( element_type const& x, vector_ptrtype& v )
298  {
299  epetra_vector_type& _v( dynamic_cast<epetra_vector_type&>( *v ) );
300  Epetra_Map v_map( _v.Map() );
301 
302  DVLOG(2) << "Local size of ublas vector" << x.localSize() << "\n";
303  DVLOG(2) << "Local size of epetra vector" << v->localSize() << "\n";
304 
305  const size_type L = v->localSize();
306 
307  for ( size_type i=0; i<L; i++ )
308  {
309  DVLOG(2) << "v[" << v_map.GID( i ) << "] = "
310  << "x(" << x.firstLocalIndex() + i << ")="
311  << x( x.firstLocalIndex() + i ) << "\n";
312 
313  v->set( v_map.GID( i ), x( x.firstLocalIndex() + i ) );
314  }
315  }
316 #endif
317 
318  template< int index, typename spaceT >
319  static Epetra_MultiVector getComponent( spaceT const& Xh, Epetra_MultiVector const& sol )
320  {
321  Epetra_Map componentMap ( epetraMap( Xh->template functionSpace<index>()->map() ) );
322  Epetra_Map globalMap ( epetraMap( Xh->map() ) );
323 
324  //DVLOG(2) << "Component map: " << componentMap << "\n";
325 
326  Epetra_MultiVector component( componentMap, 1 );
327 
328  int Length = component.MyLength();
329 
330  int shift = Xh->nDofStart( index );
331 
332  for ( int i=0; i < Length; i++ )
333  {
334  int compGlobalID = componentMap.GID( i );
335 
336  if ( compGlobalID >= 0 )
337  {
338  int compLocalID = componentMap.LID( compGlobalID );
339 
340  int localID = globalMap.LID( compGlobalID+shift );
341  // int globalID = globalMap.GID(localID);
342 
343  DVLOG(2) << "[MyBackend] Copy entry sol[" << localID << "]=" << sol[0][localID]
344  << " to component[" << compLocalID << "]\n";
345 
346  component[0][compLocalID] = sol[0][localID];
347 
348  DVLOG(2) << component[0][compLocalID] << "\n";
349  }
350  }
351 
352  return component;
353  }
354 
355 
356  template< int index, typename spaceT >
357  static void UpdateComponent( spaceT const& Xh, Epetra_MultiVector& sol, Epetra_MultiVector& comp )
358  {
359  Epetra_Map componentMap ( epetraMap( Xh->template functionSpace<index>()->map() ) );
360  Epetra_Map globalMap ( epetraMap( Xh->map() ) );
361 
362  int shift = Xh->nDofStart( index );
363 
364  int Length = comp.MyLength();
365 
366  for ( int i=0; i < Length; i++ )
367  {
368  int compGlobalID = componentMap.GID( i );
369 
370  if ( compGlobalID >= 0 )
371  {
372  int compLocalID = componentMap.LID( compGlobalID );
373 
374  int localID = globalMap.LID( compGlobalID+shift );
375  // int globalID = globalMap.GID(localID);
376 
377  DVLOG(2) << "Copy entry component[" << compLocalID << "] to sol[" << localID << "]="
378  << sol[0][localID]
379  << "]\n";
380 
381  sol[0][localID] = comp[0][compLocalID] ;
382 
383  DVLOG(2) << comp[0][compLocalID] << "\n";
384  }
385  }
386  }
387 
388 
389 
390  // -- SETTING OF OPTIONS --
391  void set_maxiter ( int maxiter );
392 
393  void set_tol ( double tol );
394 
395  void set_drop ( double drop );
396 
397  void set_overlap ( int overlap );
398 
399  void set_restart ( int restart );
400 
401  void set_fillin ( int fillin );
402 
403  void set_order ( int order );
404 
405  void set_overlap_type( std::string str );
406 
407  void set_residual ( std::string str );
408 
409  void set_localparts ( int localparts );
410 
411  void set_solver( std::string str );
412 
413  void set_prec_local_solver( std::string str );
414 
415  void set_prec_local_solver_type( std::string str );
416 
417  void set_prec( std::string str );
418 
419  void set_verbose( const std::string verb="none" );
420 
421  void set_options( list_type opts );
422 
423  list_type get_options();
424 
425  void set_symmetric( bool /*is_sym*/ ) {}
426 
427 
428  // -- LINEAR ALGEBRA INTERFACE --
429  static void
430  applyMatrix( sparse_matrix_type const& A,
431  const Vector<value_type>& x,
432  vector_type& b );
433 
434  real_type dot( const vector_type& f, const vector_type& x ) const;
435 
436  void prod( sparse_matrix_type const& A,
437  vector_type const& x,
438  vector_type& b ) const
439  {
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() );
444  }
445 
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 );
450 
451  bool converged()
452  {
453  return true;
454  }
455 
456 private:
457 
458  mpi::communicator M_comm;
459 
460  list_type M_options;
461 
462  std::string M_prec_type, M_local_solver, M_local_solver_type;
463 
464  operator_ptrtype M_Prec;
465 
466 }; // class BackendTrilinos
467 
473 po::options_description backendtrilinos_options( std::string const& prefix );
474 
475 
476 } // Feel
477 
478 #endif // FEELPP_HAS_TRILINOS_EPETRA
479 
480 
481 #endif /* _BACKENDTRILINOS_HPP_ */
482 

Generated on Fri Oct 25 2013 14:24:04 for Feel++ by doxygen 1.8.4