31 #ifndef _BACKENDEIGEN_HPP_
32 #define _BACKENDEIGEN_HPP_
34 #include <boost/program_options/variables_map.hpp>
39 #include <feel/feelalg/matrixeigendense.hpp>
40 #include <feel/feelalg/matrixeigensparse.hpp>
45 #include <Eigen/Dense>
46 #include <Eigen/Sparse>
51 namespace po = boost::program_options;
59 template<
typename T,
int _Options = 0>
66 typedef typename super::value_type value_type;
67 typedef typename super::real_type real_type;
69 static const bool IsDense = (_Options == 1);
70 static const bool IsSparse = (_Options == 0);
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;
81 typedef typename sparse_matrix_type::graph_ptrtype graph_ptrtype;
85 typedef typename super::vector_ptrtype vector_ptrtype;
87 typedef boost::shared_ptr<vector_type> eigen_vector_ptrtype;
89 typedef typename super::solve_return_type solve_return_type;
90 typedef typename super::nl_solve_return_type nl_solve_return_type;
93 typedef typename super::datamap_ptrtype datamap_ptrtype;
97 BackendEigen( WorldComm
const& worldComm=Environment::worldComm() );
99 WorldComm
const& worldComm=Environment::worldComm() );
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,
107 Context ctx( matrix_properties );
110 auto A= sparse_matrix_ptrtype(
new eigen_sparse_matrix_type( space1->nDof(), space2->nDof() ) );
111 A->setMatrixProperties( matrix_properties );
116 sparse_matrix_ptrtype
125 sparse_matrix_ptrtype mat(
new eigen_sparse_matrix_type( m,n ) );
126 mat->setMatrixProperties( matrix_properties );
130 sparse_matrix_ptrtype
135 graph_ptrtype
const & graph,
138 sparse_matrix_ptrtype mat(
new eigen_sparse_matrix_type( m,n ) );
139 mat->setMatrixProperties( matrix_properties );
143 sparse_matrix_ptrtype
146 auto A = sparse_matrix_ptrtype(
new eigen_sparse_matrix_type( d1->nGlobalElements(), d2->nGlobalElements() ) );
147 A->setMatrixProperties( matrix_properties );
151 sparse_matrix_ptrtype
157 auto A = sparse_matrix_ptrtype(
new eigen_sparse_matrix_type( m, n ) );
163 sparse_matrix_ptrtype
164 newZeroMatrix( datamap_ptrtype
const& d1, datamap_ptrtype
const& d2 )
166 auto A = sparse_matrix_ptrtype(
new eigen_sparse_matrix_type( d1->nGlobalElements(), d2->nGlobalElements() ) );
171 template<
typename SpaceT>
172 static vector_ptrtype newVector( boost::shared_ptr<SpaceT>
const& space )
174 return vector_ptrtype(
new eigen_vector_type( space->nDof() ) );
177 template<
typename SpaceT>
178 static vector_ptrtype newVector( SpaceT
const& space )
180 return vector_ptrtype(
new eigen_vector_type( space.nDof() ) );
200 eigen_sparse_matrix_type
const& _A =
dynamic_cast<eigen_sparse_matrix_type const&
>( A );
203 _b.
vec() = _A.mat()*_x.
vec();
206 solve_return_type solve( sparse_matrix_type
const& A,
210 return this->solve( A, x, b, mpl::bool_<IsDense>() );
213 solve_return_type solve( sparse_matrix_ptrtype
const& A,
215 const vector_ptrtype& b )
217 return this->solve( *A, *x, *b );
220 solve_return_type
solve( sparse_matrix_ptrtype
const& A,
221 sparse_matrix_ptrtype
const& P,
223 const vector_ptrtype& b )
225 return this->solve( *A, *x, *b );
234 return _f.
vec().dot( _x.
vec() );
238 solve_return_type solve( sparse_matrix_type
const& A,
242 solve_return_type solve( sparse_matrix_type
const& A,
250 template<
typename T,
int _Options>
251 BackendEigen<T,_Options>::BackendEigen( WorldComm
const& )
256 template<
typename T,
int _Options>
257 BackendEigen<T,_Options>::BackendEigen( po::variables_map
const& vm, std::string
const& prefix, WorldComm
const& )
261 std::string _prefix =
prefix;
263 if ( !_prefix.empty() )
269 template<
typename T,
int _Options>
270 typename BackendEigen<T,_Options>::solve_return_type
271 BackendEigen<T,_Options>::solve( sparse_matrix_type
const& _A,
273 const vector_type& _b,
276 bool reusePC = ( this->precMatrixStructure() == SAME_PRECONDITIONER );
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());
283 return boost::make_tuple(
true,1,1e-10);
286 template<
typename T,
int _Options>
287 typename BackendEigen<T,_Options>::solve_return_type
288 BackendEigen<T,_Options>::solve( sparse_matrix_type
const& _A,
290 const vector_type& _b,
293 bool reusePC = ( this->precMatrixStructure() == SAME_PRECONDITIONER );
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 ) );
300 Eigen::SimplicialLDLT<typename eigen_sparse_matrix_type::matrix_type> solver;
301 solver.compute(A.mat());
302 x.vec() = solver.solve(b.vec());
308 return boost::make_tuple(
true,1,1e-10);;