45 using namespace Eigen;
46 Matrix<int,Dynamic,1> nnz_row (a.nrow());
48 for (
size_type i = 0, q = 0, n = a.nrow(); i < n; ++i) {
49 nnz_row[i] = ia[i+1] - ia[i];
51 SparseMatrix<T> a_tmp (a.nrow(),a.ncol());
53 a_tmp.reserve (nnz_row);
55 for (
size_type i = 0, n = a.nrow(); i < n; ++i) {
57 a_tmp.insert (i, (*p).first) = (*p).second;
60 a_tmp.makeCompressed();
61 if (_a.is_symmetric() && _a.is_definite_positive()) {
62 _ldlt_a.compute (a_tmp);
63 check_macro (_ldlt_a.info() == Success,
"eigen ldlt factorization failed");
64 if (base::option().compute_determinant) {
65 T det_a = _ldlt_a.determinant();
66 _det.mantissa = (det_a >= 0) ? 1 : -1;
67 _det.exponant = log(fabs(det_a)) / log(
T(10));
72 _superlu_a.analyzePattern (a_tmp);
73 _superlu_a.factorize (a_tmp);
74 check_macro (_superlu_a.info() == Success,
"eigen lu factorization failed");
76 if (base::option().compute_determinant) {
77 _det.mantissa = _superlu_a.signDeterminant();
78 _det.exponant = _superlu_a.logAbsDeterminant() / log(
T(10));
87 if (b.dis_size() == 0)
return b;
89 using namespace Eigen;
90 Map<Matrix<T,Dynamic,1> > b_map ((
T*)(b.begin().operator->()), b.size()),
91 x_map ( x.begin().operator->(), x.size());
92 if (_a.is_symmetric() && _a.is_definite_positive()) {
93 x_map = _ldlt_a.solve (b_map);
95 x_map = _superlu_a.solve (b_map);
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")