Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
ilut.cc
Go to the documentation of this file.
1
21// ilut, seq implementations
22//
23#include "rheolef/ilut.h"
24
25#ifdef _RHEOLEF_HAVE_EIGEN
26namespace rheolef {
27using namespace std;
28
29// =========================================================================
30// the class interface
31// =========================================================================
32template<class T, class M>
33void
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 _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");
57 }
58}
59template<class T, class M>
62{
63 if (b.dis_size() == 0) return b; // empty matrix
64 vec<T,M> x(b.ownership());
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);
69 return x;
70}
71template<class T, class M>
74{
75 if (b.dis_size() == 0) return b; // empty matrix
76 vec<T,M> x(b.ownership());
77 fatal_macro ("ilut trans_solve: not implemented, sorry"); // TODO
78 return x;
79}
80// ----------------------------------------------------------------------------
81// instanciation in library
82// ----------------------------------------------------------------------------
83
85
86#ifdef _RHEOLEF_HAVE_MPI
88#endif // _RHEOLEF_HAVE_MPI
89
90} // namespace rheolef
91#endif // HAVE_EIGEN
see the csr page for the full documentation
Definition csr.h:317
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:73
vec< T, M > solve(const vec< T, M > &rhs) const
Definition ilut.cc:61
see the vec page for the full documentation
Definition vec.h:79
#define fatal_macro(message)
Definition dis_macros.h:33
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.
STL namespace.
Definition sphere.icc:25