rheolef  6.5
solver_pastix.h
Go to the documentation of this file.
1 #ifndef _RHEOLEF_SOLVER_PASTIX_H
2 #define _RHEOLEF_SOLVER_PASTIX_H
3 #include "rheolef/config.h"
4 
5 #ifdef _RHEOLEF_HAVE_PASTIX
6 
7 #include "rheolef/solver.h"
8 extern "C" {
9 #define COMPLEXFLOAT_ /* workarroud a compile problem here */
10 #define COMPLEXDOUBLE_
11 #ifndef _RHEOLEF_HAVE_MPI
12 #define FORCE_NOMPI
13 #define MPI_Comm int
14 #endif // _RHEOLEF_HAVE_MPI
15 #include "pastix.h"
16 #include "cscd_utils.h"
17 #undef COMPLEXFLOAT_
18 #undef COMPLEXDOUBLE_
19 #undef FORCE_NOMPI
20 }
21 namespace rheolef {
22 
23 template<class T, class M>
24 class solver_pastix_base_rep : public solver_abstract_rep<T,M> {
25 public:
26 
27 
28  solver_pastix_base_rep ();
29  explicit solver_pastix_base_rep (const csr<T,M>& a, const solver_option_type& opt = solver_option_type());
30  void update_values (const csr<T,M>& a);
31  ~solver_pastix_base_rep ();
32 
33 
34  vec<T,M> trans_solve (const vec<T,M>& rhs) const;
35  vec<T,M> solve (const vec<T,M>& rhs) const;
36 
37 protected:
38 
39  void load (const csr<T,M>& a, const solver_option_type& opt = solver_option_type());
40  bool is_symmetric () const { return _is_sym; }
41  void resize (pastix_int_t n, pastix_int_t nnz);
42  void load_symmetric (const csr<T,M>& a);
43  void load_unsymmetric (const csr<T,M>& a);
44  void load_both_continued (const csr<T,M>& a);
45  void check () const;
46  void symbolic_factorization ();
47  void numeric_factorization ();
48 
49 public: // TODO: protected
50  static const pastix_int_t _base = 1;
51  pastix_int_t _n;
52  pastix_int_t _nnz;
53  mutable std::vector<pastix_int_t> _ptr;
54  mutable std::vector<pastix_int_t> _idx; // row index, in csc format: dis_i = idx[p], p=0..nnz-1
55  mutable std::vector<T> _val;
56  bool _is_sym;
57  size_t _pattern_dimension;
58  mutable pastix_data_t* _pastix_data; // Pointer to a storage structure needed by pastix
59  mutable pastix_int_t _iparm [IPARM_SIZE]; // integer parameters for pastix
60  mutable T _dparm [DPARM_SIZE]; // floating parameters for pastix
61  distributor _csr_row_ownership;
62  distributor _csr_col_ownership;
63  solver_option_type _opt;
64  mutable std::vector<T> _new_rhs;
65  mutable std::vector<pastix_int_t> _new_i2dis_i_base;
66  mutable std::vector<pastix_int_t> _i2new_dis_i; // permutation
67  bool _have_pastix_bug_small_matrix; // circumvent when np < a.dis_nrow
68  csr<T,M> _a_when_bug; // use it when pastix bug (too small)
69 
70  pastix_data_t** pp_data() const { return &_pastix_data; }
71 
72 };
73 template<class T,class M>
74 class solver_pastix_rep {};
75 
76 template<class T>
77 class solver_pastix_rep<T,sequential> : public solver_pastix_base_rep<T,sequential> {
78 public:
79  typedef solver_pastix_base_rep<T,sequential> base;
80 
81 
82  solver_pastix_rep () : base() {}
83  explicit solver_pastix_rep (const csr<T,sequential>& a, const solver_option_type& opt = solver_option_type())
84  : base(a,opt) {}
85  void update_values (const csr<T,sequential>& a) { base::update_values(a); }
86 
87 
88  vec<T,sequential> trans_solve (const vec<T,sequential>& rhs) const { return base::trans_solve(rhs); }
89  vec<T,sequential> solve (const vec<T,sequential>& rhs) const { return base::solve(rhs); }
90 
91 protected:
92 
93  void load (const csr<T,sequential>& a, const solver_option_type& opt = solver_option_type()) { load(a,opt); }
94  bool is_symmetric () const { return base::is_symmetric(); }
95  void resize (pastix_int_t n, pastix_int_t nnz) { resize (n,nnz); }
96  void load_symmetric (const csr<T,sequential>& a) { load_symmetric(a); }
97  void load_unsymmetric (const csr<T,sequential>& a) { load_unsymmetric(a); }
98  void load_both_continued (const csr<T,sequential>& a) { load_both_continued(a); }
99  void check () const { check(); }
100  void symbolic_factorization () { symbolic_factorization(); }
101  void numeric_factorization () { numeric_factorization(); }
102 };
103 #ifdef _RHEOLEF_HAVE_MPI
104 template<class T>
105 class solver_pastix_rep<T,distributed> : public solver_pastix_base_rep<T,distributed> {
106  typedef solver_pastix_base_rep<T,distributed> base;
107 
108 public:
109 
110 
111  solver_pastix_rep ();
112  explicit solver_pastix_rep (const csr<T,distributed>& a, const solver_option_type& opt = solver_option_type());
113  void update_values (const csr<T,distributed>& a);
114  ~solver_pastix_rep ();
115 
116 
117  const communicator& comm () const { return _comm; }
118  vec<T,distributed> trans_solve (const vec<T,distributed>& rhs) const;
119  vec<T,distributed> solve (const vec<T,distributed>& rhs) const;
120 
121 protected:
122 
123  void load (const csr<T,distributed>& a, const solver_option_type& opt = solver_option_type());
124  bool is_symmetric () const { return base::_is_sym; }
125  void resize (pastix_int_t n, pastix_int_t nnz);
126  void load_symmetric (const csr<T,distributed>& a);
127  void load_unsymmetric (const csr<T,distributed>& a);
128  void load_both_continued (const csr<T,distributed>& a);
129  void check () const;
130  void symbolic_factorization ();
131  void numeric_factorization ();
132 
133 protected:
134  communicator _comm;
135  std::vector<pastix_int_t> _i2dis_i_base; // dis_j = i2dis_i_base[j] - base, j=0..n-1
136  pastix_int_t _new_n; // new re-ordering
137  pastix_int_t* _new_ptr;
138  pastix_int_t* _new_idx;
139  T* _new_val;
140 };
141 #endif // _RHEOLEF_HAVE_MPI
142 
143 } // namespace rheolef
144 #endif // _RHEOLEF_HAVE_PASTIX
145 #endif // _RHEOLEF_SOLVER_PASTIX_H
146