Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
solver_cholmod.cc
Go to the documentation of this file.
1
21// 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)
30namespace rheolef {
31using namespace std;
32
33template<class T, class M>
34void
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}
67template<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}
79template<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
see the csr page for the full documentation
Definition csr.h:317
void update_values(const csr< T, M > &a)
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
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