Rheolef  7.1
an efficient C++ finite element environment
solver_eigen.cc
Go to the documentation of this file.
1 // 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 
34 namespace rheolef {
35 using namespace std;
36 
37 // =========================================================================
38 // the class interface
39 // =========================================================================
40 template<class T, class M>
41 void
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 }
83 template<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 }
99 template<class T, class M>
100 vec<T,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
Eigen
Definition: field_eigen.h:35
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
solver_eigen.h
rheolef::solver_eigen_rep::solve
vec< T, M > solve(const vec< T, M > &rhs) const
Definition: solver_eigen.cc:85
rheolef::vec
see the vec page for the full documentation
Definition: vec.h:79
p
Definition: sphere.icc:25
a
Definition: diffusion_isotropic.h:25
rheolef::solver_eigen_rep::trans_solve
vec< T, M > trans_solve(const vec< T, M > &rhs) const
Definition: solver_eigen.cc:101
rheolef::csr
see the csr page for the full documentation
Definition: asr.h:31
fatal_macro
#define fatal_macro(message)
Definition: dis_macros.h:33
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
rheolef::solver_abstract_rep::size_type
csr< T, M >::size_type size_type
Definition: solver.h:193
rheolef::solve
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
mkgeo_ball.b
b
Definition: mkgeo_ball.sh:152
mkgeo_ball.a
a
Definition: mkgeo_ball.sh:151
rheolef::solver_eigen_rep::update_values
void update_values(const csr< T, M > &a)
Definition: solver_eigen.cc:42
mkgeo_ball.n
n
Definition: mkgeo_ball.sh:150
rheolef::solver_eigen_rep
Definition: solver_eigen.h:42
T
Expr1::float_type T
Definition: field_expr.h:218