Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
solver_eigen.cc
Go to the documentation of this file.
1
21// direct solver eigen, seq implementations
22//
23// Note : why a dis implementation based on eigen ?
24// Because when dis_ext_nnz == 0, then the matrix is block diagonal.
25// in that case the eigen could be better than mumps that initialize stuff
26// for the distributed case.
27// Is could appends e.g. for block-diagonal mass matrix "Pkd"
28// This also occurs when nproc==1.
29//
30#include "rheolef/config.h"
31#ifdef _RHEOLEF_HAVE_EIGEN
32#include "solver_eigen.h"
33
34namespace rheolef {
35using namespace std;
36
37// =========================================================================
38// the class interface
39// =========================================================================
40template<class T, class M>
41void
43{
44 _a = a;
45 using namespace Eigen;
46 Matrix<int,Dynamic,1> nnz_row (a.nrow());
47 typename csr<T,M>::const_iterator ia = a.begin();
48 for (size_type i = 0, q = 0, n = a.nrow(); i < n; ++i) {
49 nnz_row[i] = ia[i+1] - ia[i];
50 }
51 SparseMatrix<T> a_tmp (a.nrow(),a.ncol());
52 if (a.nrow() != 0) {
53 a_tmp.reserve (nnz_row);
54 }
55 for (size_type i = 0, n = a.nrow(); i < n; ++i) {
56 for (typename csr<T,M>::const_data_iterator p = ia[i]; p < ia[i+1]; ++p) {
57 a_tmp.insert (i, (*p).first) = (*p).second;
58 }
59 }
60 a_tmp.makeCompressed();
61 if (_a.is_symmetric() && _a.is_definite_positive()) {
62 _ldlt_a.compute (a_tmp);
63 check_macro (_ldlt_a.info() == Success, "eigen ldlt factorization failed");
64 if (base::option().compute_determinant) {
65 T det_a = _ldlt_a.determinant(); // TODO: wait for eigen::ldlt to support logAbsDeterminant & signDeterminant
66 _det.mantissa = (det_a >= 0) ? 1 : -1;
67 _det.exponant = log(fabs(det_a)) / log(T(10));
68 _det.base = 10;
69 }
70 } else {
71 if (a.nrow() != 0) {
72 _superlu_a.analyzePattern (a_tmp);
73 _superlu_a.factorize (a_tmp);
74 check_macro (_superlu_a.info() == Success, "eigen lu factorization failed");
75 }
76 if (base::option().compute_determinant) {
77 _det.mantissa = _superlu_a.signDeterminant();
78 _det.exponant = _superlu_a.logAbsDeterminant() / log(T(10));
79 _det.base = 10;
80 }
81 }
82}
83template<class T, class M>
86{
87 if (b.dis_size() == 0) return b; // empty matrix
88 vec<T,M> x (b.ownership());
89 using namespace Eigen;
90 Map<Matrix<T,Dynamic,1> > b_map ((T*)(b.begin().operator->()), b.size()),
91 x_map ( x.begin().operator->(), x.size());
92 if (_a.is_symmetric() && _a.is_definite_positive()) {
93 x_map = _ldlt_a.solve (b_map);
94 } else {
95 x_map = _superlu_a.solve (b_map);
96 }
97 return x;
98}
99template<class T, class M>
102{
103 if (_a.is_symmetric()) return solve(b);
104 if (b.dis_size() == 0) return b; // empty matrix
105 vec<T,M> x(b.ownership());
106 fatal_macro ("eigen superlu trans_solve: not yet implemented, sorry");
107#ifdef TODO
108 x_map = _superlu_a.colssPermutation()*b_map;
109 _superlu_a.matrixU().transSolveInPlace (x_map);
110 _superlu_a.matrixL().transSolveInPlace (x_map); // L & U trans_solve not implemented
111 x_map = _superlu_a.rowsPermutation()*x_map;
112#endif // TODO
113 return x;
114}
115// ----------------------------------------------------------------------------
116// instanciation in library
117// ----------------------------------------------------------------------------
118
120
121#ifdef _RHEOLEF_HAVE_MPI
123#endif // _RHEOLEF_HAVE_MPI
124
125} // namespace rheolef
126#endif // HAVE_EIGEN
see the csr page for the full documentation
Definition csr.h:317
void update_values(const csr< T, M > &a)
base::size_type size_type
vec< T, M > trans_solve(const vec< T, M > &rhs) const
vec< T, M > solve(const vec< T, M > &rhs) const
see the vec page for the full documentation
Definition vec.h:79
#define fatal_macro(message)
Definition dis_macros.h:33
Expr1::float_type T
Definition field_expr.h:230
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
This file is part of Rheolef.
void solve(tiny_matrix< T > &a, tiny_vector< size_t > &piv, const tiny_vector< T > &b, tiny_vector< T > &x)
Definition tiny_lu.h:92
STL namespace.
Definition sphere.icc:25