Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
solver_pastix_seq.cc
Go to the documentation of this file.
1
21// 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
29namespace rheolef {
30using namespace std;
31
32template<class T, class M>
33void
34solver_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}
42template<class T, class M>
43void
44solver_pastix_base_rep<T,M>::check () const
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}
63template<class T, class M>
64void
65solver_pastix_base_rep<T,M>::load_both_continued (const csr<T,M>& a)
66{
67warning_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");
87warning_macro ("load both done");
88}
89template<class T, class M>
90void
91solver_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}
98template<class T, class M>
99void
100solver_pastix_base_rep<T,M>::load_symmetric (const csr<T,M>& a)
101{
102warning_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 }
116warning_macro ("load sym(1): nrow="<<a.nrow()<<"...");
117warning_macro ("load sym(1): nnz="<<a.nnz()<<"...");
118warning_macro ("load sym(1): nnz_upper_dia="<<nnz_upper_dia<<"...");
119 resize (a.nrow(), nnz_upper_dia);
120warning_macro ("load sym(2)...");
121 load_both_continued (a);
122warning_macro ("load sym done");
123}
124template<class T, class M>
125void
126solver_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}
209template<class T, class M>
210void
211solver_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}
237template<class T, class M>
238void
239solver_pastix_base_rep<T,M>::load (const csr<T,M>& a, const solver_option& opt)
240{
241warning_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 }
263warning_macro ("symbfact...");
264 symbolic_factorization ();
265warning_macro ("numbfact...");
266 numeric_factorization();
267warning_macro ("load done");
268}
269template<class T, class M>
270solver_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}
292template<class T, class M>
293solver_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}
316template<class T, class M>
317void
318solver_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}
350template<class T, class M>
351vec<T,M>
352solver_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}
392template<class T, class M>
393vec<T,M>
394solver_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}
403template<class T, class M>
404solver_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
432template class solver_pastix_base_rep<double,sequential>;
433
434#ifdef _RHEOLEF_HAVE_MPI
435template class solver_pastix_base_rep<double,distributed>;
436#endif // _RHEOLEF_HAVE_MPI
437
438} // namespace rheolef
439#endif // _RHEOLEF_HAVE_PASTIX
static const long int decide
#define trace_macro(message)
Definition dis_macros.h:111
#define warning_macro(message)
Definition dis_macros.h:53
Expr1::float_type T
Definition field_expr.h:230
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
This file is part of Rheolef.
STL namespace.
void check(Float p, const field &uh)
void load(idiststream &in, Float &p, field &uh)
Expr1::memory_type M