Loading...
Searching...
No Matches
solver_pastix_mpi.cc
Go to the documentation of this file.
// direct solver based on pastix
//
#include "rheolef/config.h"
#if defined(_RHEOLEF_HAVE_PASTIX) && defined(_RHEOLEF_HAVE_MPI)
#include "solver_pastix.h"
#include "rheolef/cg.h"
#include "rheolef/eye.h"
namespace rheolef {
using namespace std;
template<class T>
solver_pastix_rep<T,distributed>::solver_pastix_rep (const csr<T>& a, const solver_option& opt)
: base(),
_comm(),
_i2dis_i_base(),
_new_n(0),
_new_ptr(0),
_new_idx(0),
_new_val(0)
{
warning_macro ("call load...");
load (a, opt);
warning_macro ("return load...");
}
template<class T>
void
solver_pastix_rep<T,distributed>::load (const csr<T>& a, const solver_option& opt)
{
warning_macro ("load...");
base::_is_sym = a.is_symmetric();
base::_pattern_dimension = a.pattern_dimension();
_comm = a.row_ownership().comm();
base::_csr_row_ownership = a.row_ownership();
base::_csr_col_ownership = a.col_ownership();
base::_opt = opt;
check_macro (a.nrow() == a.ncol(), "factorization: only square matrix are supported");
// there is a bug in pastix related to "np < a.dis_nrow" : use an alternative solver...
if (a.nrow() == 0) {
base::_have_pastix_bug_small_matrix = true;
}
// ask if other procs have a bug ?
base::_have_pastix_bug_small_matrix = mpi::all_reduce (comm(), base::_have_pastix_bug_small_matrix, mpi::maximum<bool>());
if (base::_have_pastix_bug_small_matrix) {
// one of the procs at least have a too small matrix for pastix !
base::_a_when_bug = a;
return;
}
if (base::_is_sym) {
load_symmetric(a);
} else {
load_unsymmetric(a);
}
if (base::_opt.do_check) {
check ();
}
symbolic_factorization ();
numeric_factorization();
warning_macro ("load done");
}
template<class T>
void
solver_pastix_rep<T,distributed>::load_unsymmetric (const csr<T>& a)
{
size_t n = a.nrow();
size_t nnz = a.nnz() + a.ext_nnz();
resize (n, nnz);
load_both_continued (a);
}
// compute how many a.dia (dis_i,dis_j) have dis_i > dis_j
template<class T>
static
size_t
compute_csr_upper_dia_nnz (const csr<T>& a)
{
size_t csr_upper_dia_nnz = 0;
size_t first_dis_i = a.row_ownership().first_index();
size_t first_dis_j = a.col_ownership().first_index();
for (size_t i = 0; i < a.nrow(); i++) {
size_t dis_i = first_dis_i + i;
typename csr<T>::const_iterator ia = a.begin();
for (typename csr<T>::const_data_iterator ap = ia[i]; ap < ia[i+1]; ap++) {
size_t j = (*ap).first;
size_t dis_j = first_dis_j + j;
if (dis_i <= dis_j) csr_upper_dia_nnz++;
}
}
return csr_upper_dia_nnz;
}
// compute how many a.ext (dis_i,pa_j) have dis_i > dis_j
template<class T>
static
size_t
compute_csr_upper_ext_nnz (const csr<T>& a)
{
size_t csr_upper_ext_nnz = 0;
size_t first_dis_i = a.row_ownership().first_index();
size_t first_dis_j = a.col_ownership().first_index();
for (size_t i = 0; i < a.nrow(); i++) {
size_t dis_i = first_dis_i + i;
typename csr<T>::const_iterator ext_ia = a.ext_begin();
for (typename csr<T>::const_data_iterator ap = ext_ia[i]; ap < ext_ia[i+1]; ap++) {
const size_t& jext = (*ap).first;
size_t dis_j = a.jext2dis_j (jext);
if (dis_i <= dis_j) csr_upper_ext_nnz++;
}
}
return csr_upper_ext_nnz;
}
template<class T>
void
solver_pastix_rep<T,distributed>::load_symmetric (const csr<T>& a)
{
// conserve only the lower part of the csc(pastix) = the upper past of the csr(rheolef)
size_t n = a.nrow();
size_t csr_upper_dia_nnz = compute_csr_upper_dia_nnz(a);
size_t csr_upper_ext_nnz = compute_csr_upper_ext_nnz(a);
size_t nnz = csr_upper_dia_nnz + csr_upper_ext_nnz;
resize (n, nnz);
load_both_continued (a);
}
template<class T>
void
solver_pastix_rep<T,distributed>::load_both_continued (const csr<T>& a)
{
size_t first_dis_i = a.row_ownership().first_index();
size_t first_dis_j = a.col_ownership().first_index();
typename csr<T>::const_iterator aptr = a.begin();
pastix_int_t bp = 0;
base::_ptr [0] = bp + base::_base;
for (size_t i = 0; i < a.nrow(); i++) {
size_t dis_i = first_dis_i + i;
_i2dis_i_base [i] = dis_i + base::_base;
for (typename csr<T>::const_data_iterator ap = aptr[i]; ap < aptr[i+1]; ap++) {
size_t j = (*ap).first;
const T& val = (*ap).second;
size_t dis_j = first_dis_j + j;
if (base::_is_sym && dis_i > dis_j) continue;
base::_val [bp] = val;
base::_idx [bp] = dis_j + base::_base;
bp++;
}
// note: pas tries par j croissant : dans a.ext, il i a des dis_j < first_dis_j
typename csr<T>::const_iterator ext_ia = a.ext_begin();
for (typename csr<T>::const_data_iterator ap = ext_ia[i]; ap < ext_ia[i+1]; ap++) {
size_t jext = (*ap).first;
const T& val = (*ap).second;
size_t dis_j = a.jext2dis_j (jext);
if (base::_is_sym && dis_i > dis_j) continue;
base::_val [bp] = val;
base::_idx [bp] = dis_j + base::_base;
bp++;
}
base::_ptr [i+1] = bp + base::_base;
}
check_macro (bp == base::_nnz, "factorization: invalid count: nnz="<<base::_nnz<<", count="<<bp);
}
template<class T>
void
solver_pastix_rep<T,distributed>::update_values (const csr<T>& a)
{
if (base::_have_pastix_bug_small_matrix) {
base::_a_when_bug = a;
return;
}
check_macro (size_t(base::_n) == a.nrow() && size_t(base::_n) == a.ncol(),
"update: local input matrix size distribution mismatch: ("<<a.nrow()<<","<<a.ncol()<<"), expect ("
<< base::_n << "," << base::_n << ")");
size_t nnz_a;
if (!base::_is_sym) {
nnz_a = a.nnz() + a.ext_nnz();
} else {
// conserve only the lower part of the csc(pastix) = the upper past of the csr(rheolef)
size_t csr_upper_dia_nnz = compute_csr_upper_dia_nnz(a);
size_t csr_upper_ext_nnz = compute_csr_upper_ext_nnz(a);
nnz_a = csr_upper_dia_nnz + csr_upper_ext_nnz;
}
check_macro (size_t(base::_nnz) == nnz_a,
"update: local input matrix nnz distribution mismatch: nnz="<<nnz_a<<", expect nnz="<<base::_nnz);
size_t first_dis_i = a.row_ownership().first_index();
size_t first_dis_j = a.col_ownership().first_index();
pastix_int_t bp = 0;
typename csr<T>::const_iterator aptr = a.begin();
for (size_t i = 0; i < a.nrow(); i++) {
size_t dis_i = first_dis_i + i;
for (typename csr<T>::const_data_iterator ap = aptr[i]; ap < aptr[i+1]; ap++) {
size_t j = (*ap).first;
size_t dis_j = first_dis_j + j;
if (base::_is_sym && dis_i > dis_j) continue;
base::_val [bp] = (*ap).second;
bp++;
}
typename csr<T>::const_iterator ext_ia = a.ext_begin();
for (typename csr<T>::const_data_iterator ap = ext_ia[i]; ap < ext_ia[i+1]; ap++) {
size_t jext = (*ap).first;
size_t dis_j = a.jext2dis_j (jext);
if (base::_is_sym && dis_i > dis_j) continue;
base::_val [bp] = (*ap).second;
bp++;
}
}
numeric_factorization();
}
template<class T>
void
solver_pastix_rep<T,distributed>::resize (pastix_int_t n, pastix_int_t nnz)
{
base::_n = n;
base::_nnz = nnz;
base::_ptr.resize(n+1);
base::_idx.resize(nnz);
base::_val.resize(nnz);
_i2dis_i_base.resize(n);
}
template<class T>
void
solver_pastix_rep<T,distributed>::check () const
{
pastix_int_t symmetry = (is_symmetric() ? API_SYM_YES : API_SYM_NO);
pastix_int_t *ptr_begin = (pastix_int_t*) base::_ptr.begin().operator->();
pastix_int_t *idx_begin = (pastix_int_t*) base::_idx.begin().operator->();
T *val_begin = (T*) base::_val.begin().operator->();
pastix_int_t *i2dis_i_base_begin = (pastix_int_t*) _i2dis_i_base.begin().operator->();
pastix_int_t status
= d_pastix_checkMatrix(_comm, base::_opt.verbose_level, symmetry, API_YES,
base::_n, &ptr_begin, &idx_begin, &val_begin, &i2dis_i_base_begin,
1);
check_macro (status == 0, "pastix check returns error status = " << status);
}
template<class T>
void
solver_pastix_rep<T,distributed>::symbolic_factorization ()
{
if (base::_have_pastix_bug_small_matrix) {
return;
}
const pastix_int_t nbthread = 1; // threads are not yet very well supported with openmpi & scotch
const pastix_int_t ordering = 0; // scotch
// tasks :
// 0. set params to default values
base::_iparm[IPARM_START_TASK] = API_TASK_INIT;
base::_iparm[IPARM_END_TASK] = API_TASK_INIT;
base::_iparm[IPARM_MODIFY_PARAMETER] = API_NO;
d_dpastix (base::pp_data(),
_comm,
base::_n,
base::_ptr.begin().operator->(),
base::_idx.begin().operator->(),
NULL, // _val.begin().operator->(),
_i2dis_i_base.begin().operator->(),
NULL,
NULL,
NULL,
1,
base::_iparm,
base::_dparm);
// Customize some parameters
base::_iparm[IPARM_THREAD_NBR] = nbthread;
if (is_symmetric()) {
base::_iparm[IPARM_SYM] = API_SYM_YES;
base::_iparm[IPARM_FACTORIZATION] = API_FACT_LDLT;
} else {
base::_iparm[IPARM_SYM] = API_SYM_NO;
base::_iparm[IPARM_FACTORIZATION] = API_FACT_LU;
}
base::_iparm[IPARM_MATRIX_VERIFICATION] = API_NO;
base::_iparm[IPARM_VERBOSE] = base::_opt.verbose_level;
base::_iparm[IPARM_ORDERING] = ordering;
bool do_incomplete;
if (base::_opt.iterative != solver_option::decide) {
do_incomplete = (base::_opt.iterative != 0);
} else {
do_incomplete = (base::_pattern_dimension > 2); // 3d-pattern => iterative & IC(k) precond
}
base::_iparm[IPARM_INCOMPLETE] = (do_incomplete ? 1 : 0);
base::_iparm[IPARM_OOC_LIMIT] = base::_opt.ooc;
if (base::_opt.iterative == 1) {
base::_dparm[DPARM_EPSILON_REFINEMENT] = base::_opt.tol;
}
base::_iparm[IPARM_LEVEL_OF_FILL] = base::_opt.level_of_fill;
base::_iparm[IPARM_AMALGAMATION_LEVEL] = base::_opt.amalgamation;
base::_iparm[IPARM_RHS_MAKING] = API_RHS_B;
// 1) get the i2new_dis_i array: its indexes are in the [0:dis_n[ range
// i2new_dis_i[i] = i2new_dis_i [i]
base::_i2new_dis_i.resize (base::_n);
pastix_int_t itmp = 0; // pastix does not manage NULL pointers: when n=0, vector::begin() retuns a NULL...
T vtmp = 0; // pastix does not manage NULL pointers: when n=0, vector::begin() retuns a NULL...
// tasks :
// 1. ordering
// 2. symbolic factorization
// 3. tasks mapping and scheduling
base::_iparm[IPARM_START_TASK] = API_TASK_ORDERING;
base::_iparm[IPARM_END_TASK] = API_TASK_ANALYSE;
d_dpastix (base::pp_data(),
_comm,
base::_n,
base::_ptr.begin().operator->(),
(base::_idx.begin().operator->() != NULL) ? base::_idx.begin().operator->() : &itmp,
NULL,
(_i2dis_i_base.begin().operator->() != NULL) ? _i2dis_i_base.begin().operator->() : &itmp,
(base::_i2new_dis_i.begin().operator->() != NULL) ? base::_i2new_dis_i.begin().operator->() : &itmp,
NULL,
NULL,
1,
base::_iparm,
base::_dparm);
_new_n = d_pastix_getLocalNodeNbr (base::pp_data());
// 2) buil new local to global column correspondance
// new_i2dis_i_base[new_i] - base = new_i2dis_i [new_i] with a base=1
base::_new_i2dis_i_base.resize (_new_n);
d_pastix_getLocalNodeLst (base::pp_data(), base::_new_i2dis_i_base.begin().operator->());
}
template<class T>
void
solver_pastix_rep<T,distributed>::numeric_factorization ()
{
if (base::_have_pastix_bug_small_matrix) {
trace_macro ("num_fact: bug, circumvent !");
return;
}
// 3) redistributing the matrix in the new numbering
pastix_int_t status2 = d_cscd_redispatch (
base::_n,
base::_ptr.begin().operator->(),
base::_idx.begin().operator->(),
base::_val.begin().operator->(),
NULL,
_i2dis_i_base.begin().operator->(),
_new_n,
&(_new_ptr),
&(_new_idx),
&(_new_val),
NULL,
base::_new_i2dis_i_base.begin().operator->(),
_comm);
check_macro (status2 == EXIT_SUCCESS, "d_csd_redispatch failed");
// pastix tasks:
// 4. numerical factorization
base::_iparm[IPARM_START_TASK] = API_TASK_NUMFACT;
base::_iparm[IPARM_END_TASK] = API_TASK_NUMFACT;
d_dpastix (base::pp_data(),
_comm,
_new_n,
_new_ptr,
_new_idx,
_new_val,
base::_new_i2dis_i_base.begin().operator->(),
base::_i2new_dis_i.begin().operator->(),
NULL,
NULL,
0,
base::_iparm,
base::_dparm);
}
template<class T>
vec<T>
solver_pastix_rep<T,distributed>::trans_solve (const vec<T>& rhs) const
{
if (base::_have_pastix_bug_small_matrix) {
error_macro ("trans_solve: bug, circumvent not yet !");
}
// ================================================================
// solve
// ================================================================
check_macro (rhs.size() == size_t(base::_n), "invalid rhs size="<<rhs.size()<<": expect size="<<base::_n);
// redispatch rhs separately: formally: new_rhs [new_i] := rhs [i2new_i[i]]
base::_new_rhs.resize (_new_n);
std::set<size_t> dis_i_set;
std::map<size_t,size_t> dis_i2new_i;
size_t first_dis_i = base::_csr_row_ownership.first_index();
size_t last_dis_i = base::_csr_row_ownership.last_index();
for (pastix_int_t new_i = 0; new_i < _new_n; new_i++) {
size_t dis_i = base::_new_i2dis_i_base [new_i] - base::_base;
if (dis_i >= first_dis_i && dis_i < last_dis_i) {
// value is locally available
size_t i = dis_i - first_dis_i;
base::_new_rhs [new_i] = rhs [i];
} else {
// value is owned by another processor
dis_i_set.insert (dis_i);
dis_i2new_i.insert (std::make_pair(dis_i,new_i));
}
}
// external terms: TODO: scatter begin&end in the symbolic step
std::map<size_t,T> dis_i_map;
rhs.get_dis_entry (dis_i_set, dis_i_map);
for (typename map<size_t,T>::const_iterator
iter = dis_i_map.begin(),
last = dis_i_map.end(); iter != last; iter++) {
size_t dis_i = (*iter).first;
const T& val = (*iter).second;
size_t new_i = dis_i2new_i [dis_i]; // find
base::_new_rhs [new_i] = val;
}
// tasks:
// 5. numerical solve
// 6. numerical refinement
// 7. clean
T vtmp = 0; // pastix does not manage NULL rhs pointer: when new_n=0, but vector::begin() retuns a NULL...
base::_iparm[IPARM_START_TASK] = API_TASK_SOLVE;
base::_iparm[IPARM_END_TASK] = API_TASK_REFINE;
d_dpastix (base::pp_data(),
_comm,
_new_n,
_new_ptr,
_new_idx,
_new_val,
base::_new_i2dis_i_base.begin().operator->(),
base::_i2new_dis_i.begin().operator->(),
NULL,
(_new_n != 0) ? base::_new_rhs.begin().operator->() : &vtmp,
1,
base::_iparm,
base::_dparm);
// new_i2dis_i_base [new_i] - base = new_i2dis_i [new_i]
vec<T> x (base::_csr_row_ownership);
for (pastix_int_t new_i = 0; new_i < _new_n; new_i++) {
size_t dis_i = base::_new_i2dis_i_base [new_i] - base::_base;
x.dis_entry (dis_i) = base::_new_rhs[new_i];
}
x.dis_entry_assembly();
return x;
}
template<class T>
vec<T>
solver_pastix_rep<T,distributed>::solve (const vec<T>& rhs) const
{
if (base::_have_pastix_bug_small_matrix) {
T tol = std::numeric_limits<T>::epsilon();
size_t max_iter = 1000; // TODO: use _opt to get tol & max_iter
vec<T> x (base::_a_when_bug.col_ownership(), 0);
int status = cg (base::_a_when_bug, x, rhs, eye<T>(), max_iter, tol, &derr);
check_macro (tol < std::numeric_limits<T>::epsilon(), "solver(cg): machine precision not reached: "<<tol);
return x;
}
// TODO: make a csc<T> wrapper around csr<T> and use csc<T> in form & form_assembly
// => avoid the transposition
if (! base::_is_sym) {
warning_macro ("solve: not yet supported for unsymmetric matrix");
}
return trans_solve (rhs);
}
template<class T>
solver_pastix_rep<T,distributed>::~solver_pastix_rep()
{
// tasks:
// 7. clean
base::_iparm[IPARM_START_TASK] = API_TASK_CLEAN;
base::_iparm[IPARM_END_TASK] = API_TASK_CLEAN;
d_dpastix (base::pp_data(),
_comm,
_new_n,
_new_ptr,
_new_idx,
_new_val,
NULL,
NULL,
NULL,
base::_new_rhs.begin().operator->(),
1,
base::_iparm,
base::_dparm);
// was allocated by d_cscd_redispatch()
if (_new_ptr != NULL) free (_new_ptr);
if (_new_idx != NULL) free (_new_idx);
if (_new_val != NULL) free (_new_val);
base::_pastix_data = 0;
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
// TODO: code is only valid here for T=double (d_pastix etc)
template class solver_pastix_rep<double,distributed>;
} // namespace rheolef
#endif // _RHEOLEF_HAVE_PASTIX
trace_macro
#define trace_macro(message)
Definition dis_macros.h:111
error_macro
#define error_macro(message)
Definition dis_macros.h:49
warning_macro
#define warning_macro(message)
Definition dis_macros.h:53
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.n
int n
Definition mkgeo_ball.sh:150
mkgeo_ball.a
int a
Definition mkgeo_ball.sh:151
mkgeo_contraction.status
status
Definition mkgeo_contraction.sh:290
rheolef
This file is part of Rheolef.
Definition compiler_eigen.h:39
rheolef::cg
int cg(const Matrix &A, Vector &x, const Vector2 &Mb, const Preconditioner &M, const solver_option &sopt=solver_option())
Definition cg.h:94
std
STL namespace.
check
void check(Float p, const field &uh)
Definition p_laplacian_post.cc:49
load
void load(idiststream &in, Float &p, field &uh)
Definition p_laplacian_post.cc:64
solver_pastix.h