Rheolef  7.2
an efficient C++ finite element environment
solver_pastix.h
Go to the documentation of this file.
1 #ifndef _RHEOLEF_SOLVER_PASTIX_H
2 #define _RHEOLEF_SOLVER_PASTIX_H
23 // solver pastix implementation: interface
24 //
25 #include "rheolef/config.h"
26 
27 #ifdef _RHEOLEF_HAVE_PASTIX
28 
29 #include "rheolef/solver.h"
30 extern "C" {
31 #define COMPLEXFLOAT_ /* workarroud a compile problem here */
32 #define COMPLEXDOUBLE_
33 #ifndef _RHEOLEF_HAVE_MPI
34 #define FORCE_NOMPI
35 #define MPI_Comm int
36 #endif // _RHEOLEF_HAVE_MPI
37 #include "pastix.h"
38 #include "cscd_utils.h"
39 #undef COMPLEXFLOAT_
40 #undef COMPLEXDOUBLE_
41 #undef FORCE_NOMPI
42 }
43 namespace rheolef {
44 
45 // =======================================================================
46 // pastix_base_rep
47 // =======================================================================
48 template<class T, class M>
49 class solver_pastix_base_rep : public solver_abstract_rep<T,M> {
50 public:
51 
52 // allocator:
53 
54  solver_pastix_base_rep ();
55  explicit solver_pastix_base_rep (const csr<T,M>& a, const solver_option& opt = solver_option());
56  bool initialized() const { return true; }
57  void update_values (const csr<T,M>& a);
58  ~solver_pastix_base_rep ();
59 
60 // accessors:
61 
62  vec<T,M> trans_solve (const vec<T,M>& rhs) const;
63  vec<T,M> solve (const vec<T,M>& rhs) const;
64 
65 // internal accessors & modifiers:
66 protected:
67 
68  void load (const csr<T,M>& a, const solver_option& opt = solver_option());
69  bool is_symmetric () const { return _is_sym; }
70  void resize (pastix_int_t n, pastix_int_t nnz);
71  void load_symmetric (const csr<T,M>& a);
72  void load_unsymmetric (const csr<T,M>& a);
73  void load_both_continued (const csr<T,M>& a);
74  void check () const;
75  void symbolic_factorization ();
76  void numeric_factorization ();
77 
78 // data:
79 public: // TODO: protected
80 //protected:
81  static const pastix_int_t _base = 1;
82  pastix_int_t _n;
83  pastix_int_t _nnz;
84  mutable std::vector<pastix_int_t> _ptr;
85  mutable std::vector<pastix_int_t> _idx; // row index, in csc format: dis_i = idx[p], p=0..nnz-1
86  mutable std::vector<T> _val;
87  bool _is_sym;
88  size_t _pattern_dimension;
89  mutable pastix_data_t* _pastix_data; // Pointer to a storage structure needed by pastix
90  mutable pastix_int_t _iparm [IPARM_SIZE]; // integer parameters for pastix
91  mutable T _dparm [DPARM_SIZE]; // floating parameters for pastix
92  distributor _csr_row_ownership;
93  distributor _csr_col_ownership;
94  solver_option _opt;
95  mutable std::vector<T> _new_rhs;
96  mutable std::vector<pastix_int_t> _new_i2dis_i_base;
97  mutable std::vector<pastix_int_t> _i2new_dis_i; // permutation
98  bool _have_pastix_bug_small_matrix; // circumvent when np < a.dis_nrow
99  csr<T,M> _a_when_bug; // use it when pastix bug (too small)
100 
101 // internal:
102  pastix_data_t** pp_data() const { return &_pastix_data; }
103 
104 };
105 // =======================================================================
106 // pastix_rep
107 // =======================================================================
108 template<class T,class M>
109 class solver_pastix_rep {};
110 
111 // ====================================================================
112 // sequential pastix interface
113 // ====================================================================
114 template<class T>
115 class solver_pastix_rep<T,sequential> : public solver_pastix_base_rep<T,sequential> {
116 public:
117  typedef solver_pastix_base_rep<T,sequential> base;
118 
119 // allocator:
120 
121  solver_pastix_rep () : base() {}
122  explicit solver_pastix_rep (const csr<T,sequential>& a, const solver_option& opt = solver_option())
123  : base(a,opt) {}
124  void update_values (const csr<T,sequential>& a) { base::update_values(a); }
125 
126 // accessors:
127 
128  vec<T,sequential> trans_solve (const vec<T,sequential>& rhs) const { return base::trans_solve(rhs); }
129  vec<T,sequential> solve (const vec<T,sequential>& rhs) const { return base::solve(rhs); }
130 
131 // internal accessors & modifiers:
132 protected:
133 
134  void load (const csr<T,sequential>& a, const solver_option& opt = solver_option()) { load(a,opt); }
135  bool is_symmetric () const { return base::is_symmetric(); }
136  void resize (pastix_int_t n, pastix_int_t nnz) { resize (n,nnz); }
137  void load_symmetric (const csr<T,sequential>& a) { load_symmetric(a); }
138  void load_unsymmetric (const csr<T,sequential>& a) { load_unsymmetric(a); }
139  void load_both_continued (const csr<T,sequential>& a) { load_both_continued(a); }
140  void check () const { check(); }
141  void symbolic_factorization () { symbolic_factorization(); }
142  void numeric_factorization () { numeric_factorization(); }
143 };
144 // ====================================================================
145 // distributed pastix interface
146 // ====================================================================
147 #ifdef _RHEOLEF_HAVE_MPI
148 template<class T>
149 class solver_pastix_rep<T,distributed> : public solver_pastix_base_rep<T,distributed> {
150  typedef solver_pastix_base_rep<T,distributed> base;
151 
152 public:
153 
154 // allocator:
155 
156  solver_pastix_rep ();
157  explicit solver_pastix_rep (const csr<T,distributed>& a, const solver_option& opt = solver_option());
158  void update_values (const csr<T,distributed>& a);
159  ~solver_pastix_rep ();
160 
161 // accessors:
162 
163  const communicator& comm () const { return _comm; }
164  vec<T,distributed> trans_solve (const vec<T,distributed>& rhs) const;
165  vec<T,distributed> solve (const vec<T,distributed>& rhs) const;
166 
167 // internal accessors & modifiers:
168 protected:
169 
170  void load (const csr<T,distributed>& a, const solver_option& opt = solver_option());
171  bool is_symmetric () const { return base::_is_sym; }
172  void resize (pastix_int_t n, pastix_int_t nnz);
173  void load_symmetric (const csr<T,distributed>& a);
174  void load_unsymmetric (const csr<T,distributed>& a);
175  void load_both_continued (const csr<T,distributed>& a);
176  void check () const;
177  void symbolic_factorization ();
178  void numeric_factorization ();
179 
180 // data:
181 protected:
182  communicator _comm;
183  std::vector<pastix_int_t> _i2dis_i_base; // dis_j = i2dis_i_base[j] - base, j=0..n-1
184  pastix_int_t _new_n; // new re-ordering
185  pastix_int_t* _new_ptr;
186  pastix_int_t* _new_idx;
187  T* _new_val;
188 };
189 #endif // _RHEOLEF_HAVE_MPI
190 
191 } // namespace rheolef
192 #endif // _RHEOLEF_HAVE_PASTIX
193 #endif // _RHEOLEF_SOLVER_PASTIX_H
see the communicator page for the full documentation
distributed
Definition: asr.cc:228
Expr1::float_type T
Definition: field_expr.h:230
This file is part of Rheolef.
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
void check(Float p, const field &uh)
void load(idiststream &in, Float &p, field &uh)