solver_pastix.h
Go to the documentation of this file.
#ifndef _RHEOLEF_SOLVER_PASTIX_H
#define _RHEOLEF_SOLVER_PASTIX_H
// solver pastix implementation: interface
//
#include "rheolef/config.h"
#ifdef _RHEOLEF_HAVE_PASTIX
#include "rheolef/solver.h"
extern "C" {
#define COMPLEXFLOAT_ /* workarroud a compile problem here */
#define COMPLEXDOUBLE_
#ifndef _RHEOLEF_HAVE_MPI
#define FORCE_NOMPI
#define MPI_Comm int
#endif // _RHEOLEF_HAVE_MPI
#include "pastix.h"
#include "cscd_utils.h"
#undef COMPLEXFLOAT_
#undef COMPLEXDOUBLE_
#undef FORCE_NOMPI
}
namespace rheolef {
// =======================================================================
// pastix_base_rep
// =======================================================================
template<class T, class M>
class solver_pastix_base_rep : public solver_abstract_rep<T,M> {
public:
// allocator:
solver_pastix_base_rep ();
explicit solver_pastix_base_rep (const csr<T,M>& a, const solver_option& opt = solver_option());
bool initialized() const { return true; }
void update_values (const csr<T,M>& a);
~solver_pastix_base_rep ();
// accessors:
vec<T,M> trans_solve (const vec<T,M>& rhs) const;
vec<T,M> solve (const vec<T,M>& rhs) const;
// internal accessors & modifiers:
protected:
void load (const csr<T,M>& a, const solver_option& opt = solver_option());
bool is_symmetric () const { return _is_sym; }
void resize (pastix_int_t n, pastix_int_t nnz);
void load_symmetric (const csr<T,M>& a);
void load_unsymmetric (const csr<T,M>& a);
void load_both_continued (const csr<T,M>& a);
void check () const;
void symbolic_factorization ();
void numeric_factorization ();
// data:
public: // TODO: protected
//protected:
static const pastix_int_t _base = 1;
pastix_int_t _n;
pastix_int_t _nnz;
mutable std::vector<pastix_int_t> _ptr;
mutable std::vector<pastix_int_t> _idx; // row index, in csc format: dis_i = idx[p], p=0..nnz-1
mutable std::vector<T> _val;
bool _is_sym;
size_t _pattern_dimension;
mutable pastix_data_t* _pastix_data; // Pointer to a storage structure needed by pastix
mutable pastix_int_t _iparm [IPARM_SIZE]; // integer parameters for pastix
mutable T _dparm [DPARM_SIZE]; // floating parameters for pastix
distributor _csr_row_ownership;
distributor _csr_col_ownership;
solver_option _opt;
mutable std::vector<T> _new_rhs;
mutable std::vector<pastix_int_t> _new_i2dis_i_base;
mutable std::vector<pastix_int_t> _i2new_dis_i; // permutation
bool _have_pastix_bug_small_matrix; // circumvent when np < a.dis_nrow
csr<T,M> _a_when_bug; // use it when pastix bug (too small)
// internal:
pastix_data_t** pp_data() const { return &_pastix_data; }
};
// =======================================================================
// pastix_rep
// =======================================================================
template<class T,class M>
class solver_pastix_rep {};
// ====================================================================
// sequential pastix interface
// ====================================================================
template<class T>
class solver_pastix_rep<T,sequential> : public solver_pastix_base_rep<T,sequential> {
public:
typedef solver_pastix_base_rep<T,sequential> base;
// allocator:
solver_pastix_rep () : base() {}
explicit solver_pastix_rep (const csr<T,sequential>& a, const solver_option& opt = solver_option())
: base(a,opt) {}
void update_values (const csr<T,sequential>& a) { base::update_values(a); }
// accessors:
vec<T,sequential> trans_solve (const vec<T,sequential>& rhs) const { return base::trans_solve(rhs); }
vec<T,sequential> solve (const vec<T,sequential>& rhs) const { return base::solve(rhs); }
// internal accessors & modifiers:
protected:
void load (const csr<T,sequential>& a, const solver_option& opt = solver_option()) { load(a,opt); }
bool is_symmetric () const { return base::is_symmetric(); }
void resize (pastix_int_t n, pastix_int_t nnz) { resize (n,nnz); }
void load_symmetric (const csr<T,sequential>& a) { load_symmetric(a); }
void load_unsymmetric (const csr<T,sequential>& a) { load_unsymmetric(a); }
void load_both_continued (const csr<T,sequential>& a) { load_both_continued(a); }
void check () const { check(); }
void symbolic_factorization () { symbolic_factorization(); }
void numeric_factorization () { numeric_factorization(); }
};
// ====================================================================
// distributed pastix interface
// ====================================================================
#ifdef _RHEOLEF_HAVE_MPI
template<class T>
class solver_pastix_rep<T,distributed> : public solver_pastix_base_rep<T,distributed> {
typedef solver_pastix_base_rep<T,distributed> base;
public:
// allocator:
solver_pastix_rep ();
explicit solver_pastix_rep (const csr<T,distributed>& a, const solver_option& opt = solver_option());
void update_values (const csr<T,distributed>& a);
~solver_pastix_rep ();
// accessors:
const communicator& comm () const { return _comm; }
vec<T,distributed> trans_solve (const vec<T,distributed>& rhs) const;
vec<T,distributed> solve (const vec<T,distributed>& rhs) const;
// internal accessors & modifiers:
protected:
void load (const csr<T,distributed>& a, const solver_option& opt = solver_option());
bool is_symmetric () const { return base::_is_sym; }
void resize (pastix_int_t n, pastix_int_t nnz);
void load_symmetric (const csr<T,distributed>& a);
void load_unsymmetric (const csr<T,distributed>& a);
void load_both_continued (const csr<T,distributed>& a);
void check () const;
void symbolic_factorization ();
void numeric_factorization ();
// data:
protected:
communicator _comm;
std::vector<pastix_int_t> _i2dis_i_base; // dis_j = i2dis_i_base[j] - base, j=0..n-1
pastix_int_t _new_n; // new re-ordering
pastix_int_t* _new_ptr;
pastix_int_t* _new_idx;
T* _new_val;
};
#endif // _RHEOLEF_HAVE_MPI
} // namespace rheolef
#endif // _RHEOLEF_HAVE_PASTIX
#endif // _RHEOLEF_SOLVER_PASTIX_H
a
Definition: diffusion_isotropic.h:25
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
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
load
void load(idiststream &in, Float &p, field &uh)
Definition: p_laplacian_post.cc:64
mkgeo_ball.n
n
Definition: mkgeo_ball.sh:150
check
void check(Float p, const field &uh)
Definition: p_laplacian_post.cc:49
rheolef::distributed
distributed
Definition: asr.cc:228
T
Expr1::float_type T
Definition: field_expr.h:218