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
backendeigen.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@feelpp.org>
6  Christophe Prud'homme <christophe.prudhomme@feelpp.org>
7  Date: 2012-06-20
8 
9  Copyright (C) 2012 Université Joseph Fourier (Grenoble I)
10 
11  This library is free software; you can redistribute it and/or
12  modify it under the terms of the GNU Lesser General Public
13  License as published by the Free Software Foundation; either
14  version 3.0 of the License, or (at your option) any later version.
15 
16  This library is distributed in the hope that it will be useful,
17  but WITHOUT ANY WARRANTY; without even the implied warranty of
18  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19  Lesser General Public License for more details.
20 
21  You should have received a copy of the GNU Lesser General Public
22  License along with this library; if not, write to the Free Software
23  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
24 */
31 #ifndef _BACKENDEIGEN_HPP_
32 #define _BACKENDEIGEN_HPP_
33 
34 #include <boost/program_options/variables_map.hpp>
36 #undef MatType
39 #include <feel/feelalg/matrixeigendense.hpp>
40 #include <feel/feelalg/matrixeigensparse.hpp>
42 
43 
44 
45 #include <Eigen/Dense>
46 #include <Eigen/Sparse>
47 
48 #include <feel/feelalg/backend.hpp>
49 namespace Feel
50 {
51 namespace po = boost::program_options;
52 
53 
59 template<typename T, int _Options = 0>
60 class BackendEigen : public Backend<T>
61 {
62  typedef Backend<T> super;
63 public:
64 
65  // -- TYPEDEFS --
66  typedef typename super::value_type value_type;
67  typedef typename super::real_type real_type;
68 
69  static const bool IsDense = (_Options == 1);
70  static const bool IsSparse = (_Options == 0);
71 
72  /* matrix */
74  typedef typename super::sparse_matrix_ptrtype sparse_matrix_ptrtype;
75  typedef typename mpl::if_<mpl::bool_<IsDense>,
76  mpl::identity<MatrixEigenDense<value_type> >,
77  mpl::identity<MatrixEigenSparse<value_type> > >::type::type eigen_sparse_matrix_type;
78  typedef boost::shared_ptr<eigen_sparse_matrix_type> eigen_sparse_matrix_ptrtype;
79 
81  typedef typename sparse_matrix_type::graph_ptrtype graph_ptrtype;
82 
83  /* vector */
84  typedef typename super::vector_type vector_type;
85  typedef typename super::vector_ptrtype vector_ptrtype;
87  typedef boost::shared_ptr<vector_type> eigen_vector_ptrtype;
88 
89  typedef typename super::solve_return_type solve_return_type;
90  typedef typename super::nl_solve_return_type nl_solve_return_type;
91 
92  typedef typename super::datamap_type datamap_type;
93  typedef typename super::datamap_ptrtype datamap_ptrtype;
94 
95  // -- CONSTRUCTOR --
96  BackendEigen();
97  BackendEigen( WorldComm const& worldComm=Environment::worldComm() );
98  BackendEigen( po::variables_map const& vm, std::string const& prefix = "",
99  WorldComm const& worldComm=Environment::worldComm() );
100 
101  // -- FACTORY METHODS --
102  template<typename DomainSpace, typename DualImageSpace>
103  static sparse_matrix_ptrtype newMatrix( boost::shared_ptr<DomainSpace> const& space1,
104  boost::shared_ptr<DualImageSpace> const& space2,
105  size_type matrix_properties = NON_HERMITIAN )
106  {
107  Context ctx( matrix_properties );
108  //if ( ctx.test( DENSE ) )
109  {
110  auto A= sparse_matrix_ptrtype( new eigen_sparse_matrix_type( space1->nDof(), space2->nDof() ) );
111  A->setMatrixProperties( matrix_properties );
112  return A;
113  }
114  }
115 
116  sparse_matrix_ptrtype
118  const size_type n,
119  const size_type m_l,
120  const size_type n_l,
121  const size_type nnz=30,
122  const size_type noz=10,
123  size_type matrix_properties = NON_HERMITIAN )
124  {
125  sparse_matrix_ptrtype mat( new eigen_sparse_matrix_type( m,n ) );
126  mat->setMatrixProperties( matrix_properties );
127  return mat;
128  }
129 
130  sparse_matrix_ptrtype
132  const size_type n,
133  const size_type m_l,
134  const size_type n_l,
135  graph_ptrtype const & graph,
136  size_type matrix_properties = NON_HERMITIAN )
137  {
138  sparse_matrix_ptrtype mat( new eigen_sparse_matrix_type( m,n ) );
139  mat->setMatrixProperties( matrix_properties );
140  return mat;
141  }
142 
143  sparse_matrix_ptrtype
144  newMatrix( datamap_ptrtype const& d1, datamap_ptrtype const& d2, size_type matrix_properties = NON_HERMITIAN, bool init = true )
145  {
146  auto A = sparse_matrix_ptrtype( new eigen_sparse_matrix_type( d1->nGlobalElements(), d2->nGlobalElements() ) );
147  A->setMatrixProperties( matrix_properties );
148  return A;
149  }
150 
151  sparse_matrix_ptrtype
153  const size_type n,
154  const size_type m_l,
155  const size_type n_l )
156  {
157  auto A = sparse_matrix_ptrtype( new eigen_sparse_matrix_type( m, n ) );
158  //A->setMatrixProperties( matrix_properties );
159  return A;
160  }
161 
162 
163  sparse_matrix_ptrtype
164  newZeroMatrix( datamap_ptrtype const& d1, datamap_ptrtype const& d2 )
165  {
166  auto A = sparse_matrix_ptrtype( new eigen_sparse_matrix_type( d1->nGlobalElements(), d2->nGlobalElements() ) );
167  //A->setMatrixProperties( matrix_properties );
168  return A;
169  }
170 
171  template<typename SpaceT>
172  static vector_ptrtype newVector( boost::shared_ptr<SpaceT> const& space )
173  {
174  return vector_ptrtype( new eigen_vector_type( space->nDof() ) );
175  }
176 
177  template<typename SpaceT>
178  static vector_ptrtype newVector( SpaceT const& space )
179  {
180  return vector_ptrtype( new eigen_vector_type( space.nDof() ) );
181  }
182 
183  vector_ptrtype
184  newVector( datamap_ptrtype const& d )
185  {
186  return vector_ptrtype( new eigen_vector_type( d->nGlobalElements() ) );
187  }
188 
189  vector_ptrtype newVector( const size_type n, const size_type n_local )
190  {
191  return vector_ptrtype( new eigen_vector_type( n ) );
192  }
193 
194 
195  // -- LINEAR ALGEBRA INTERFACE --
196  void prod( sparse_matrix_type const& A,
197  vector_type const& x,
198  vector_type& b ) const
199  {
200  eigen_sparse_matrix_type const& _A = dynamic_cast<eigen_sparse_matrix_type const&>( A );
201  eigen_vector_type const& _x = dynamic_cast<eigen_vector_type const&>( x );
202  eigen_vector_type& _b = dynamic_cast<eigen_vector_type&>( b );
203  _b.vec() = _A.mat()*_x.vec();
204  }
205 
206  solve_return_type solve( sparse_matrix_type const& A,
207  vector_type& x,
208  const vector_type& b )
209  {
210  return this->solve( A, x, b, mpl::bool_<IsDense>() );
211  }
212 
213  solve_return_type solve( sparse_matrix_ptrtype const& A,
214  vector_ptrtype& x,
215  const vector_ptrtype& b )
216  {
217  return this->solve( *A, *x, *b );
218  }
219 
220  solve_return_type solve( sparse_matrix_ptrtype const& A,
221  sparse_matrix_ptrtype const& P,
222  vector_ptrtype& x,
223  const vector_ptrtype& b )
224  {
225  return this->solve( *A, *x, *b );
226  }
227 
228 
229  real_type dot( const vector_type& f,
230  const vector_type& x ) const
231  {
232  eigen_vector_type const& _f = dynamic_cast<eigen_vector_type const&>( f );
233  eigen_vector_type const& _x = dynamic_cast<eigen_vector_type const&>( x );
234  return _f.vec().dot( _x.vec() );
235  }
236 
237 private:
238  solve_return_type solve( sparse_matrix_type const& A,
239  vector_type& x,
240  const vector_type& b,
241  mpl::bool_<true> );
242  solve_return_type solve( sparse_matrix_type const& A,
243  vector_type& x,
244  const vector_type& b,
245  mpl::bool_<false> );
246 
247 }; // class BackendEigen
248 
249 // -- CONSTRUCTOR --
250 template<typename T, int _Options>
251 BackendEigen<T,_Options>::BackendEigen( WorldComm const& )
252  :
253  super()
254 {}
255 
256 template<typename T, int _Options>
257 BackendEigen<T,_Options>::BackendEigen( po::variables_map const& vm, std::string const& prefix, WorldComm const& )
258  :
259  super( vm, prefix )
260 {
261  std::string _prefix = prefix;
262 
263  if ( !_prefix.empty() )
264  _prefix += "-";
265 }
266 
267 
268 
269 template<typename T, int _Options>
270 typename BackendEigen<T,_Options>::solve_return_type
271 BackendEigen<T,_Options>::solve( sparse_matrix_type const& _A,
272  vector_type& _x,
273  const vector_type& _b,
274  mpl::bool_<true>)
275 {
276  bool reusePC = ( this->precMatrixStructure() == SAME_PRECONDITIONER );
277 
278  eigen_sparse_matrix_type const& A( dynamic_cast<eigen_sparse_matrix_type const&>( _A ) );
279  eigen_vector_type & x( dynamic_cast<eigen_vector_type &>( _x ) );
280  eigen_vector_type const& b( dynamic_cast<eigen_vector_type const&>( _b ) );
281  x.vec() = A.mat().lu().solve(b.vec());
282 
283  return boost::make_tuple(true,1,1e-10);
284 } // BackendEigen::solve
285 
286 template<typename T, int _Options>
287 typename BackendEigen<T,_Options>::solve_return_type
288 BackendEigen<T,_Options>::solve( sparse_matrix_type const& _A,
289  vector_type& _x,
290  const vector_type& _b,
291  mpl::bool_<false>)
292 {
293  bool reusePC = ( this->precMatrixStructure() == SAME_PRECONDITIONER );
294 
295  eigen_sparse_matrix_type const& A( dynamic_cast<eigen_sparse_matrix_type const&>( _A ) );
296  eigen_vector_type & x( dynamic_cast<eigen_vector_type &>( _x ) );
297  eigen_vector_type const& b( dynamic_cast<eigen_vector_type const&>( _b ) );
298 
299  //x.vec()=A.mat().template fullPivLu().solve(b.vec());
300  Eigen::SimplicialLDLT<typename eigen_sparse_matrix_type::matrix_type> solver;
301  solver.compute(A.mat());
302  x.vec() = solver.solve(b.vec());
303 
304  // if(solver.info()!=Eigen::Succeeded) {
305  // // solving failed
306  // return boost::make_tuple(false,1,1e-10);;
307  // }
308  return boost::make_tuple(true,1,1e-10);;
309 } // BackendEigen::solve
310 
311 
312 
313 
314 po::options_description backendeigen_options( std::string const& prefix = "" );
315 
316 } // Feel
317 
318 #endif /* _BACKENDEIGEN_HPP_ */

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