48 #include "rheolef/config.h"
49 #ifdef _RHEOLEF_HAVE_MUMPS
52 #undef DEBUG_MUMPS_SCOTCH
53 #ifdef DEBUG_MUMPS_SCOTCH
54 #undef _RHEOLEF_HAVE_MUMPS_WITH_PTSCOTCH
65 template<
class T,
class M>
68 nnz_dia_upper (
const csr<T,M>&
a)
72 typename csr<T,M>::const_iterator dia_ia =
a.begin();
74 for (
typename csr<T,M>::const_data_iterator
p = dia_ia[i];
p < dia_ia[i+1];
p++) {
76 if (j >= i) dia_nnz++;
82 template<
class T,
class M>
85 nnz_ext_upper (
const csr<T,M>&
a)
89 #ifdef _RHEOLEF_HAVE_MPI
93 nnz_ext_upper (
const csr<T,distributed>&
a)
96 size_type first_i =
a.row_ownership().first_index();
98 typename csr<T,distributed>::const_iterator ext_ia =
a.ext_begin();
100 for (
typename csr<T,distributed>::const_data_iterator
p = ext_ia[i];
p < ext_ia[i+1];
p++) {
103 if (dis_j >= dis_i) ext_nnz++;
109 template<
class T,
class M>
112 csr2mumps_struct_seq (
const csr<T,M>&
a, vector<int>& row, vector<int>& col,
bool use_symmetry)
115 typename csr<T,M>::const_iterator dia_ia =
a.begin();
116 for (
size_type i = 0,
n =
a.nrow(), q = 0; i <
n; i++) {
117 for (
typename csr<T,M>::const_data_iterator
p = dia_ia[i];
p < dia_ia[i+1];
p++) {
119 if (!use_symmetry || j >= i) {
127 template<
class T,
class M>
130 csr2mumps_values_seq (
const csr<T,M>&
a, vector<T>& val,
bool use_symmetry)
133 typename csr<T,M>::const_iterator dia_ia =
a.begin();
134 for (
size_type i = 0,
n =
a.nrow(), q = 0; i <
n; i++) {
135 for (
typename csr<T,M>::const_data_iterator
p = dia_ia[i];
p < dia_ia[i+1];
p++) {
137 if (!use_symmetry || j >= i) {
138 val[q] = (*p).second;
144 #ifdef _RHEOLEF_HAVE_MPI
148 csr2mumps_struct_dis (
const csr<T,sequential>&
a, vector<int>& row, vector<int>& col,
bool use_symmetry,
bool)
155 csr2mumps_struct_dis (
const csr<T,distributed>&
a, vector<int>& row, vector<int>& col,
bool use_symmetry,
bool drop_ext_nnz)
158 typename csr<T,distributed>::const_iterator dia_ia =
a.begin();
159 typename csr<T,distributed>::const_iterator ext_ia =
a.ext_begin();
163 size_type first_i =
a.row_ownership().first_index();
164 for (
size_type i = 0,
n =
a.nrow(), q = 0; i <
n; i++) {
165 for (
typename csr<T,distributed>::const_data_iterator
p = dia_ia[i];
p < dia_ia[i+1];
p++) {
169 assert_macro (j < nc,
"invalid matrix");
170 assert_macro (dis_i < dis_nr && dis_j < dis_nc,
"invalid matrix");
171 if (!use_symmetry || dis_j >= dis_i) {
177 if (drop_ext_nnz)
continue;
178 for (
typename csr<T,distributed>::const_data_iterator
p = ext_ia[i];
p < ext_ia[i+1];
p++) {
181 assert_macro (dis_i < dis_nr && dis_j < dis_nc,
"invalid matrix");
182 if (!use_symmetry || dis_j >= dis_i) {
193 csr2mumps_values_dis (
const csr<T,sequential>&
a, vector<T>& val,
bool use_symmetry,
bool)
200 csr2mumps_values_dis (
const csr<T,distributed>&
a, vector<T>& val,
bool use_symmetry,
bool drop_ext_nnz)
203 std:fill (val.begin(), val.end(), std::numeric_limits<T>::max());
207 size_type first_i =
a.row_ownership().first_index();
208 typename csr<T,distributed>::const_iterator dia_ia =
a.begin();
209 typename csr<T,distributed>::const_iterator ext_ia =
a.ext_begin();
210 for (
size_type i = 0,
n =
a.nrow(), q = 0; i <
n; i++) {
211 for (
typename csr<T,distributed>::const_data_iterator
p = dia_ia[i];
p < dia_ia[i+1];
p++) {
214 if (!use_symmetry || dis_j >= dis_i) {
216 check_macro(q < val.size(),
"invalid index");
218 val[q] = (*p).second;
222 if (drop_ext_nnz)
continue;
223 for (
typename csr<T,distributed>::const_data_iterator
p = ext_ia[i];
p < ext_ia[i+1];
p++) {
226 if (!use_symmetry || dis_j >= dis_i) {
228 check_macro(q < val.size(),
"invalid index");
230 val[q] = (*p).second;
236 T val_min = std::numeric_limits<T>::max(), val_max = 0;
237 for (
size_t q = 0; q < val.size(); q++) {
238 val_min = std::min (val_min, abs(val[q]));
239 val_max = std::max (val_max, abs(val[q]));
241 warning_macro (
"val_min="<<val_min<<
", val_max="<<val_max);
248 template<
class T,
class M>
249 solver_mumps_rep<T,M>::solver_mumps_rep (
const csr<T,M>&
a,
const solver_option& opt)
250 : solver_abstract_rep<
T,
M>(opt),
251 _has_mumps_instance(false),
252 _drop_ext_nnz(false),
260 _drop_ext_nnz = is_distributed<M>::value && opt.force_seq;
263 template<
class T,
class M>
265 solver_mumps_rep<T,M>::update_values (
const csr<T,M>&
a)
268 if (
a.dis_nrow() <= 1) {
271 _a00 = (*(
a.begin()[0])).second;
273 _has_mumps_instance =
false;
276 _has_mumps_instance =
true;
280 bool use_symmetry =
a.is_symmetric();
281 bool use_definite_positive =
a.is_definite_positive();
283 #ifdef _RHEOLEF_HAVE_MPI
284 _mumps_par.comm_fortran = MPI_Comm_c2f (
a.row_ownership().comm());
286 if (use_symmetry && !use_definite_positive) {
290 }
else if (use_symmetry && use_definite_positive) {
301 _mumps_par.cntl[1-1] = 0.0;
314 dmumps_c(&_mumps_par);
315 check_macro (_mumps_par.infog[1-1] == 0,
"mumps: an error has occurred");
319 if (base::option().verbose_level != 0) {
320 _mumps_par.icntl[1-1] = 1;
321 _mumps_par.icntl[2-1] = 1;
322 _mumps_par.icntl[3-1] = 1;
323 _mumps_par.icntl[4-1] = base::option().verbose_level;
326 _mumps_par.icntl[1-1] = -1;
327 _mumps_par.icntl[2-1] = -1;
328 _mumps_par.icntl[3-1] = -1;
329 _mumps_par.icntl[4-1] = 0;
331 if (base::option().compute_determinant) {
332 _mumps_par.icntl[33-1] = 1;
334 _mumps_par.icntl[5-1] = 0;
335 #if defined(_RHEOLEF_HAVE_MUMPS_WITH_PTSCOTCH) || defined(_RHEOLEF_HAVE_MUMPS_WITH_SCOTCH)
336 _mumps_par.icntl[7-1] = 3;
337 #elif defined(_RHEOLEF_HAVE_MUMPS_WITH_PARMETIS) || defined(_RHEOLEF_HAVE_MUMPS_WITH_METIS)
338 _mumps_par.icntl[7-1] = 5;
340 _mumps_par.icntl[7-1] = 7;
342 _mumps_par.icntl[14-1] = 40;
343 _mumps_par.icntl[29-1] = 0;
344 _mumps_par.icntl[22-1] = 0;
345 if (is_distributed<M>::value) {
346 _mumps_par.icntl[18-1] = 3;
348 #if defined(_RHEOLEF_HAVE_MUMPS_WITH_SCOTCH) || defined(_RHEOLEF_HAVE_MUMPS_WITH_METIS)
349 _mumps_par.icntl[28-1] = 1;
350 #elif defined(_RHEOLEF_HAVE_MUMPS_WITH_PTSCOTCH) || defined(_RHEOLEF_HAVE_MUMPS_WITH_PARMETIS)
351 _mumps_par.icntl[28-1] = 2;
354 #ifdef _RHEOLEF_HAVE_MUMPS_WITH_PTSCOTCH
355 _mumps_par.icntl[29-1] = 1;
356 #elif _RHEOLEF_HAVE_MUMPS_WITH_PARMETIS
357 _mumps_par.icntl[29-1] = 2;
359 _mumps_par.icntl[29-1] = 0;
362 _mumps_par.icntl[18-1] = 0;
365 _mumps_par.n =
a.dis_nrow();
366 size_type dia_nnz = ((use_symmetry) ? nnz_dia_upper(
a) :
a.nnz());
367 size_type ext_nnz = ((use_symmetry) ? nnz_ext_upper(
a) :
a.ext_nnz());
368 size_type dis_nnz = dia_nnz + (_drop_ext_nnz ? 0 : ext_nnz);
369 if (is_distributed<M>::value) {
370 #ifdef _RHEOLEF_HAVE_MPI
371 _row.resize (dis_nnz);
372 _col.resize (dis_nnz);
373 csr2mumps_struct_dis (
a, _row, _col, use_symmetry, _drop_ext_nnz);
374 _mumps_par.nz_loc = dis_nnz;
375 _mumps_par.irn_loc = _row.begin().operator->();
376 _mumps_par.jcn_loc = _col.begin().operator->();
379 _row.resize (dia_nnz);
380 _col.resize (dia_nnz);
381 csr2mumps_struct_seq (
a, _row, _col, use_symmetry);
382 _mumps_par.nz = dia_nnz;
383 #if (_RHEOLEF_MUMPS_VERSION_MAJOR >= 5)
384 _mumps_par.nnz = dia_nnz;
386 _mumps_par.irn = _row.begin().operator->();
387 _mumps_par.jcn = _col.begin().operator->();
390 dmumps_c(&_mumps_par);
391 trace_macro (
"symbolic: ordering effectively used = " << _mumps_par.infog[7-1]);
392 check_macro (_mumps_par.infog[1-1] == 0,
"mumps: error infog(1)="<<_mumps_par.infog[1-1]<<
" has occurred (see mumps/infog(1) documentation; infog(2)="<<_mumps_par.infog[2-1]<<
")");
393 if (
a.dis_nrow() <= 1) {
396 _a00 = (*(
a.begin()[0])).second;
403 if (is_distributed<M>::value) {
404 #ifdef _RHEOLEF_HAVE_MPI
405 _val.resize (dis_nnz);
406 csr2mumps_values_dis (
a, _val, use_symmetry, _drop_ext_nnz);
407 _mumps_par.a_loc = _val.begin().operator->();
410 _val.resize (dia_nnz);
411 csr2mumps_values_seq (
a, _val, use_symmetry);
412 _mumps_par.a = _val.begin().operator->();
415 bool finished =
false;
416 while (!finished && _mumps_par.icntl[14-1] <= 2000) {
417 dmumps_c(&_mumps_par);
418 if ((_mumps_par.infog[1-1] == -8 || _mumps_par.infog[1-1] == -9) &&
419 _mumps_par.infog[2-1] >= 0) {
422 _mumps_par.icntl[14-1] += 10;
423 if (_mumps_par.icntl[14-1] > 100) {
424 dis_warning_macro (
"numeric factorization: increase working space of "<<_mumps_par.icntl[14-1]<<
"% and retry...");
431 "solver failed: mumps error infog(1)=" <<_mumps_par.infog[1-1]
432 <<
", infog(2)=" <<_mumps_par.infog[2-1]
433 <<
", icntl(14)="<<_mumps_par.icntl[14-1]
434 <<
", icntl(23)="<<_mumps_par.icntl[23-1]);
435 if (_mumps_par.icntl[33-1] != 0) {
437 _det.mantissa = _mumps_par.rinfog[12-1];
439 _det.exponant = _mumps_par.infog[34-1];
443 template<
class T,
class M>
447 if (
b.dis_size() <= 1) {
450 return vec<T,M>(
b.ownership(),
b[0]/_a00);
452 return vec<T,M>(
b.ownership());
465 if (is_distributed<M>::value) {
467 _mumps_par.icntl[21-1] = 1;
468 _mumps_par.icntl[10-1] = 0;
469 communicator comm =
b.ownership().comm();
470 distributor host_ownership (
b.dis_size(), comm, comm.rank() == 0 ?
b.dis_size() : 0);
471 b_host.resize (host_ownership);
472 size_type first_i =
b.ownership().first_index();
475 b_host.dis_entry(dis_i) =
b[i];
477 b_host.dis_entry_assembly();
478 if (comm.rank() == 0) {
480 _mumps_par.rhs = b_host.begin().operator->();
482 sol_size = _mumps_par.info[23-1];
483 sol.resize (sol_size);
484 perm.resize (sol_size);
485 _mumps_par.lsol_loc = sol_size;
486 _mumps_par.sol_loc = (sol_size > 0) ? sol.begin().operator->() : &dummy;
487 _mumps_par.isol_loc = (sol_size > 0) ? perm.begin().operator->() : &idummy;
490 _mumps_par.icntl[21-1] = 0;
491 _mumps_par.icntl[10-1] = 0;
494 _mumps_par.rhs = x.begin().operator->();
497 _mumps_par.icntl [9-1] = 1;
499 dmumps_c(&_mumps_par);
501 "solver failed: mumps error infog(1)="<<_mumps_par.infog[1-1]
502 <<
", infog(2)="<<_mumps_par.infog[2-1]);
503 if (_mumps_par.infog[1-1] > 0) {
504 warning_macro (
"mumps warning infog(1)="<<_mumps_par.infog[1-1]
505 <<
", infog(2)="<<_mumps_par.infog[2-1]);
507 if (is_distributed<M>::value) {
509 x.resize (
b.ownership());
510 for (
size_t i = 0; i < sol_size; i++) {
512 assert_macro (dis_i < x.dis_size(),
"invalid index");
513 x.dis_entry(dis_i) = sol[i];
515 x.dis_entry_assembly();
519 template<
class T,
class M>
521 solver_mumps_rep<T,M>::trans_solve (
const vec<T,M>&
b)
const
523 error_macro (
"not yet");
526 template<
class T,
class M>
527 solver_mumps_rep<T,M>::~solver_mumps_rep ()
529 if (!_has_mumps_instance)
return;
530 _has_mumps_instance =
false;
535 dmumps_c(&_mumps_par);
536 check_macro (_mumps_par.infog[1-1] == 0,
"mumps: an error has occurred");
543 template class solver_mumps_rep<double,sequential>;
545 #ifdef _RHEOLEF_HAVE_MPI
546 template class solver_mumps_rep<double,distributed>;
field::size_type size_type
#define dis_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 solve(tiny_matrix< T > &a, tiny_vector< size_t > &piv, const tiny_vector< T > &b, tiny_vector< T > &x)