Loading...
Searching...
No Matches
mic.cc
Go to the documentation of this file.
// preconditioner eigen, seq implementations
//
#include "rheolef/mic.h"
#ifdef _RHEOLEF_HAVE_EIGEN
namespace rheolef {
using namespace std;
// =========================================================================
// the class interface
// =========================================================================
template<class T, class M>
void
solver_mic_rep<T,M>::update_values (const csr<T,M>& 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.nrow() != 0) {
_mic_a.setInitialShift (_shift);
_mic_a.compute (a_tmp);
check_macro (_mic_a.info() == Success, "eigen MIC incomplete factorization failed");
}
}
template<class T, class M>
vec<T,M>
solver_mic_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())), b.size()),
x_map ( &(*x.begin()), x.size());
x_map = _mic_a.solve (b_map);
return x;
}
template<class T, class M>
vec<T,M>
solver_mic_rep<T,M>::trans_solve (const vec<T,M>& b) const
{
return solve(b);
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
template class solver_mic_rep<Float,sequential>;
#ifdef _RHEOLEF_HAVE_MPI
template class solver_mic_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_mic_rep
Definition mic.h:93
rheolef::solver_mic_rep::update_values
void update_values(const csr< T, M > &a)
Definition mic.cc:34
rheolef::solver_mic_rep::size_type
base::size_type size_type
Definition mic.h:98
rheolef::solver_mic_rep::trans_solve
vec< T, M > trans_solve(const vec< T, M > &rhs) const
Definition mic.cc:72
rheolef::solver_mic_rep::solve
vec< T, M > solve(const vec< T, M > &rhs) const
Definition mic.cc:60
rheolef::solver_mic_rep::_shift
T _shift
Definition mic.h:117
rheolef::solver_mic_rep::_mic_a
Eigen::IncompleteCholesky< T, Eigen::Lower, Eigen::AMDOrdering< int > > _mic_a
Definition mic.h:118
rheolef::vec
see the vec page for the full documentation
Definition vec.h:79
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
p
Definition sphere.icc:25