Loading...
Searching...
No Matches
solver_trilinos_ifpack.cc
Go to the documentation of this file.
// preconditioner based on trilinos/ifpack, both seq & mpi implementations
// note: the library requires mpi to be present => no seq-only implementation
//
#include "rheolef/config.h"
#if defined(_RHEOLEF_HAVE_TRILINOS) && defined(_RHEOLEF_HAVE_MPI)
#include "solver_trilinos_ifpack.h"
#include "Ifpack.h"
#include "Epetra_Vector.h"
namespace rheolef {
using namespace std;
// convert from rheolef to trilinos
static
Epetra_CrsMatrix*
csr2petra (const csr<double,sequential>& a, Epetra_Map*& petra_ownership_ptr)
{
error_macro ("not yet");
return 0;
}
static
Epetra_CrsMatrix*
csr2petra (const csr<double,distributed>& a, Epetra_Map*& petra_ownership_ptr)
{
typedef csr<double,distributed>::size_type size_type;
distributor ownership = a.row_ownership();
Epetra_MpiComm petra_comm (ownership.comm()); // or Epetra_SerialComm
petra_ownership_ptr = new_macro (Epetra_Map (ownership.dis_size(), ownership.size(), 0, petra_comm));
std::vector<int> nz_by_row (a.nrow());
csr<double,distributed>::const_iterator dia_ia = a.begin();
csr<double,distributed>::const_iterator ext_ia = a.ext_begin();
int max_nz_by_row = 0;
int nnz = 0;
for (size_type i = 0, n = a.nrow(); i < n; i++) {
nz_by_row[i] = (dia_ia[i+1] - dia_ia[i]) + (ext_ia[i+1] - ext_ia[i]);
max_nz_by_row = std::max(max_nz_by_row, nz_by_row[i]);
nnz += nz_by_row[i];
}
bool static_storage = true;
Epetra_CrsMatrix* a_petra_ptr = new_macro (Epetra_CrsMatrix (Copy, *petra_ownership_ptr, nz_by_row.begin().operator->(), static_storage));
Epetra_CrsMatrix& a_petra = *a_petra_ptr;
std::vector<int> jdx (max_nz_by_row);
std::vector<double> val (max_nz_by_row);
size_type first_i = ownership.first_index();
for (size_type i = 0, n = a.nrow(); i < n; i++) {
size_type q = 0;
for (csr<double,distributed>::const_data_iterator p = dia_ia[i]; p < dia_ia[i+1]; p++, q++) {
jdx[q] = (*p).first + first_i;
val[q] = (*p).second;
}
for (csr<double,distributed>::const_data_iterator p = ext_ia[i]; p < ext_ia[i+1]; p++, q++) {
jdx[q] = a.jext2dis_j ((*p).first);
val[q] = (*p).second;
}
size_type dis_i = first_i + i;
check_macro (int(q) == nz_by_row[i], "unexpected");
a_petra.InsertGlobalValues(dis_i, nz_by_row[i], val.begin().operator->(), jdx.begin().operator->());
}
bool optimize = true;
a_petra.FillComplete(optimize);
return a_petra_ptr;
}
template<class T, class M>
void
solver_trilinos_ifpack_rep<T,M>::destroy_values ()
{
if (_ilu_ptr) { delete_macro(_ilu_ptr); _ilu_ptr = 0; }
if (_petra_ownership_ptr) { delete_macro(_petra_ownership_ptr); _petra_ownership_ptr = 0; }
}
template<class T, class M>
void
solver_trilinos_ifpack_rep<T,M>::update_values (const csr<T,M>& a)
{
destroy_values ();
#ifdef TO_CLEAN
// ILUT, with thresold, cannot exploit the constant sparsity pattern of a
// thus, the factorization may be recomputed completely
double k = 1; // TODO: in options
double drop_tol = 0;
int nnz = a.dis_nnz();
int n = a.dis_nrow();
int level_fill = int (k*nnz/(2.*n) + 1);
check_macro (a.is_symmetric(), "ict(k,e): unsymmetric matrix not supported");
_ilu_ptr = new_macro (Ifpack_CrsIct (*a_petra_ptr, drop_tol, level_fill));
_ilu_ptr->InitValues(*a_petra_ptr);
_ilu_ptr->Factor();
#endif // TO_CLEAN
string type = (a.is_symmetric() ? "IC" : "ILU"); // TODO: ICT & ILUT
Ifpack factory;
Epetra_CrsMatrix* a_petra_ptr = csr2petra (a, _petra_ownership_ptr);
_ilu_ptr = factory.Create (type, a_petra_ptr);
Teuchos::ParameterList params;
// TODO: fact parameters
// TODO: overlap with domains
// fact: level-of-fill for IC and ILU.
// fact: ict level-of-fill [double] Level-of-fill for ICT.
// fact: ilut level-of-fill [double] Level-of-fill for ILUT.
// fact: relax value [double] Relaxation value.
// fact: absolute threshold [double] Value of alpha in equation (16).
// fact: level-of-fill relative threshold [double] Value of rho in equation (16).
// B = alpha*sgn(A) + rho*A (16)
params.set("fact: level-of-fill", 0);
_ilu_ptr -> SetParameters (params);
_ilu_ptr -> Initialize();
_ilu_ptr -> Compute();
#ifdef TO_CLEAN
delete_macro(a_petra_ptr); // deleted by _ilu_ptr ?
#endif // TO_CLEAN
}
template<class T, class M>
solver_trilinos_ifpack_rep<T,M>::~solver_trilinos_ifpack_rep ()
{
destroy_values ();
}
template<class T, class M>
solver_trilinos_ifpack_rep<T,M>::solver_trilinos_ifpack_rep (const csr<T,M>& a, const solver_option& opt)
: solver_abstract_rep<T,M>(opt),
_ilu_ptr(0),
_petra_ownership_ptr(0)
{
// TODO: (k, drop_tol) in options
update_values (a);
}
template<class T, class M>
void
solver_trilinos_ifpack_rep<T,M>::solve (const vec<T,M>& b, bool do_transpose, vec<T,M>& x) const
{
const double* b_values = b.begin().operator->();
double* x_values = x.begin().operator->();
check_macro (int(x.ownership().size()) == _petra_ownership_ptr -> NumMyElements(),
"incomplete choleski preconditioner: incompatible right-hand size");
Epetra_Vector b_petra (View, *_petra_ownership_ptr, const_cast<double*>(b_values));
Epetra_Vector x_petra (View, *_petra_ownership_ptr, x_values);
_ilu_ptr -> SetUseTranspose (do_transpose);
_ilu_ptr -> ApplyInverse (b_petra, x_petra);
}
template<class T, class M>
vec<T,M>
solver_trilinos_ifpack_rep<T,M>::solve (const vec<T,M>& b) const
{
vec<T,M> x (b.ownership());
solve (b, false, x);
return x;
}
template<class T, class M>
vec<T,M>
solver_trilinos_ifpack_rep<T,M>::trans_solve (const vec<T,M>& b) const
{
vec<T,M> x (b.ownership());
solve (b, true, x);
return x;
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
// TODO: code is only valid here for T=double
template class solver_trilinos_ifpack_rep<double,sequential>;
#ifdef _RHEOLEF_HAVE_MPI
template class solver_trilinos_ifpack_rep<double,distributed>;
#endif // _RHEOLEF_HAVE_MPI
} // namespace rheolef
#endif // _RHEOLEF_HAVE_TRILINOS && _RHEOLEF_HAVE_MPI
rheolef::csr
see the csr page for the full documentation
Definition csr.h:317
rheolef::distributor
see the distributor page for the full documentation
Definition distributor.h:69
rheolef::solver_abstract_rep
Definition solver.h:191
rheolef::solver_option
see the solver_option page for the full documentation
Definition solver_option.h:155
error_macro
#define error_macro(message)
Definition dis_macros.h:49
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)")
mkgeo_ball.b
int b
Definition mkgeo_ball.sh:152
mkgeo_ball.n
int n
Definition mkgeo_ball.sh:150
mkgeo_ball.a
int a
Definition mkgeo_ball.sh:151
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
solver_trilinos_ifpack.h
rheolef::type
M
Expr1::memory_type M
Definition vec_expr_v2.h:416