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
65template<
class T,
class M>
67typename csr<T,M>::size_type
68nnz_dia_upper (
const csr<T,M>& a)
70 typedef typename csr<T,M>::size_type
size_type;
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++;
82template<
class T,
class M>
84typename csr<T,M>::size_type
85nnz_ext_upper (
const csr<T,M>& a)
89#ifdef _RHEOLEF_HAVE_MPI
92typename csr<T,distributed>::size_type
93nnz_ext_upper (
const csr<T,distributed>& a)
95 typedef typename csr<T,distributed>::size_type
size_type;
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++;
109template<
class T,
class M>
112csr2mumps_struct_seq (
const csr<T,M>& a, vector<int>& row, vector<int>& col,
bool use_symmetry)
114 typedef typename csr<T,M>::size_type
size_type;
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) {
127template<
class T,
class M>
130csr2mumps_values_seq (
const csr<T,M>& a, vector<T>& val,
bool use_symmetry)
132 typedef typename csr<T,M>::size_type
size_type;
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
148csr2mumps_struct_dis (
const csr<T,sequential>& a, vector<int>& row, vector<int>& col,
bool use_symmetry,
bool)
155csr2mumps_struct_dis (
const csr<T,distributed>& a, vector<int>& row, vector<int>& col,
bool use_symmetry,
bool drop_ext_nnz)
157 typedef typename csr<T,distributed>::size_type
size_type;
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++) {
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) {
193csr2mumps_values_dis (
const csr<T,sequential>& a, vector<T>& val,
bool use_symmetry,
bool)
200csr2mumps_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());
206 typedef typename csr<T,distributed>::size_type
size_type;
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) {
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) {
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]));
248template<
class T,
class M>
251 _has_mumps_instance(false),
252 _drop_ext_nnz(false),
263template<
class T,
class M>
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;
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);
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;
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];
443template<
class T,
class M>
447 if (b.dis_size() <= 1) {
450 return vec<T,M>(b.ownership(), b[0]/_a00);
467 _mumps_par.icntl[21-1] = 1;
468 _mumps_par.icntl[10-1] = 0;
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();
473 for (
size_type i = 0, n = b.size(); i < n; i++) {
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];
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]);
510 for (
size_t i = 0; i < sol_size; i++) {
513 x.dis_entry(dis_i) = sol[i];
515 x.dis_entry_assembly();
519template<
class T,
class M>
526template<
class T,
class M>
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");
545#ifdef _RHEOLEF_HAVE_MPI
field::size_type size_type
see the communicator page for the full documentation
see the csr page for the full documentation
see the distributor page for the full documentation
void update_values(const csr< T, M > &a)
base::size_type size_type
solver_mumps_rep(const csr< T, M > &a, const solver_option &opt=solver_option())
vec< T, M > trans_solve(const vec< T, M > &rhs) const
vec< T, M > solve(const vec< T, M > &rhs) const
see the solver_option page for the full documentation
see the vec page for the full documentation
void resize(const distributor &ownership, const T &init_val=std::numeric_limits< T >::max())
#define trace_macro(message)
#define assert_macro(ok_condition, message)
#define error_macro(message)
#define dis_warning_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.