23#include "rheolef/ilut.h"
25#ifdef _RHEOLEF_HAVE_EIGEN
32template<
class T,
class M>
36 using namespace Eigen;
37 Matrix<int,Dynamic,1> nnz_row (a.nrow());
39 for (
size_type i = 0, q = 0, n = a.nrow(); i < n; ++i) {
40 nnz_row[i] = ia[i+1] - ia[i];
42 SparseMatrix<T> a_tmp (a.nrow(),a.ncol());
44 a_tmp.reserve (nnz_row);
46 for (
size_type i = 0, n = a.nrow(); i < n; ++i) {
48 a_tmp.insert (i, (*p).first) = (*p).second;
51 a_tmp.makeCompressed();
53 _ilut_a.setFillfactor (
int(_fill_factor));
54 _ilut_a.setDroptol (
T(_drop_tol));
55 _ilut_a.compute (a_tmp);
56 check_macro (_ilut_a.info() == Success,
"building ILUT preconditioner failed");
59template<
class T,
class M>
63 if (b.dis_size() == 0)
return b;
65 using namespace Eigen;
66 Map<Matrix<T,Dynamic,1> > b_map ((
T*)(&(*b.begin())), b.size()),
67 x_map ( &(*x.begin()), x.size());
68 x_map = _ilut_a.solve (b_map);
71template<
class T,
class M>
75 if (b.dis_size() == 0)
return b;
77 fatal_macro (
"ilut trans_solve: not implemented, sorry");
86#ifdef _RHEOLEF_HAVE_MPI
see the csr page for the full documentation
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
#define fatal_macro(message)
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.