Rheolef  7.1
an efficient C++ finite element environment
solver_pastix_seq.cc
Go to the documentation of this file.
1 // c++ & boost::mpi adaptation of pastix "simple.c" sequential example
22 //
23 #include "rheolef/config.h"
24 #ifdef _RHEOLEF_HAVE_PASTIX
25 #include "solver_pastix.h"
26 #include "rheolef/cg.h"
27 #include "rheolef/eye.h"
28 
29 namespace rheolef {
30 using namespace std;
31 
32 template<class T, class M>
33 void
34 solver_pastix_base_rep<T,M>::resize (pastix_int_t n, pastix_int_t nnz)
35 {
36  _n = n;
37  _nnz = nnz;
38  _ptr.resize(n+1);
39  _idx.resize(nnz);
40  _val.resize(nnz);
41 }
42 template<class T, class M>
43 void
45 {
46  trace_macro ("check...");
47  if (_n == 0) return;
54  pastix_int_t symmetry = (is_symmetric() ? API_SYM_YES : API_SYM_NO);
55  pastix_int_t *ptr_begin = (pastix_int_t*) _ptr.begin().operator->();
56  pastix_int_t *idx_begin = (pastix_int_t*) _idx.begin().operator->();
57  T *val_begin = (T*) _val.begin().operator->();
58  pastix_int_t status = d_pastix_checkMatrix (0, _opt.verbose_level, symmetry, API_YES,
59  _n, &ptr_begin, &idx_begin, &val_begin, NULL, 1);
60  check_macro (status == 0, "pastix check returns error status = " << status);
61  trace_macro ("pastix matrix check: ok");
62 }
63 template<class T, class M>
64 void
65 solver_pastix_base_rep<T,M>::load_both_continued (const csr<T,M>& a)
66 {
67 warning_macro ("load both...");
68  size_t first_dis_i = a.row_ownership().first_index();
69  size_t first_dis_j = a.col_ownership().first_index();
70  typename csr<T,M>::const_iterator aptr = a.begin();
71  pastix_int_t bp = 0;
72  _ptr [0] = bp + _base;
73  for (size_t i = 0; i < a.nrow(); i++) {
74  size_t dis_i = first_dis_i + i;
75  for (typename csr<T,M>::const_data_iterator ap = aptr[i]; ap < aptr[i+1]; ap++) {
76  size_t j = (*ap).first;
77  const T& val = (*ap).second;
78  size_t dis_j = first_dis_j + j;
79  if (_is_sym && dis_i > dis_j) continue;
80  _val [bp] = val;
81  _idx [bp] = dis_j + _base;
82  bp++;
83  }
84  _ptr [i+1] = bp + _base;
85  }
86  check_macro (bp == _nnz, "factorization: invalid nnz count");
87 warning_macro ("load both done");
88 }
89 template<class T, class M>
90 void
91 solver_pastix_base_rep<T,M>::load_unsymmetric (const csr<T,M>& a)
92 {
93  size_t n = a.nrow();
94  size_t nnz = a.nnz();
95  resize (n, nnz);
96  load_both_continued (a);
97 }
98 template<class T, class M>
99 void
100 solver_pastix_base_rep<T,M>::load_symmetric (const csr<T,M>& a)
101 {
102 warning_macro ("load sym...");
103  // count nz entries in the upper+dia part
104  size_t nnz_upper_dia = 0;
105  size_t first_dis_i = a.row_ownership().first_index();
106  size_t first_dis_j = a.col_ownership().first_index();
107  typename csr<T,M>::const_iterator aptr = a.begin();
108  for (size_t i = 0, n = a.nrow(); i < n; i++) {
109  size_t dis_i = first_dis_i + i;
110  for (typename csr<T,M>::const_data_iterator ap = aptr[i]; ap < aptr[i+1]; ap++) {
111  size_t j = (*ap).first;
112  size_t dis_j = first_dis_j + j;
113  if (dis_i <= dis_j) nnz_upper_dia++;
114  }
115  }
116 warning_macro ("load sym(1): nrow="<<a.nrow()<<"...");
117 warning_macro ("load sym(1): nnz="<<a.nnz()<<"...");
118 warning_macro ("load sym(1): nnz_upper_dia="<<nnz_upper_dia<<"...");
119  resize (a.nrow(), nnz_upper_dia);
120 warning_macro ("load sym(2)...");
121  load_both_continued (a);
122 warning_macro ("load sym done");
123 }
124 template<class T, class M>
125 void
126 solver_pastix_base_rep<T,M>::symbolic_factorization ()
127 {
128  if (_n == 0) return;
129  if (_have_pastix_bug_small_matrix) {
130  trace_macro ("sym_fact: bug, circumvent !");
131  return;
132  }
133  const pastix_int_t nbthread = 1; // threads are not yet very well supported with openmpi & scotch
134  const pastix_int_t ordering = 0; // scotch
135 
136  // tasks :
137  // 0. set params to default values
138  _iparm[IPARM_START_TASK] = API_TASK_INIT;
139  _iparm[IPARM_END_TASK] = API_TASK_INIT;
140  _iparm[IPARM_MODIFY_PARAMETER] = API_NO;
141  d_pastix (&_pastix_data,
142  0,
143  _n,
144  _ptr.begin().operator->(),
145  _idx.begin().operator->(),
146  NULL, // _val.begin().operator->(),
147  NULL,
148  NULL,
149  NULL,
150  0,
151  _iparm,
152  _dparm);
153 
154  // Customize some parameters
155  _iparm[IPARM_THREAD_NBR] = nbthread;
156  if (is_symmetric()) {
157  _iparm[IPARM_SYM] = API_SYM_YES;
158  _iparm[IPARM_FACTORIZATION] = API_FACT_LDLT;
159  } else {
160  _iparm[IPARM_SYM] = API_SYM_NO;
161  _iparm[IPARM_FACTORIZATION] = API_FACT_LU;
162  }
163  _iparm[IPARM_MATRIX_VERIFICATION] = API_NO;
164  _iparm[IPARM_VERBOSE] = _opt.verbose_level;
165  _iparm[IPARM_ORDERING] = ordering;
166  bool do_incomplete;
167  if (_opt.iterative != solver_option::decide) {
168  do_incomplete = (_opt.iterative != 0);
169  } else {
170  do_incomplete = (_pattern_dimension > 2); // 3d-pattern => iterative & IC(k) precond
171  }
172  _iparm[IPARM_INCOMPLETE] = (do_incomplete ? 1 : 0);
173  _iparm[IPARM_OOC_LIMIT] = _opt.ooc;
174  if (_opt.iterative == 1) {
175  _dparm[DPARM_EPSILON_REFINEMENT] = _opt.tol;
176  }
177  _iparm[IPARM_LEVEL_OF_FILL] = _opt.level_of_fill;
178  _iparm[IPARM_AMALGAMATION_LEVEL] = _opt.amalgamation;
179  _iparm[IPARM_RHS_MAKING] = API_RHS_B;
180  // 1) get the i2new_dis_i array: its indexes are in the [0:dis_n[ range
181  // i2new_dis_i[i] = i2new_dis_i [i]
182  _i2new_dis_i.resize (_n);
183 
184  pastix_int_t itmp = 0; // pastix does not manage NULL pointers: when n=0, vector::begin() retuns a NULL...
185  T vtmp = 0; // pastix does not manage NULL pointers: when n=0, vector::begin() retuns a NULL...
186 
187  // tasks :
188  // 1. ordering
189  // 2. symbolic factorization
190  // 3. tasks mapping and scheduling
191  _iparm[IPARM_START_TASK] = API_TASK_ORDERING;
192  _iparm[IPARM_END_TASK] = API_TASK_ANALYSE;
193 
194  std::vector<pastix_int_t> new_i2i (_n);
195  std::vector<double> dummy_rhs (_n);
196  d_pastix (&_pastix_data,
197  0,
198  _n,
199  _ptr.begin().operator->(),
200  (_idx.begin().operator->() != NULL) ? _idx.begin().operator->() : &itmp,
201  0,
202  (_i2new_dis_i.begin().operator->() != NULL) ? _i2new_dis_i.begin().operator->() : &itmp,
203  (new_i2i.begin().operator->() != NULL) ? new_i2i.begin().operator->() : &itmp,
204  0,
205  0,
206  _iparm,
207  _dparm);
208 }
209 template<class T, class M>
210 void
211 solver_pastix_base_rep<T,M>::numeric_factorization ()
212 {
213  if (_n == 0) return;
214  if (_have_pastix_bug_small_matrix) {
215  trace_macro ("num_fact: bug, circumvent !");
216  return;
217  }
218  // pastix tasks:
219  // 4. numerical factorization
220  pastix_int_t itmp = 0; // pastix does not manage NULL pointers: when n=0, vector::begin() retuns a NULL...
221  T vtmp = 0; // pastix does not manage NULL rhs pointer: when new_n=0, but vector::begin() retuns a NULL...
222  _iparm[IPARM_START_TASK] = API_TASK_NUMFACT;
223  _iparm[IPARM_END_TASK] = API_TASK_NUMFACT;
224  d_pastix (&_pastix_data,
225  0,
226  _n,
227  _ptr.begin().operator->(),
228  (_idx.begin().operator->() != NULL) ? _idx.begin().operator->() : &itmp,
229  (_val.begin().operator->() != NULL) ? _val.begin().operator->() : &vtmp,
230  (_i2new_dis_i.begin().operator->() != NULL) ? _i2new_dis_i.begin().operator->() : &itmp,
231  NULL,
232  NULL,
233  0,
234  _iparm,
235  _dparm);
236 }
237 template<class T, class M>
238 void
239 solver_pastix_base_rep<T,M>::load (const csr<T,M>& a, const solver_option& opt)
240 {
241 warning_macro ("load...");
242 #define _RHEOLEF_SEQ_PASTIX_HAVE_SYM_BUG
243 #ifdef _RHEOLEF_SEQ_PASTIX_HAVE_SYM_BUG
244  _is_sym = false;
245 #else // ! _RHEOLEF_SEQ_PASTIX_HAVE_SYM_BUG
246  _is_sym = a.is_symmetric();
247 #endif // ! _RHEOLEF_SEQ_PASTIX_HAVE_SYM_BUG
248  _pattern_dimension = a.pattern_dimension();
249  _csr_row_ownership = a.row_ownership();
250  _csr_col_ownership = a.col_ownership();
251  _opt = opt;
252 
253  check_macro (a.nrow() == a.ncol(), "factorization: only square matrix are supported");
254 
255  if (_is_sym) {
256  load_symmetric(a);
257  } else {
258  load_unsymmetric(a);
259  }
260  if (_opt.do_check) {
261  check ();
262  }
263 warning_macro ("symbfact...");
264  symbolic_factorization ();
265 warning_macro ("numbfact...");
266  numeric_factorization();
267 warning_macro ("load done");
268 }
269 template<class T, class M>
270 solver_pastix_base_rep<T,M>::solver_pastix_base_rep ()
271  : solver_abstract_rep<T,M>(solver_option()),
272  _n(),
273  _nnz(),
274  _ptr(),
275  _idx(),
276  _val(),
277  _is_sym(false),
278  _pattern_dimension(0),
279  _pastix_data(0),
280  _iparm(),
281  _dparm(),
282  _csr_row_ownership(),
283  _csr_col_ownership(),
284  _opt(),
285  _new_rhs(),
286  _new_i2dis_i_base(),
287  _i2new_dis_i(),
288  _have_pastix_bug_small_matrix(false),
289  _a_when_bug()
290 {
291 }
292 template<class T, class M>
293 solver_pastix_base_rep<T,M>::solver_pastix_base_rep (const csr<T,M>& a, const solver_option& opt)
294  : solver_abstract_rep<T,M>(opt),
295  _n(),
296  _nnz(),
297  _ptr(),
298  _idx(),
299  _val(),
300  _is_sym(false),
301  _pattern_dimension(0),
302  _pastix_data(0),
303  _iparm(),
304  _dparm(),
305  _csr_row_ownership(),
306  _csr_col_ownership(),
307  _opt(),
308  _new_rhs(),
309  _new_i2dis_i_base(),
310  _i2new_dis_i(),
311  _have_pastix_bug_small_matrix(false),
312  _a_when_bug()
313 {
314  load (a, opt);
315 }
316 template<class T, class M>
317 void
318 solver_pastix_base_rep<T,M>::update_values (const csr<T,M>& a)
319 {
320  if (_n == 0) return;
321  check_macro (size_t(_n) == a.nrow() && size_t(_n) == a.ncol(),
322  "update: local input matrix size distribution mismatch: ("<<a.nrow()<<","<<a.ncol()<<"), expect ("
323  << _n << "," << _n << ")");
324  size_t nnz_a;
325  if (!_is_sym) {
326  nnz_a = a.nnz();
327  } else {
328  // conserve only the lower part of the csc(pastix) = the upper past of the csr(rheolef)
329  nnz_a = a.nnz() - (a.nnz() - a.nrow())/2; // TODO: fix when zero on diag
330  }
331  check_macro (size_t(_nnz) == nnz_a,
332  "update: local input matrix nnz distribution mismatch: nnz="<<nnz_a<<", expect nnz="<<_nnz);
333 
334  size_t first_dis_i = a.row_ownership().first_index();
335  size_t first_dis_j = a.col_ownership().first_index();
336  pastix_int_t bp = 0;
337  typename csr<T,M>::const_iterator aptr = a.begin();
338  for (size_t i = 0; i < a.nrow(); i++) {
339  size_t dis_i = first_dis_i + i;
340  for (typename csr<T,M>::const_data_iterator ap = aptr[i]; ap < aptr[i+1]; ap++) {
341  size_t j = (*ap).first;
342  size_t dis_j = first_dis_j + j;
343  if (_is_sym && dis_i > dis_j) continue;
344  _val [bp] = (*ap).second;
345  bp++;
346  }
347  }
348  numeric_factorization();
349 }
350 template<class T, class M>
351 vec<T,M>
352 solver_pastix_base_rep<T,M>::trans_solve (const vec<T,M>& rhs) const
353 {
354  if (_n == 0) return rhs;
355  // ================================================================
356  // solve
357  // ================================================================
358  check_macro (rhs.size() == size_t(_n), "invalid rhs size="<<rhs.size()<<": expect size="<<_n);
359 
360  _new_rhs.resize (_n);
361  for (pastix_int_t i = 0; i < _n; i++) {
362  _new_rhs [i] = rhs [i];
363  }
364  // tasks:
365  // 5. numerical solve
366  // 6. numerical refinement
367  // 7. clean
368  pastix_int_t itmp = 0; // pastix does not manage NULL pointers: when n=0, vector::begin() retuns a NULL...
369  T vtmp = 0; // pastix does not manage NULL rhs pointer: when new_n=0, but vector::begin() retuns a NULL...
370  _iparm[IPARM_START_TASK] = API_TASK_SOLVE;
371  _iparm[IPARM_END_TASK] = API_TASK_REFINE;
372  d_pastix (&_pastix_data,
373  0,
374  _n,
375  _ptr.begin().operator->(),
376  (_idx.begin().operator->() != NULL) ? _idx.begin().operator->() : &itmp,
377  (_val.begin().operator->() != NULL) ? _val.begin().operator->() : &vtmp,
378  (_i2new_dis_i.begin().operator->() != NULL) ? _i2new_dis_i.begin().operator->() : &itmp,
379  NULL,
380  (_n != 0) ? _new_rhs.begin().operator->() : &vtmp,
381  1,
382  _iparm,
383  _dparm);
384 
385  // new_i2dis_i_base [new_i] - base = new_i2dis_i [new_i]
386  vec<T,M> x (_csr_row_ownership);
387  for (pastix_int_t i = 0; i < _n; i++) {
388  x [i] = _new_rhs [i];
389  }
390  return x;
391 }
392 template<class T, class M>
393 vec<T,M>
394 solver_pastix_base_rep<T,M>::solve (const vec<T,M>& rhs) const
395 {
396  // TODO: make a csc<T> wrapper around csr<T> and use csc<T> in form & form_assembly
397  // => avoid the transposition
398  if (! _is_sym) {
399  warning_macro ("solve: not yet supported for unsymmetric matrix");
400  }
401  return trans_solve (rhs);
402 }
403 template<class T, class M>
404 solver_pastix_base_rep<T,M>::~solver_pastix_base_rep()
405 {
406  if (_pastix_data == 0) return;
407  // tasks:
408  // 7. clean
409  _iparm[IPARM_START_TASK] = API_TASK_CLEAN;
410  _iparm[IPARM_END_TASK] = API_TASK_CLEAN;
411  d_pastix (&_pastix_data,
412  0,
413  _n,
414  0,
415  0,
416  0,
417  NULL,
418  NULL,
419  _new_rhs.begin().operator->(),
420  1,
421  _iparm,
422  _dparm);
423 
424  // was allocated by d_cscd_redispatch()
425  _pastix_data == 0;
426 }
427 // ----------------------------------------------------------------------------
428 // instanciation in library
429 // ----------------------------------------------------------------------------
430 // TODO: code is only valid here for T=double (d_pastix etc)
431 
432 template class solver_pastix_base_rep<double,sequential>;
433 
434 #ifdef _RHEOLEF_HAVE_MPI
435 template class solver_pastix_base_rep<double,distributed>;
436 #endif // _RHEOLEF_HAVE_MPI
437 
438 } // namespace rheolef
439 #endif // _RHEOLEF_HAVE_PASTIX
mkgeo_ball.n
int n
Definition: mkgeo_ball.sh:150
warning_macro
#define warning_macro(message)
Definition: dis_macros.h:53
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
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
check
void check(Float p, const field &uh)
Definition: p_laplacian_post.cc:49
mkgeo_contraction.status
status
Definition: mkgeo_contraction.sh:290
rheolef::solver_option::decide
static const long int decide
Definition: solver_option.h:158
M
Expr1::memory_type M
Definition: vec_expr_v2.h:416
T
Expr1::float_type T
Definition: field_expr.h:261
trace_macro
#define trace_macro(message)
Definition: dis_macros.h:111