23#include "rheolef/config.h"
24#ifdef _RHEOLEF_HAVE_PASTIX
26#include "rheolef/cg.h"
27#include "rheolef/eye.h"
32template<
class T,
class M>
34solver_pastix_base_rep<T,M>::resize (pastix_int_t n, pastix_int_t nnz)
42template<
class T,
class M>
44solver_pastix_base_rep<T,M>::check ()
const
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);
63template<
class T,
class M>
65solver_pastix_base_rep<T,M>::load_both_continued (
const csr<T,M>& a)
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();
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;
81 _idx [bp] = dis_j + _base;
84 _ptr [i+1] = bp + _base;
86 check_macro (bp == _nnz,
"factorization: invalid nnz count");
89template<
class T,
class M>
91solver_pastix_base_rep<T,M>::load_unsymmetric (
const csr<T,M>& a)
96 load_both_continued (a);
98template<
class T,
class M>
100solver_pastix_base_rep<T,M>::load_symmetric (
const csr<T,M>& a)
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++;
118warning_macro (
"load sym(1): nnz_upper_dia="<<nnz_upper_dia<<
"...");
119 resize (
a.nrow(), nnz_upper_dia);
121 load_both_continued (a);
124template<
class T,
class M>
126solver_pastix_base_rep<T,M>::symbolic_factorization ()
129 if (_have_pastix_bug_small_matrix) {
133 const pastix_int_t nbthread = 1;
134 const pastix_int_t ordering = 0;
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,
144 _ptr.begin().operator->(),
145 _idx.begin().operator->(),
155 _iparm[IPARM_THREAD_NBR] = nbthread;
156 if (is_symmetric()) {
157 _iparm[IPARM_SYM] = API_SYM_YES;
158 _iparm[IPARM_FACTORIZATION] = API_FACT_LDLT;
160 _iparm[IPARM_SYM] = API_SYM_NO;
161 _iparm[IPARM_FACTORIZATION] = API_FACT_LU;
163 _iparm[IPARM_MATRIX_VERIFICATION] = API_NO;
164 _iparm[IPARM_VERBOSE] = _opt.verbose_level;
165 _iparm[IPARM_ORDERING] = ordering;
168 do_incomplete = (_opt.iterative != 0);
170 do_incomplete = (_pattern_dimension > 2);
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;
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;
182 _i2new_dis_i.resize (_n);
184 pastix_int_t itmp = 0;
191 _iparm[IPARM_START_TASK] = API_TASK_ORDERING;
192 _iparm[IPARM_END_TASK] = API_TASK_ANALYSE;
194 std::vector<pastix_int_t> new_i2i (_n);
195 std::vector<double> dummy_rhs (_n);
196 d_pastix (&_pastix_data,
199 _ptr.begin().operator->(),
200 (_idx.begin().operator->() != NULL) ? _idx.begin().operator->() : &itmp,
202 (_i2new_dis_i.begin().operator->() != NULL) ? _i2new_dis_i.begin().operator->() : &itmp,
203 (new_i2i.begin().operator->() != NULL) ? new_i2i.begin().operator->() : &itmp,
209template<
class T,
class M>
211solver_pastix_base_rep<T,M>::numeric_factorization ()
214 if (_have_pastix_bug_small_matrix) {
220 pastix_int_t itmp = 0;
222 _iparm[IPARM_START_TASK] = API_TASK_NUMFACT;
223 _iparm[IPARM_END_TASK] = API_TASK_NUMFACT;
224 d_pastix (&_pastix_data,
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,
237template<
class T,
class M>
239solver_pastix_base_rep<T,M>::load (
const csr<T,M>& a,
const solver_option& opt)
242#define _RHEOLEF_SEQ_PASTIX_HAVE_SYM_BUG
243#ifdef _RHEOLEF_SEQ_PASTIX_HAVE_SYM_BUG
246 _is_sym =
a.is_symmetric();
248 _pattern_dimension =
a.pattern_dimension();
249 _csr_row_ownership =
a.row_ownership();
250 _csr_col_ownership =
a.col_ownership();
253 check_macro (
a.nrow() ==
a.ncol(),
"factorization: only square matrix are supported");
264 symbolic_factorization ();
266 numeric_factorization();
269template<
class T,
class M>
270solver_pastix_base_rep<T,M>::solver_pastix_base_rep ()
271 : solver_abstract_rep<
T,
M>(solver_option()),
278 _pattern_dimension(0),
282 _csr_row_ownership(),
283 _csr_col_ownership(),
288 _have_pastix_bug_small_matrix(false),
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),
301 _pattern_dimension(0),
305 _csr_row_ownership(),
306 _csr_col_ownership(),
311 _have_pastix_bug_small_matrix(false),
316template<
class T,
class M>
318solver_pastix_base_rep<T,M>::update_values (
const csr<T,M>& a)
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 <<
")");
329 nnz_a =
a.nnz() - (
a.nnz() -
a.nrow())/2;
332 "update: local input matrix nnz distribution mismatch: nnz="<<nnz_a<<
", expect nnz="<<_nnz);
334 size_t first_dis_i =
a.row_ownership().first_index();
335 size_t first_dis_j =
a.col_ownership().first_index();
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;
348 numeric_factorization();
350template<
class T,
class M>
352solver_pastix_base_rep<T,M>::trans_solve (
const vec<T,M>& rhs)
const
354 if (_n == 0)
return rhs;
358 check_macro (rhs.size() ==
size_t(_n),
"invalid rhs size="<<rhs.size()<<
": expect size="<<_n);
360 _new_rhs.resize (_n);
361 for (pastix_int_t i = 0; i < _n; i++) {
362 _new_rhs [i] = rhs [i];
368 pastix_int_t itmp = 0;
370 _iparm[IPARM_START_TASK] = API_TASK_SOLVE;
371 _iparm[IPARM_END_TASK] = API_TASK_REFINE;
372 d_pastix (&_pastix_data,
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,
380 (_n != 0) ? _new_rhs.begin().operator->() : &vtmp,
386 vec<T,M> x (_csr_row_ownership);
387 for (pastix_int_t i = 0; i < _n; i++) {
388 x [i] = _new_rhs [i];
392template<
class T,
class M>
394solver_pastix_base_rep<T,M>::solve (
const vec<T,M>& rhs)
const
399 warning_macro (
"solve: not yet supported for unsymmetric matrix");
401 return trans_solve (rhs);
403template<
class T,
class M>
404solver_pastix_base_rep<T,M>::~solver_pastix_base_rep()
406 if (_pastix_data == 0)
return;
409 _iparm[IPARM_START_TASK] = API_TASK_CLEAN;
410 _iparm[IPARM_END_TASK] = API_TASK_CLEAN;
411 d_pastix (&_pastix_data,
419 _new_rhs.begin().operator->(),
432template class solver_pastix_base_rep<double,sequential>;
434#ifdef _RHEOLEF_HAVE_MPI
435template class solver_pastix_base_rep<double,distributed>;
static const long int decide
#define trace_macro(message)
#define warning_macro(message)
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.
void check(Float p, const field &uh)
void load(idiststream &in, Float &p, field &uh)