Loading...
Searching...
No Matches
solver_pastix_seq.cc
Go to the documentation of this file.
// c++ & boost::mpi adaptation of pastix "simple.c" sequential example
//
#include "rheolef/config.h"
#ifdef _RHEOLEF_HAVE_PASTIX
#include "solver_pastix.h"
#include "rheolef/cg.h"
#include "rheolef/eye.h"
namespace rheolef {
using namespace std;
template<class T, class M>
void
solver_pastix_base_rep<T,M>::resize (pastix_int_t n, pastix_int_t nnz)
{
_n = n;
_nnz = nnz;
_ptr.resize(n+1);
_idx.resize(nnz);
_val.resize(nnz);
}
template<class T, class M>
void
solver_pastix_base_rep<T,M>::check () const
{
trace_macro ("check...");
if (_n == 0) return;
pastix_int_t symmetry = (is_symmetric() ? API_SYM_YES : API_SYM_NO);
pastix_int_t *ptr_begin = (pastix_int_t*) _ptr.begin().operator->();
pastix_int_t *idx_begin = (pastix_int_t*) _idx.begin().operator->();
T *val_begin = (T*) _val.begin().operator->();
pastix_int_t status = d_pastix_checkMatrix (0, _opt.verbose_level, symmetry, API_YES,
_n, &ptr_begin, &idx_begin, &val_begin, NULL, 1);
check_macro (status == 0, "pastix check returns error status = " << status);
trace_macro ("pastix matrix check: ok");
}
template<class T, class M>
void
solver_pastix_base_rep<T,M>::load_both_continued (const csr<T,M>& a)
{
warning_macro ("load both...");
size_t first_dis_i = a.row_ownership().first_index();
size_t first_dis_j = a.col_ownership().first_index();
typename csr<T,M>::const_iterator aptr = a.begin();
pastix_int_t bp = 0;
_ptr [0] = bp + _base;
for (size_t i = 0; i < a.nrow(); i++) {
size_t dis_i = first_dis_i + i;
for (typename csr<T,M>::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 (_is_sym && dis_i > dis_j) continue;
_val [bp] = val;
_idx [bp] = dis_j + _base;
bp++;
}
_ptr [i+1] = bp + _base;
}
check_macro (bp == _nnz, "factorization: invalid nnz count");
warning_macro ("load both done");
}
template<class T, class M>
void
solver_pastix_base_rep<T,M>::load_unsymmetric (const csr<T,M>& a)
{
size_t n = a.nrow();
size_t nnz = a.nnz();
resize (n, nnz);
load_both_continued (a);
}
template<class T, class M>
void
solver_pastix_base_rep<T,M>::load_symmetric (const csr<T,M>& a)
{
warning_macro ("load sym...");
// count nz entries in the upper+dia part
size_t nnz_upper_dia = 0;
size_t first_dis_i = a.row_ownership().first_index();
size_t first_dis_j = a.col_ownership().first_index();
typename csr<T,M>::const_iterator aptr = a.begin();
for (size_t i = 0, n = a.nrow(); i < n; i++) {
size_t dis_i = first_dis_i + i;
for (typename csr<T,M>::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 (dis_i <= dis_j) nnz_upper_dia++;
}
}
warning_macro ("load sym(1): nrow="<<a.nrow()<<"...");
warning_macro ("load sym(1): nnz="<<a.nnz()<<"...");
warning_macro ("load sym(1): nnz_upper_dia="<<nnz_upper_dia<<"...");
resize (a.nrow(), nnz_upper_dia);
warning_macro ("load sym(2)...");
load_both_continued (a);
warning_macro ("load sym done");
}
template<class T, class M>
void
solver_pastix_base_rep<T,M>::symbolic_factorization ()
{
if (_n == 0) return;
if (_have_pastix_bug_small_matrix) {
trace_macro ("sym_fact: bug, circumvent !");
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
_iparm[IPARM_START_TASK] = API_TASK_INIT;
_iparm[IPARM_END_TASK] = API_TASK_INIT;
_iparm[IPARM_MODIFY_PARAMETER] = API_NO;
d_pastix (&_pastix_data,
0,
_n,
_ptr.begin().operator->(),
_idx.begin().operator->(),
NULL, // _val.begin().operator->(),
NULL,
NULL,
NULL,
0,
_iparm,
_dparm);
// Customize some parameters
_iparm[IPARM_THREAD_NBR] = nbthread;
if (is_symmetric()) {
_iparm[IPARM_SYM] = API_SYM_YES;
_iparm[IPARM_FACTORIZATION] = API_FACT_LDLT;
} else {
_iparm[IPARM_SYM] = API_SYM_NO;
_iparm[IPARM_FACTORIZATION] = API_FACT_LU;
}
_iparm[IPARM_MATRIX_VERIFICATION] = API_NO;
_iparm[IPARM_VERBOSE] = _opt.verbose_level;
_iparm[IPARM_ORDERING] = ordering;
bool do_incomplete;
if (_opt.iterative != solver_option::decide) {
do_incomplete = (_opt.iterative != 0);
} else {
do_incomplete = (_pattern_dimension > 2); // 3d-pattern => iterative & IC(k) precond
}
_iparm[IPARM_INCOMPLETE] = (do_incomplete ? 1 : 0);
_iparm[IPARM_OOC_LIMIT] = _opt.ooc;
if (_opt.iterative == 1) {
_dparm[DPARM_EPSILON_REFINEMENT] = _opt.tol;
}
_iparm[IPARM_LEVEL_OF_FILL] = _opt.level_of_fill;
_iparm[IPARM_AMALGAMATION_LEVEL] = _opt.amalgamation;
_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]
_i2new_dis_i.resize (_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
_iparm[IPARM_START_TASK] = API_TASK_ORDERING;
_iparm[IPARM_END_TASK] = API_TASK_ANALYSE;
std::vector<pastix_int_t> new_i2i (_n);
std::vector<double> dummy_rhs (_n);
d_pastix (&_pastix_data,
0,
_n,
_ptr.begin().operator->(),
(_idx.begin().operator->() != NULL) ? _idx.begin().operator->() : &itmp,
0,
(_i2new_dis_i.begin().operator->() != NULL) ? _i2new_dis_i.begin().operator->() : &itmp,
(new_i2i.begin().operator->() != NULL) ? new_i2i.begin().operator->() : &itmp,
0,
0,
_iparm,
_dparm);
}
template<class T, class M>
void
solver_pastix_base_rep<T,M>::numeric_factorization ()
{
if (_n == 0) return;
if (_have_pastix_bug_small_matrix) {
trace_macro ("num_fact: bug, circumvent !");
return;
}
// pastix tasks:
// 4. numerical factorization
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 rhs pointer: when new_n=0, but vector::begin() retuns a NULL...
_iparm[IPARM_START_TASK] = API_TASK_NUMFACT;
_iparm[IPARM_END_TASK] = API_TASK_NUMFACT;
d_pastix (&_pastix_data,
0,
_n,
_ptr.begin().operator->(),
(_idx.begin().operator->() != NULL) ? _idx.begin().operator->() : &itmp,
(_val.begin().operator->() != NULL) ? _val.begin().operator->() : &vtmp,
(_i2new_dis_i.begin().operator->() != NULL) ? _i2new_dis_i.begin().operator->() : &itmp,
NULL,
NULL,
0,
_iparm,
_dparm);
}
template<class T, class M>
void
solver_pastix_base_rep<T,M>::load (const csr<T,M>& a, const solver_option& opt)
{
warning_macro ("load...");
#define _RHEOLEF_SEQ_PASTIX_HAVE_SYM_BUG
#ifdef _RHEOLEF_SEQ_PASTIX_HAVE_SYM_BUG
_is_sym = false;
#else // ! _RHEOLEF_SEQ_PASTIX_HAVE_SYM_BUG
_is_sym = a.is_symmetric();
#endif // ! _RHEOLEF_SEQ_PASTIX_HAVE_SYM_BUG
_pattern_dimension = a.pattern_dimension();
_csr_row_ownership = a.row_ownership();
_csr_col_ownership = a.col_ownership();
_opt = opt;
check_macro (a.nrow() == a.ncol(), "factorization: only square matrix are supported");
if (_is_sym) {
load_symmetric(a);
} else {
load_unsymmetric(a);
}
if (_opt.do_check) {
check ();
}
warning_macro ("symbfact...");
symbolic_factorization ();
warning_macro ("numbfact...");
numeric_factorization();
warning_macro ("load done");
}
template<class T, class M>
solver_pastix_base_rep<T,M>::solver_pastix_base_rep ()
: solver_abstract_rep<T,M>(solver_option()),
_n(),
_nnz(),
_ptr(),
_idx(),
_val(),
_is_sym(false),
_pattern_dimension(0),
_pastix_data(0),
_iparm(),
_dparm(),
_csr_row_ownership(),
_csr_col_ownership(),
_opt(),
_new_rhs(),
_new_i2dis_i_base(),
_i2new_dis_i(),
_have_pastix_bug_small_matrix(false),
_a_when_bug()
{
}
template<class T, class M>
solver_pastix_base_rep<T,M>::solver_pastix_base_rep (const csr<T,M>& a, const solver_option& opt)
: solver_abstract_rep<T,M>(opt),
_n(),
_nnz(),
_ptr(),
_idx(),
_val(),
_is_sym(false),
_pattern_dimension(0),
_pastix_data(0),
_iparm(),
_dparm(),
_csr_row_ownership(),
_csr_col_ownership(),
_opt(),
_new_rhs(),
_new_i2dis_i_base(),
_i2new_dis_i(),
_have_pastix_bug_small_matrix(false),
_a_when_bug()
{
load (a, opt);
}
template<class T, class M>
void
solver_pastix_base_rep<T,M>::update_values (const csr<T,M>& a)
{
if (_n == 0) return;
check_macro (size_t(_n) == a.nrow() && size_t(_n) == a.ncol(),
"update: local input matrix size distribution mismatch: ("<<a.nrow()<<","<<a.ncol()<<"), expect ("
<< _n << "," << _n << ")");
size_t nnz_a;
if (!_is_sym) {
nnz_a = a.nnz();
} else {
// conserve only the lower part of the csc(pastix) = the upper past of the csr(rheolef)
nnz_a = a.nnz() - (a.nnz() - a.nrow())/2; // TODO: fix when zero on diag
}
check_macro (size_t(_nnz) == nnz_a,
"update: local input matrix nnz distribution mismatch: nnz="<<nnz_a<<", expect nnz="<<_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,M>::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,M>::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 (_is_sym && dis_i > dis_j) continue;
_val [bp] = (*ap).second;
bp++;
}
}
numeric_factorization();
}
template<class T, class M>
vec<T,M>
solver_pastix_base_rep<T,M>::trans_solve (const vec<T,M>& rhs) const
{
if (_n == 0) return rhs;
// ================================================================
// solve
// ================================================================
check_macro (rhs.size() == size_t(_n), "invalid rhs size="<<rhs.size()<<": expect size="<<_n);
_new_rhs.resize (_n);
for (pastix_int_t i = 0; i < _n; i++) {
_new_rhs [i] = rhs [i];
}
// tasks:
// 5. numerical solve
// 6. numerical refinement
// 7. clean
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 rhs pointer: when new_n=0, but vector::begin() retuns a NULL...
_iparm[IPARM_START_TASK] = API_TASK_SOLVE;
_iparm[IPARM_END_TASK] = API_TASK_REFINE;
d_pastix (&_pastix_data,
0,
_n,
_ptr.begin().operator->(),
(_idx.begin().operator->() != NULL) ? _idx.begin().operator->() : &itmp,
(_val.begin().operator->() != NULL) ? _val.begin().operator->() : &vtmp,
(_i2new_dis_i.begin().operator->() != NULL) ? _i2new_dis_i.begin().operator->() : &itmp,
NULL,
(_n != 0) ? _new_rhs.begin().operator->() : &vtmp,
1,
_iparm,
_dparm);
// new_i2dis_i_base [new_i] - base = new_i2dis_i [new_i]
vec<T,M> x (_csr_row_ownership);
for (pastix_int_t i = 0; i < _n; i++) {
x [i] = _new_rhs [i];
}
return x;
}
template<class T, class M>
vec<T,M>
solver_pastix_base_rep<T,M>::solve (const vec<T,M>& rhs) const
{
// TODO: make a csc<T> wrapper around csr<T> and use csc<T> in form & form_assembly
// => avoid the transposition
if (! _is_sym) {
warning_macro ("solve: not yet supported for unsymmetric matrix");
}
return trans_solve (rhs);
}
template<class T, class M>
solver_pastix_base_rep<T,M>::~solver_pastix_base_rep()
{
if (_pastix_data == 0) return;
// tasks:
// 7. clean
_iparm[IPARM_START_TASK] = API_TASK_CLEAN;
_iparm[IPARM_END_TASK] = API_TASK_CLEAN;
d_pastix (&_pastix_data,
0,
_n,
0,
0,
0,
NULL,
NULL,
_new_rhs.begin().operator->(),
1,
_iparm,
_dparm);
// was allocated by d_cscd_redispatch()
_pastix_data == 0;
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
// TODO: code is only valid here for T=double (d_pastix etc)
template class solver_pastix_base_rep<double,sequential>;
#ifdef _RHEOLEF_HAVE_MPI
template class solver_pastix_base_rep<double,distributed>;
#endif // _RHEOLEF_HAVE_MPI
} // namespace rheolef
#endif // _RHEOLEF_HAVE_PASTIX
rheolef::solver_option::decide
static const long int decide
Definition solver_option.h:158
trace_macro
#define trace_macro(message)
Definition dis_macros.h:111
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
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
M
Expr1::memory_type M
Definition vec_expr_v2.h:416