Rheolef  7.1
an efficient C++ finite element environment
ilut.cc
Go to the documentation of this file.
1 // ilut, seq implementations
22 //
23 #include "rheolef/ilut.h"
24 
25 #ifdef _RHEOLEF_HAVE_EIGEN
26 namespace rheolef {
27 using namespace std;
28 
29 // =========================================================================
30 // the class interface
31 // =========================================================================
32 template<class T, class M>
33 void
35 {
36  using namespace Eigen;
37  Matrix<int,Dynamic,1> nnz_row (a.nrow());
38  typename csr<T,M>::const_iterator ia = a.begin();
39  for (size_type i = 0, q = 0, n = a.nrow(); i < n; ++i) {
40  nnz_row[i] = ia[i+1] - ia[i];
41  }
42  SparseMatrix<T> a_tmp (a.nrow(),a.ncol());
43  if (a.nrow() != 0) {
44  a_tmp.reserve (nnz_row);
45  }
46  for (size_type i = 0, n = a.nrow(); i < n; ++i) {
47  for (typename csr<T,M>::const_data_iterator p = ia[i]; p < ia[i+1]; ++p) {
48  a_tmp.insert (i, (*p).first) = (*p).second;
49  }
50  }
51  a_tmp.makeCompressed();
52  if (a.nrow() != 0) {
53 warning_macro ("fill_factor="<<_fill_factor<<" drop_tol="<<_drop_tol);
54  _ilut_a.setFillfactor (int(_fill_factor));
55  _ilut_a.setDroptol (T(_drop_tol));
56  _ilut_a.compute (a_tmp);
57  check_macro (_ilut_a.info() == Success, "building ILUT preconditioner failed");
58  }
59 }
60 template<class T, class M>
63 {
64  if (b.dis_size() == 0) return b; // empty matrix
65  vec<T,M> x(b.ownership());
66  using namespace Eigen;
67  Map<Matrix<T,Dynamic,1> > b_map ((T*)(&(*b.begin())), b.size()),
68  x_map ( &(*x.begin()), x.size());
69  x_map = _ilut_a.solve (b_map);
70  return x;
71 }
72 template<class T, class M>
75 {
76  if (b.dis_size() == 0) return b; // empty matrix
77  vec<T,M> x(b.ownership());
78  fatal_macro ("ilut trans_solve: not implemented, sorry"); // TODO
79  return x;
80 }
81 // ----------------------------------------------------------------------------
82 // instanciation in library
83 // ----------------------------------------------------------------------------
84 
86 
87 #ifdef _RHEOLEF_HAVE_MPI
89 #endif // _RHEOLEF_HAVE_MPI
90 
91 } // namespace rheolef
92 #endif // HAVE_EIGEN
void update_values(const csr< T, M > &a)
Definition: ilut.cc:34
base::size_type size_type
Definition: ilut.h:103
vec< T, M > trans_solve(const vec< T, M > &rhs) const
Definition: ilut.cc:74
vec< T, M > solve(const vec< T, M > &rhs) const
Definition: ilut.cc:62
Expr1::float_type T
Definition: field_expr.h:261
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.
Definition: sphere.icc:25