Loading...
Searching...
No Matches
solver_eigen.cc
Go to the documentation of this file.
// direct solver eigen, seq implementations
//
// Note : why a dis implementation based on eigen ?
// Because when dis_ext_nnz == 0, then the matrix is block diagonal.
// in that case the eigen could be better than mumps that initialize stuff
// for the distributed case.
// Is could appends e.g. for block-diagonal mass matrix "Pkd"
// This also occurs when nproc==1.
//
#include "rheolef/config.h"
#ifdef _RHEOLEF_HAVE_EIGEN
#include "solver_eigen.h"
namespace rheolef {
using namespace std;
// =========================================================================
// the class interface
// =========================================================================
template<class T, class M>
void
solver_eigen_rep<T,M>::update_values (const csr<T,M>& a)
{
_a = a;
using namespace Eigen;
Matrix<int,Dynamic,1> nnz_row (a.nrow());
typename csr<T,M>::const_iterator ia = a.begin();
for (size_type i = 0, q = 0, n = a.nrow(); i < n; ++i) {
nnz_row[i] = ia[i+1] - ia[i];
}
SparseMatrix<T> a_tmp (a.nrow(),a.ncol());
if (a.nrow() != 0) {
a_tmp.reserve (nnz_row);
}
for (size_type i = 0, n = a.nrow(); i < n; ++i) {
for (typename csr<T,M>::const_data_iterator p = ia[i]; p < ia[i+1]; ++p) {
a_tmp.insert (i, (*p).first) = (*p).second;
}
}
a_tmp.makeCompressed();
if (_a.is_symmetric() && _a.is_definite_positive()) {
_ldlt_a.compute (a_tmp);
check_macro (_ldlt_a.info() == Success, "eigen ldlt factorization failed");
if (base::option().compute_determinant) {
T det_a = _ldlt_a.determinant(); // TODO: wait for eigen::ldlt to support logAbsDeterminant & signDeterminant
_det.mantissa = (det_a >= 0) ? 1 : -1;
_det.exponant = log(fabs(det_a)) / log(T(10));
_det.base = 10;
}
} else {
if (a.nrow() != 0) {
_superlu_a.analyzePattern (a_tmp);
_superlu_a.factorize (a_tmp);
check_macro (_superlu_a.info() == Success, "eigen lu factorization failed");
}
if (base::option().compute_determinant) {
_det.mantissa = _superlu_a.signDeterminant();
_det.exponant = _superlu_a.logAbsDeterminant() / log(T(10));
_det.base = 10;
}
}
}
template<class T, class M>
vec<T,M>
solver_eigen_rep<T,M>::solve (const vec<T,M>& b) const
{
if (b.dis_size() == 0) return b; // empty matrix
vec<T,M> x (b.ownership());
using namespace Eigen;
Map<Matrix<T,Dynamic,1> > b_map ((T*)(b.begin().operator->()), b.size()),
x_map ( x.begin().operator->(), x.size());
if (_a.is_symmetric() && _a.is_definite_positive()) {
x_map = _ldlt_a.solve (b_map);
} else {
x_map = _superlu_a.solve (b_map);
}
return x;
}
template<class T, class M>
vec<T,M>
solver_eigen_rep<T,M>::trans_solve (const vec<T,M>& b) const
{
if (_a.is_symmetric()) return solve(b);
if (b.dis_size() == 0) return b; // empty matrix
vec<T,M> x(b.ownership());
fatal_macro ("eigen superlu trans_solve: not yet implemented, sorry");
#ifdef TODO
x_map = _superlu_a.colssPermutation()*b_map;
_superlu_a.matrixU().transSolveInPlace (x_map);
_superlu_a.matrixL().transSolveInPlace (x_map); // L & U trans_solve not implemented
x_map = _superlu_a.rowsPermutation()*x_map;
#endif // TODO
return x;
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
template class solver_eigen_rep<Float,sequential>;
#ifdef _RHEOLEF_HAVE_MPI
template class solver_eigen_rep<Float,distributed>;
#endif // _RHEOLEF_HAVE_MPI
} // namespace rheolef
#endif // HAVE_EIGEN
rheolef::csr
see the csr page for the full documentation
Definition csr.h:317
rheolef::solver_eigen_rep
Definition solver_eigen.h:42
rheolef::solver_eigen_rep::update_values
void update_values(const csr< T, M > &a)
Definition solver_eigen.cc:42
rheolef::solver_eigen_rep::size_type
base::size_type size_type
Definition solver_eigen.h:47
rheolef::solver_eigen_rep::trans_solve
vec< T, M > trans_solve(const vec< T, M > &rhs) const
Definition solver_eigen.cc:101
rheolef::solver_eigen_rep::solve
vec< T, M > solve(const vec< T, M > &rhs) const
Definition solver_eigen.cc:85
rheolef::vec
see the vec page for the full documentation
Definition vec.h:79
fatal_macro
#define fatal_macro(message)
Definition dis_macros.h:33
T
Expr1::float_type T
Definition field_expr.h:230
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
Eigen
Definition field_eigen.h:35
rheolef
This file is part of Rheolef.
Definition compiler_eigen.h:39
rheolef::solve
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
std
STL namespace.
solver_eigen.h
p
Definition sphere.icc:25