HermesCommon  3.0
mumps_solver.cpp
Go to the documentation of this file.
1 // This file is part of HermesCommon
2 //
3 // Copyright (c) 2009 hp-FEM group at the University of Nevada, Reno (UNR).
4 // Email: hpfem-group@unr.edu, home page: http://www.hpfem.org/.
5 //
6 // Hermes2D is free software; you can redistribute it and/or modify
7 // it under the terms of the GNU General Public License as published
8 // by the Free Software Foundation; either version 2 of the License,
9 // or (at your option) any later version.
10 //
11 // Hermes2D is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 // GNU General Public License for more details.
15 //
16 // You should have received a copy of the GNU General Public License
17 // along with Hermes2D; if not, write to the Free Software
18 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
22 #include "config.h"
23 #ifdef WITH_MUMPS
24 #include "mumps_solver.h"
25 #include "callstack.h"
26 #include "util/memory_handling.h"
27 
28 namespace Hermes
29 {
30  namespace Algebra
31  {
32  inline double mumps_to_Scalar(double x)
33  {
34  return x;
35  }
36 
37  inline std::complex<double> mumps_to_Scalar(ZMUMPS_COMPLEX x)
38  {
39  return std::complex<double>(x.r, x.i);
40  }
41 
42  double inline real(double x)
43  {
44  return x;
45  }
46 
47  double inline imag(double x)
48  {
49  return 0;
50  }
51 
52  double inline real(ZMUMPS_COMPLEX x)
53  {
54  return mumps_to_Scalar(x).real();
55  }
56 
57  double inline imag(ZMUMPS_COMPLEX x)
58  {
59  return mumps_to_Scalar(x).imag();
60  }
61 
62  extern "C"
63  {
64 #ifndef _WINDOWS
65  extern void dmumps_c(DMUMPS_STRUC_C *mumps_param_ptr);
66  extern void zmumps_c(ZMUMPS_STRUC_C *mumps_param_ptr);
67 #endif
68  }
69 
70 #define USE_COMM_WORLD -987654
71 
72  inline ZMUMPS_COMPLEX& operator +=(ZMUMPS_COMPLEX &a, double b)
73  {
74  a.i += b;
75  return a;
76  }
77 
78  inline ZMUMPS_COMPLEX& operator +=(ZMUMPS_COMPLEX &a, std::complex<double> b)
79  {
80  a.r += b.real();
81  a.i += b.imag();
82  return a;
83  }
84 
85  inline ZMUMPS_COMPLEX& operator +=(ZMUMPS_COMPLEX &a, ZMUMPS_COMPLEX b)
86  {
87  a.r += b.r;
88  a.i += b.i;
89  return a;
90  }
91 
92  inline ZMUMPS_COMPLEX& operator *=(ZMUMPS_COMPLEX &a, std::complex<double> b)
93  {
94  std::complex<double> a_c = std::complex<double>(a.r, a.i);
95  std::complex<double> result = a_c * b;
96  a.r = result.real();
97  a.i = result.imag();
98  return a;
99  }
100 
101  inline void mumps_assign_Scalar(ZMUMPS_COMPLEX & a, std::complex<double> b)
102  {
103  a.r = b.real();
104  a.i = b.imag();
105  }
106 
107  inline void mumps_assign_Scalar(double & a, double b)
108  {
109  a = b;
110  }
111 
112  template<typename Scalar>
113  MumpsMatrix<Scalar>::MumpsMatrix() : CSCMatrix<Scalar>(), irn(nullptr), jcn(nullptr), Ax(nullptr)
114  {
115  }
116 
117  template<typename Scalar>
119  {
120  this->free();
121  }
122 
123  template<typename Scalar>
125  {
126  this->Ax = calloc_with_check<MumpsMatrix<Scalar>, typename mumps_type<Scalar>::mumps_Scalar>(this->nnz, this);
127 
128  irn = malloc_with_check<MumpsMatrix<Scalar>, int>(this->nnz, this);
129  jcn = malloc_with_check<MumpsMatrix<Scalar>, int>(this->nnz, this);
130  for (unsigned int i = 0; i < this->nnz; i++)
131  {
132  irn[i] = 1;
133  jcn[i] = 1;
134  }
135  }
136 
137  template<typename Scalar>
139  {
141  free_with_check(Ax);
142  free_with_check(irn);
143  free_with_check(jcn);
144  }
145 
146  template<typename Scalar>
147  Scalar MumpsMatrix<Scalar>::get(unsigned int m, unsigned int n) const
148  {
149  // Find m-th row in the n-th column.
150  int mid = CSMatrix<Scalar>::find_position(this->Ai + this->Ap[n], this->Ap[n + 1] - this->Ap[n], m);
151  // Return 0 if the entry has not been found.
152  if (mid < 0)
153  return 0.0;
154  else
155  mid += this->Ap[n];
156 
157  return mumps_to_Scalar(Ax[mid]);
158  }
159 
160  template<typename Scalar>
162  {
163  memset(this->Ax, 0, sizeof(typename mumps_type<Scalar>::mumps_Scalar) * this->Ap[this->size]);
164  }
165 
166  template<>
167  void MumpsMatrix<double>::add(unsigned int m, unsigned int n, double v)
168  {
169  // Find m-th row in the n-th column.
170  int pos = CSMatrix<double>::find_position(this->Ai + this->Ap[n], this->Ap[n + 1] - this->Ap[n], m);
171  // Make sure we are adding to an existing non-zero entry.
172  if (pos < 0)
173  throw Hermes::Exceptions::Exception("Sparse matrix entry not found");
174  // Add offset to the n-th column.
175  pos += this->Ap[n];
176 #pragma omp atomic
177  Ax[pos] += v;
178  // MUMPS is indexing from 1
179  irn[pos] = m + 1;
180  jcn[pos] = n + 1;
181  }
182 
183  template<>
184  void MumpsMatrix<std::complex<double> >::add(unsigned int m, unsigned int n, std::complex<double> v)
185  {
186  // Find m-th row in the n-th column.
187  int pos = CSMatrix<std::complex<double> >::find_position(this->Ai + this->Ap[n], this->Ap[n + 1] - this->Ap[n], m);
188  // Make sure we are adding to an existing non-zero entry.
189  if (pos < 0)
190  throw Hermes::Exceptions::Exception("Sparse matrix entry not found");
191  // Add offset to the n-th column.
192  pos += this->Ap[n];
193 #pragma omp critical (MumpsMatrix_add)
194  Ax[pos] += v;
195  // MUMPS is indexing from 1
196  irn[pos] = m + 1;
197  jcn[pos] = n + 1;
198  }
199 
200  template<typename Scalar>
201  void MumpsMatrix<Scalar>::export_to_file(const char *filename, const char *var_name, MatrixExportFormat fmt, char* number_format)
202  {
203  bool invert_storage = false;
204  switch (fmt)
205  {
207  {
208  FILE* file = fopen(filename, "w");
209  if (!file)
210  throw Exceptions::IOException(Exceptions::IOException::Write, filename);
212  fprintf(file, "%%%%MatrixMarket matrix coordinate real general\n");
213  else
214  fprintf(file, "%%%%MatrixMarket matrix coordinate complex general\n");
215 
216  fprintf(file, "%d %d %d\n", this->size, this->size, this->nnz);
217 
218  if (invert_storage)
219  this->switch_orientation();
220  for (unsigned int j = 0; j < this->size; j++)
221  {
222  for (int i = this->Ap[j]; i < this->Ap[j + 1]; i++)
223  {
224  Hermes::Helpers::fprint_coordinate_num(file, this->Ai[i] + 1, j + 1, mumps_to_Scalar(this->Ax[i]), number_format);
225  fprintf(file, "\n");
226  }
227  }
228  if (invert_storage)
229  this->switch_orientation();
230 
231  fclose(file);
232  }
233  break;
234 
236  {
237 #ifdef WITH_MATIO
238  mat_sparse_t sparse;
239  sparse.nzmax = this->nnz;
240  if (invert_storage)
241  this->switch_orientation();
242 
243  sparse.nir = this->nnz;
244  sparse.ir = this->Ai;
245  sparse.njc = this->size + 1;
246  sparse.jc = (int *)this->Ap;
247  sparse.ndata = this->nnz;
248 
249  size_t dims[2];
250  dims[0] = this->size;
251  dims[1] = this->size;
252 
253  mat_t *mat = Mat_CreateVer(filename, "", MAT_FT_MAT5);
254 
255  matvar_t *matvar;
256 
257  // For complex. No allocation here.
258  double* Ax_re = nullptr;
259  double* Ax_im = nullptr;
260 
261  // For real.
263  {
264  sparse.data = this->Ax;
265  matvar = Mat_VarCreate(var_name, MAT_C_SPARSE, MAT_T_DOUBLE, 2, dims, &sparse, MAT_F_DONT_COPY_DATA);
266  }
267  else
268  {
269  // For complex.
270  Ax_re = malloc_with_check<CSMatrix<Scalar>, double>(this->nnz, this);
271  Ax_im = malloc_with_check<CSMatrix<Scalar>, double>(this->nnz, this);
272  struct mat_complex_split_t z = { Ax_re, Ax_im };
273 
274  for (int i = 0; i < this->nnz; i++)
275  {
276  Ax_re[i] = ((std::complex<double>)(mumps_to_Scalar(this->Ax[i]))).real();
277  Ax_im[i] = ((std::complex<double>)(mumps_to_Scalar(this->Ax[i]))).imag();
278  sparse.data = &z;
279  }
280  matvar = Mat_VarCreate(var_name, MAT_C_SPARSE, MAT_T_DOUBLE, 2, dims, &sparse, MAT_F_DONT_COPY_DATA | MAT_F_COMPLEX);
281  }
282 
283  if (matvar)
284  {
285  Mat_VarWrite(mat, matvar, MAT_COMPRESSION_ZLIB);
286  Mat_VarFree(matvar);
287  }
288  if (invert_storage)
289  this->switch_orientation();
290  free_with_check(Ax_re);
291  free_with_check(Ax_im);
292  Mat_Close(mat);
293 
294  if (!matvar)
295  throw Exceptions::IOException(Exceptions::IOException::Write, filename);
296 #endif
297  }
298  break;
299 
301  {
302  FILE* file = fopen(filename, "w");
303  if (!file)
304  throw Exceptions::IOException(Exceptions::IOException::Write, filename);
305 
306  if (invert_storage)
307  this->switch_orientation();
308  for (unsigned int j = 0; j < this->size; j++)
309  {
310  for (int i = this->Ap[j]; i < this->Ap[j + 1]; i++)
311  {
312  Helpers::fprint_coordinate_num(file, this->Ai[i], j, mumps_to_Scalar(this->Ax[i]), number_format);
313  fprintf(file, "\n");
314  }
315  }
316  if (invert_storage)
317  this->switch_orientation();
318 
319  fclose(file);
320  }
321  break;
322 
323 #ifdef WITH_BSON
324  case EXPORT_FORMAT_BSON:
325  {
326  // Init bson
327  bson bw;
328  bson_init(&bw);
329 
330  // Matrix size.
331  bson_append_int(&bw, "size", this->size);
332  // Nonzeros.
333  bson_append_int(&bw, "nnz", this->nnz);
334 
335  if (invert_storage)
336  this->switch_orientation();
337 
338  bson_append_start_array(&bw, "Ap");
339  for (unsigned int i = 0; i < this->size; i++)
340  bson_append_int(&bw, "p", this->Ap[i]);
341  bson_append_finish_array(&bw);
342 
343  bson_append_start_array(&bw, "Ai");
344  for (unsigned int i = 0; i < this->nnz; i++)
345  bson_append_int(&bw, "i", this->Ai[i]);
346  bson_append_finish_array(&bw);
347 
348  bson_append_start_array(&bw, "Ax");
349  for (unsigned int i = 0; i < this->nnz; i++)
350  bson_append_double(&bw, "x", real(this->Ax[i]));
351  bson_append_finish_array(&bw);
352 
354  {
355  bson_append_start_array(&bw, "Ax-imag");
356  for (unsigned int i = 0; i < this->nnz; i++)
357  bson_append_double(&bw, "x-i", imag(this->Ax[i]));
358  bson_append_finish_array(&bw);
359  }
360  bson_append_finish_array(&bw);
361 
362  if (invert_storage)
363  this->switch_orientation();
364 
365  // Done.
366  bson_finish(&bw);
367 
368  // Write to disk.
369  FILE *fpw;
370  fpw = fopen(filename, "wb");
371  const char *dataw = (const char *)bson_data(&bw);
372  fwrite(dataw, bson_size(&bw), 1, fpw);
373  fclose(fpw);
374 
375  bson_destroy(&bw);
376  }
377  break;
378 #endif
379  }
380  }
381 
382  template<typename Scalar>
383  void MumpsMatrix<Scalar>::import_from_file(const char *filename, const char *var_name, MatrixExportFormat fmt)
384  {
385  bool invert_storage = false;
386  switch (fmt)
387  {
389  throw Exceptions::MethodNotOverridenException("CSMatrix<Scalar>::import_from_file - Matrix Market");
390  break;
392  {
393 #ifdef WITH_MATIO
394  mat_t *matfp;
395  matvar_t *matvar;
396 
397  matfp = Mat_Open(filename, MAT_ACC_RDONLY);
398 
399  if (!matfp)
400  throw Exceptions::IOException(Exceptions::IOException::Read, filename);
401 
402  matvar = Mat_VarRead(matfp, var_name);
403 
404  if (matvar)
405  {
406  mat_sparse_t *sparse = (mat_sparse_t *)matvar->data;
407 
408  this->nnz = sparse->nir;
409  this->size = sparse->njc;
410 
411  this->Ap = malloc_with_check<MumpsMatrix<Scalar>, int>(this->size + 1, this);
412  this->Ai = malloc_with_check<MumpsMatrix<Scalar>, int>(this->nnz, this);
413  this->Ax = malloc_with_check<MumpsMatrix<Scalar>, typename mumps_type<Scalar>::mumps_Scalar>(this->nnz, this);
414  this->irn = malloc_with_check<MumpsMatrix<Scalar>, int>(this->nnz, this);
415  this->jcn = malloc_with_check<MumpsMatrix<Scalar>, int>(this->nnz, this);
416 
417  void* data = nullptr;
419  data = sparse->data;
420  else
421  {
422  std::complex<double>* complex_data = malloc_with_check<CSMatrix<Scalar>, std::complex<double> >(this->nnz, this);
423  double* real_array = (double*)((mat_complex_split_t*)sparse->data)->Re;
424  double* imag_array = (double*)((mat_complex_split_t*)sparse->data)->Im;
425  for (int i = 0; i < this->nnz; i++)
426  complex_data[i] = std::complex<double>(real_array[i], imag_array[i]);
427  data = (void*)complex_data;
428  }
429  memcpy(this->Ax, data, this->nnz * sizeof(Scalar));
431  free_with_check(data);
432  memcpy(this->Ap, sparse->jc, this->size * sizeof(int));
433  this->Ap[this->size] = this->nnz;
434  memcpy(this->Ai, sparse->ir, this->nnz * sizeof(int));
435 
436  for (unsigned int i = 0; i < this->nnz; i++)
437  this->irn[i] = this->Ai[i] + 1;
438  for (unsigned int i = 0; i < this->size; i++)
439  {
440  for (int j = this->Ap[i]; j < this->Ap[i + 1]; j++)
441  jcn[j] = i + 1;
442  }
443 
444  if (invert_storage)
445  this->switch_orientation();
446  }
447 
448  Mat_Close(matfp);
449 
450  if (!matvar)
451  throw Exceptions::IOException(Exceptions::IOException::Read, filename);
452 #else
453  throw Exceptions::Exception("MATIO not included.");
454 #endif
455  }
456  break;
457 
459  throw Exceptions::MethodNotOverridenException("CSMatrix<Scalar>::import_from_file - Simple format");
460  break;
461 #ifdef WITH_BSON
462  case EXPORT_FORMAT_BSON:
463  {
464  FILE *fpr;
465  fpr = fopen(filename, "rb");
466 
467  // file size:
468  fseek(fpr, 0, SEEK_END);
469  int size = ftell(fpr);
470  rewind(fpr);
471 
472  // allocate memory to contain the whole file:
473  char *datar = malloc_with_check<char>(size);
474  fread(datar, size, 1, fpr);
475  fclose(fpr);
476 
477  bson br;
478  bson_init_finished_data(&br, datar, 0);
479 
480  bson_iterator it;
481  bson sub;
482  bson_find(&it, &br, "size");
483  this->size = bson_iterator_int(&it);
484  bson_find(&it, &br, "nnz");
485  this->nnz = bson_iterator_int(&it);
486 
487  this->Ap = malloc_with_check<MumpsMatrix<Scalar>, int>(this->size + 1, this);
488  this->Ai = malloc_with_check<MumpsMatrix<Scalar>, int>(this->nnz, this);
489  this->Ax = malloc_with_check<MumpsMatrix<Scalar>, typename mumps_type<Scalar>::mumps_Scalar>(this->nnz, this);
490  this->irn = malloc_with_check<MumpsMatrix<Scalar>, int>(this->nnz, this);
491  this->jcn = malloc_with_check<MumpsMatrix<Scalar>, int>(this->nnz, this);
492 
493  // coeffs
494  bson_iterator it_coeffs;
495  bson_find(&it_coeffs, &br, "Ap");
496  bson_iterator_subobject_init(&it_coeffs, &sub, 0);
497  bson_iterator_init(&it, &sub);
498  int index_coeff = 0;
499  while (bson_iterator_next(&it))
500  this->Ap[index_coeff++] = bson_iterator_int(&it);
501  this->Ap[this->size] = this->nnz;
502 
503  bson_find(&it_coeffs, &br, "Ai");
504  bson_iterator_subobject_init(&it_coeffs, &sub, 0);
505  bson_iterator_init(&it, &sub);
506  index_coeff = 0;
507  while (bson_iterator_next(&it))
508  this->Ai[index_coeff++] = bson_iterator_int(&it);
509 
510  bson_find(&it_coeffs, &br, "Ax");
511  bson_iterator_subobject_init(&it_coeffs, &sub, 0);
512  bson_iterator_init(&it, &sub);
513  index_coeff = 0;
514  while (bson_iterator_next(&it))
515  mumps_assign_Scalar(this->Ax[index_coeff++], bson_iterator_double(&it));
516 
518  {
519  bson_find(&it_coeffs, &br, "Ax-imag");
520  bson_iterator_subobject_init(&it_coeffs, &sub, 0);
521  bson_iterator_init(&it, &sub);
522  index_coeff = 0;
523  while (bson_iterator_next(&it))
524  this->Ax[index_coeff++] += bson_iterator_double(&it);
525  }
526 
527  for (unsigned int i = 0; i < this->nnz; i++)
528  this->irn[i] = this->Ai[i] + 1;
529  for (unsigned int i = 0; i < this->size; i++)
530  {
531  for (int j = this->Ap[i]; j < this->Ap[i + 1]; j++)
532  jcn[j] = i + 1;
533  }
534 
535  if (invert_storage)
536  this->switch_orientation();
537 
538  bson_destroy(&br);
539  free_with_check(datar);
540  }
541 #endif
542  }
543  }
544 
545  template<typename Scalar>
546  void MumpsMatrix<Scalar>::multiply_with_vector(Scalar* vector_in, Scalar*& vector_out, bool vector_out_initialized) const
547  {
548  if (!vector_out_initialized)
549  vector_out = malloc_with_check<Scalar>(this->size);
550 
551  memset(vector_out, 0, sizeof(Scalar)* this->size);
552 
553  Scalar a;
554  for (unsigned int i = 0; i < this->nnz; i++)
555  {
556  a = mumps_to_Scalar(Ax[i]);
557  vector_out[jcn[i] - 1] += vector_in[irn[i] - 1] * a;
558  }
559  }
560 
561  template<typename Scalar>
563  {
564  for (int i = 0; i < this->nnz; i++)
565  {
566  Ax[i] *= value;
567  }
568  }
569 
570  template<typename Scalar>
571  void MumpsMatrix<Scalar>::create(unsigned int size, unsigned int nnz_, int* ap, int* ai, Scalar* ax)
572  {
573  this->nnz = nnz_;
574  this->size = size;
575  this->Ap = malloc_with_check<MumpsMatrix<Scalar>, int>(this->size + 1, this);
576  this->Ai = malloc_with_check<MumpsMatrix<Scalar>, int>(this->nnz, this);
577  this->Ax = malloc_with_check<MumpsMatrix<Scalar>, typename mumps_type<Scalar>::mumps_Scalar>(this->nnz, this);
578  irn = malloc_with_check<MumpsMatrix<Scalar>, int>(this->nnz, this);
579  jcn = malloc_with_check<MumpsMatrix<Scalar>, int>(this->nnz, this);
580 
581  for (unsigned int i = 0; i < this->size; i++)
582  {
583  this->Ap[i] = ap[i];
584  for (int j = ap[i]; j < ap[i + 1]; j++)
585  jcn[j] = i;
586  }
587  this->Ap[this->size] = ap[this->size];
588  for (unsigned int i = 0; i < this->nnz; i++)
589  {
590  mumps_assign_Scalar(this->Ax[i], ax[i]);
591  this->Ai[i] = ai[i];
592  irn[i] = ai[i];
593  }
594  }
595 
596  template<typename Scalar>
598  {
600 
601  nmat->nnz = this->nnz;
602  nmat->size = this->size;
603  nmat->Ap = malloc_with_check<MumpsMatrix<Scalar>, int>(this->size + 1, nmat);
604  nmat->Ai = malloc_with_check<MumpsMatrix<Scalar>, int>(this->nnz, nmat);
605  nmat->Ax = malloc_with_check<MumpsMatrix<Scalar>, typename mumps_type<Scalar>::mumps_Scalar>(this->nnz, nmat);
606  nmat->irn = malloc_with_check<MumpsMatrix<Scalar>, int>(this->nnz, nmat);
607  nmat->jcn = malloc_with_check<MumpsMatrix<Scalar>, int>(this->nnz, nmat);
608  for (unsigned int i = 0; i < this->nnz; i++)
609  {
610  nmat->Ai[i] = this->Ai[i];
611  nmat->Ax[i] = Ax[i];
612  nmat->irn[i] = irn[i];
613  nmat->jcn[i] = jcn[i];
614  }
615  for (unsigned int i = 0; i < this->size + 1; i++)
616  {
617  nmat->Ap[i] = this->Ap[i];
618  }
619  return nmat;
620  }
621 
622  template class HERMES_API MumpsMatrix < double > ;
623  template class HERMES_API MumpsMatrix < std::complex<double> > ;
624  }
625 
626  namespace Solvers
627  {
629 #define ICNTL(I) icntl[(I)-1]
630 #define MUMPS_INFO(param, I) (param).infog[(I)-1]
631 #define INFOG(I) infog[(I)-1]
632 
634 #define JOB_INIT -1
635 #define JOB_END -2
636 #define JOB_ANALYZE_FACTORIZE_SOLVE 6
637 #define JOB_FACTORIZE_SOLVE 5
638 #define JOB_SOLVE 3
639 
640  template<typename Scalar>
642  DirectSolver<Scalar>(m, rhs), m(m), rhs(rhs), icntl_14(init_icntl_14)
643  {
644  inited = false;
645 
646  // Initial values for some fields of the MUMPS_STRUC structure that may be accessed
647  // before MUMPS has been initialized.
648  param.rhs = nullptr;
649  // see the case HERMES_REUSE_MATRIX_REORDERING_AND_SCALING
650  param.INFOG(33) = -999;
651  // in setup_factorization()
652  }
653 
654  template<typename Scalar>
656  {
657  free();
658  }
659 
660  template<typename Scalar>
662  {
663  // Terminate the current instance of MUMPS.
664  if (inited)
665  {
666  param.job = JOB_END;
667  mumps_c(&param);
668  }
669 
670  if (param.rhs != nullptr)
671  free_with_check(param.rhs);
672  }
673 
674  template<>
676  {
677  dmumps_c(param);
678  }
679 
680  template<>
681  void MumpsSolver<std::complex<double> >::mumps_c(mumps_type<std::complex<double> >::mumps_struct * param)
682  {
683  zmumps_c(param);
684  }
685 
686  template<typename Scalar>
688  {
689  switch (param.INFOG(1))
690  {
691  // no error
692  case 0: return true;
693  case -1: throw Hermes::Exceptions::LinearMatrixSolverException("Error occured on processor %d", MUMPS_INFO(param, 2)); break;
694  case -2: throw Hermes::Exceptions::LinearMatrixSolverException("Number of nonzeros (NNZ) is out of range."); break;
695  case -3: throw Hermes::Exceptions::LinearMatrixSolverException("MUMPS called with an invalid option for JOB."); break;
696  case -5: throw Hermes::Exceptions::LinearMatrixSolverException("Problem of REAL or COMPLEX workspace allocation of size %i during analysis.", param.INFOG(2)); break;
697  case -6: throw Hermes::Exceptions::LinearMatrixSolverException("Matrix is singular in structure."); break;
698  case -7: throw Hermes::Exceptions::LinearMatrixSolverException("Problem of INTEGER workspace allocation of size %i during analysis.", param.INFOG(2)); break;
699  case -8:
700  case -9: return false;
701  case -10: throw Hermes::Exceptions::LinearMatrixSolverException("Numerically singular matrix."); break;
702  default: Hermes::Exceptions::LinearMatrixSolverException("Non-detailed exception in MUMPS: INFOG(1) = %d", param.INFOG(1)); break;
703  }
704  return false;
705  }
706 
707  template<typename Scalar>
709  {
710  if (inited)
711  {
712  // If there is already an instance of MUMPS running,
713  // terminate it.
714  param.job = JOB_END;
715  mumps_c(&param);
716  }
717 
718  param.job = JOB_INIT;
719  // host also performs calculations
720  param.par = 1;
721  // 0 = unsymmetric
722  param.sym = 0;
723  param.comm_fortran = USE_COMM_WORLD;
724 
725  mumps_c(&param);
726  inited = check_status();
727 
728  if (inited)
729  {
730  // No printings.
731  param.ICNTL(1) = -1;
732  param.ICNTL(2) = -1;
733  param.ICNTL(3) = -1;
734  param.ICNTL(4) = 0;
735 
736  // =/ both centralized assembled matrix
737  param.ICNTL(5) = 0;
738  // =\ both centralized assembled matrix
739  param.ICNTL(18) = 0;
740  // centralized dense RHS
741  param.ICNTL(20) = 0;
742  // centralized dense solution
743  param.ICNTL(21) = 0;
744 
745  // Fixing the memory problems - this parameter specifies the maximum
746  // extra fill-in.
747  // Extract from the docs (4.10, page 27):
748  // ICNTL(14) is accessed by the host both during the analysis and the factorization phases. It corresponds
749  // to the percentage increase in the estimated working space. When significant extra fill-in is caused
750  // by numerical pivoting, increasing ICNTL(14) may help. Except in special cases, the default value
751  // is 20 (which corresponds to a 20 % increase).
752  param.ICNTL(14) = 100 * this->icntl_14;
753 
754  // Specify the matrix.
755  param.n = m->size;
756  param.nz = m->nnz;
757  param.irn = m->irn;
758  param.jcn = m->jcn;
759  param.a = m->Ax;
760  }
761 
762  return inited;
763  }
764 
765  template<typename Scalar>
767  {
768  return m->size;
769  }
770 
771  template<typename Scalar>
773  {
774  assert(m != nullptr);
775  assert(rhs != nullptr);
776 
777  this->tick();
778 
779  // Prepare the MUMPS data structure with input for the solver driver
780  // (according to the chosen factorization reuse strategy), as well as
781  // the system matrix.
782  if (!setup_factorization())
783  throw Hermes::Exceptions::LinearMatrixSolverException("LU factorization could not be completed.");
784 
785  // Specify the right-hand side (will be replaced by the solution).
786  param.rhs = malloc_with_check<MumpsSolver<Scalar>, typename mumps_type<Scalar>::mumps_Scalar>(m->size, this);
787  memcpy(param.rhs, rhs->v, m->size * sizeof(typename mumps_type<Scalar>::mumps_Scalar));
788 
789  // Do the jobs specified in setup_factorization().
790  mumps_c(&param);
791 
792  // Throws appropriate exception.
793  if (check_status())
794  {
795  free_with_check(this->sln);
796  this->sln = malloc_with_check<MumpsSolver<Scalar>, Scalar>(m->size, this);
797  for (unsigned int i = 0; i < rhs->get_size(); i++)
798  this->sln[i] = mumps_to_Scalar(param.rhs[i]);
799  }
800  else
801  {
802  free_with_check(param.rhs);
803 
804  icntl_14 *= 2;
805  if (icntl_14 > max_icntl_14)
806  throw Hermes::Exceptions::LinearMatrixSolverException("MUMPS memory overflow - potentially singular matrix");
807  else
808  {
809  this->reinit();
810  this->solve();
811  }
812  return;
813  /* From the MUMPS docs - these two cases are forwarded from check_status().
814  case 8: throw Hermes::Exceptions::LinearMatrixSolverException("Main internal integer workarray IS too small for factorization. This may happen, for example, if
815  numerical pivoting leads to significantly more fill-in than was predicted by the analysis. The user
816  should increase the value of ICNTL(14) before calling the factorization again (JOB=2).
817  case 9: Main internal real/complex workarray S too small. If INFO(2) is positive, then the number of entries
818  that are missing in S at the moment when the error is raised is available in INFO(2). If INFO(2) is
819  negative, then its absolute value should be multiplied by 1 million. If an error 9 occurs, the user
820  should increase the value of ICNTL(14) before calling the factorization (JOB=2) again, except if
821  ICNTL(23) is provided, in which case ICNTL(23) should be increased.
822  */
823  }
824 
825  this->tick();
826  this->time = this->accumulated();
827 
828  free_with_check(param.rhs);
829  param.rhs = nullptr;
830  }
831 
832  template<typename Scalar>
834  {
835  // When called for the first time, all three phases (analysis, factorization,
836  // solution) must be performed.
837  int eff_fact_scheme = this->reuse_scheme;
838  if (!inited)
839  if (this->reuse_scheme == HERMES_REUSE_MATRIX_REORDERING || this->reuse_scheme == HERMES_REUSE_MATRIX_STRUCTURE_COMPLETELY)
840  eff_fact_scheme = HERMES_CREATE_STRUCTURE_FROM_SCRATCH;
841 
842  switch (eff_fact_scheme)
843  {
845  // (Re)initialize new_ instance.
846  reinit();
847 
848  // Let MUMPS decide when and how to compute matrix reordering and scaling.
849  param.ICNTL(6) = 7;
850  param.ICNTL(8) = 77;
851  param.job = JOB_ANALYZE_FACTORIZE_SOLVE;
852 
853  break;
855  // Let MUMPS reuse results of the symbolic analysis and perform
856  // scaling during each factorization (values 1-8 may be set here,
857  // corresponding to different scaling algorithms during factorization;
858  // see the MUMPS documentation for details).
859  param.ICNTL(8) = 7;
860  param.job = JOB_FACTORIZE_SOLVE;
861 
862  break;
864  // Perform scaling along with reordering during the symbolic analysis phase
865  // and then reuse it during subsequent factorizations. new_ instance of MUMPS
866  // has to be created before the analysis phase.
867  if (param.INFOG(33) != -2)
868  {
869  reinit();
870  param.ICNTL(6) = 5;
871  param.job = JOB_ANALYZE_FACTORIZE_SOLVE;
872  // After analysis is done, INFOG(33) will be set to -2 by MUMPS.
873  }
874  else
875  {
876  param.job = JOB_FACTORIZE_SOLVE;
877  }
878  break;
880  param.job = JOB_SOLVE;
881  break;
882  }
883 
884  return true;
885  }
886 
887  template class HERMES_API MumpsSolver < double > ;
888  template class HERMES_API MumpsSolver < std::complex<double> > ;
889  }
890 }
891 #endif
unsigned int nnz
Number of non-zero entries ( = Ap[size]).
Definition: cs_matrix.h:123
General namespace for the Hermes library.
Exception interface Basically a std::exception, but with a constructor with string and with print_msg...
Definition: exceptions.h:49
unsigned int size
matrix size
Definition: matrix.h:97
int * Ap
Index to Ax/Ai, where each column / row starts.
Definition: cs_matrix.h:121
File containing common definitions, and basic global enums etc. for HermesCommon. ...
MatrixExportFormat
Format of file matrix and vector output.
Plain ascii file lines contains row column and value.
Vector used with MUMPS solver.
Definition: vector.h:116
int * jcn
Column indices.
Definition: mumps_solver.h:116
General CSC Matrix class. (can be used in umfpack, in that case use the CSCMatrix subclass...
Definition: cs_matrix.h:131
IO exception. Internal. Exception occurs when something fails to be written to / read from the disk...
Definition: exceptions.h:76
mumps_type< Scalar >::mumps_struct param
MUMPS specific structure with solver parameters.
Definition: mumps_solver.h:152
Matrix used with MUMPS solver.
Definition: mumps_solver.h:77
#define JOB_INIT
Job definitions according to MUMPS documentation.
File containing functionality for investigating call stack.
int * Ai
Row / Column indices of values in Ax.
Definition: cs_matrix.h:119
Type deductors.
Definition: common.h:276
MUMPS solver interface.
int * irn
Row indices.
Definition: mumps_solver.h:114
pattern as the one already factorized.
Method is not overriden and should be.
Definition: exceptions.h:201
General CS Matrix class. Either row- or column- specific (see subclassses).
Definition: cs_matrix.h:35
mumps_type< Scalar >::mumps_Scalar * Ax
Matrix entries (column-wise).
Definition: mumps_solver.h:118
virtual int get_matrix_size()
Get size of matrix.
Matrix Market which can be read by pysparse library.