Rheolef  7.1
an efficient C++ finite element environment
solver_cholmod.cc
Go to the documentation of this file.
1 // direct solver cholmod, seq implementations
22 // use the eigen interface, for convenience
23 //
24 // - cholmod LLt supernodal : when matrix is symmetric definite positive
25 // TODO:
26 // - cholmod simplicial could be considered when matrix is sym. indefinite
27 //
28 #include "solver_cholmod.h"
29 #if defined(_RHEOLEF_HAVE_CHOLMOD)
30 namespace rheolef {
31 using namespace std;
32 
33 template<class T, class M>
34 void
36 {
37  // copy into eigen; then eigen->cholmod without copy: use pointers to arrays (ia,ja,a)
38  _a = a;
39  if (_a.nnz() == 0) {
40  return; // empty matrix
41  }
42  using namespace Eigen;
43  Matrix<int,Dynamic,1> nnz_row (a.nrow());
44  typename csr<T,M>::const_iterator ia = a.begin();
45  for (size_type i = 0, q = 0, n = a.nrow(); i < n; ++i) {
46  nnz_row[i] = ia[i+1] - ia[i];
47  }
48  SparseMatrix<T> a_tmp (a.nrow(),a.ncol());
49  a_tmp.reserve (nnz_row);
50  for (size_type i = 0, n = a.nrow(); i < n; ++i) {
51  for (typename csr<T,M>::const_data_iterator p = ia[i]; p < ia[i+1]; ++p) {
52  a_tmp.insert (i, (*p).first) = (*p).second;
53  }
54  }
55  a_tmp.makeCompressed();
56  _llt_a.compute (a_tmp);
57  check_macro (_llt_a.info() == Success, "cholmod LLt factorization failed: non-positive definite matrix");
58  // TODO: note: when matrix is distributed and block diagonal, each proc returns its block determinant
59  // => do the product
60  if (base::option().compute_determinant) {
61  T det_a = _llt_a.determinant();
62  _det.mantissa = 1; // sign, always positive
63  _det.exponant = _llt_a.logDeterminant() / log(T(10));
64  _det.base = 10;
65  }
66 }
67 template<class T, class M>
70 {
71  vec<T,M> x(b.ownership());
72  if (_a.nnz() == 0) return x; // empty matrix
73  using namespace Eigen;
74  Map<Matrix<T,Dynamic,1> > b_map ((T*)(&(*b.begin())), b.size()),
75  x_map ( &(*x.begin()), x.size());
76  x_map = _llt_a.solve (b_map);
77  return x;
78 }
79 template<class T, class M>
82 {
83  return solve(b);
84 }
85 // ----------------------------------------------------------------------------
86 // instanciation in library
87 // ----------------------------------------------------------------------------
88 
90 
91 #ifdef _RHEOLEF_HAVE_MPI
93 #endif // _RHEOLEF_HAVE_MPI
94 
95 } // namespace rheolef
96 #endif // _RHEOLEF_HAVE_CHOLMOD
mkgeo_ball.n
int n
Definition: mkgeo_ball.sh:150
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)")
mkgeo_ball.b
int b
Definition: mkgeo_ball.sh:152
rheolef::vec
see the vec page for the full documentation
Definition: vec.h:79
p
Definition: sphere.icc:25
rheolef::solver_cholmod_rep::solve
vec< T, M > solve(const vec< T, M > &rhs) const
Definition: solver_cholmod.cc:69
a
Definition: diffusion_isotropic.h:25
rheolef::csr
see the csr page for the full documentation
Definition: csr.h:317
rheolef::solver_cholmod_rep::trans_solve
vec< T, M > trans_solve(const vec< T, M > &rhs) const
Definition: solver_cholmod.cc:81
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
solver_cholmod.h
rheolef::solver_cholmod_rep::update_values
void update_values(const csr< T, M > &a)
Definition: solver_cholmod.cc:35
mkgeo_ball.a
int a
Definition: mkgeo_ball.sh:151
rheolef::solver_cholmod_rep
Definition: solver_cholmod.h:40
T
Expr1::float_type T
Definition: field_expr.h:261