GetFEM  5.4.2
gmm_MUMPS_interface.h
Go to the documentation of this file.
1 /* -*- c++ -*- (enables emacs c++ mode) */
2 /*===========================================================================
3 
4  Copyright (C) 2003-2020 Yves Renard, Julien Pommier
5 
6  This file is a part of GetFEM
7 
8  GetFEM is free software; you can redistribute it and/or modify it
9  under the terms of the GNU Lesser General Public License as published
10  by the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version along with the GCC Runtime Library
12  Exception either version 3.1 or (at your option) any later version.
13  This program is distributed in the hope that it will be useful, but
14  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16  License and GCC Runtime Library Exception for more details.
17  You should have received a copy of the GNU Lesser General Public License
18  along with this program; if not, write to the Free Software Foundation,
19  Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
20 
21  As a special exception, you may use this file as it is a part of a free
22  software library without restriction. Specifically, if other files
23  instantiate templates or use macros or inline functions from this file,
24  or you compile this file and link it with other files to produce an
25  executable, this file does not by itself cause the resulting executable
26  to be covered by the GNU Lesser General Public License. This exception
27  does not however invalidate any other reasons why the executable file
28  might be covered by the GNU Lesser General Public License.
29 
30 ===========================================================================*/
31 
32 /**@file gmm_MUMPS_interface.h
33  @author Yves Renard <Yves.Renard@insa-lyon.fr>,
34  @author Julien Pommier <Julien.Pommier@insa-toulouse.fr>
35  @date December 8, 2005.
36  @brief Interface with MUMPS (LU direct solver for sparse matrices).
37 */
38 #if defined(GMM_USES_MUMPS) || defined(HAVE_DMUMPS_C_H)
39 
40 #ifndef GMM_MUMPS_INTERFACE_H
41 #define GMM_MUMPS_INTERFACE_H
42 
43 #include "gmm_kernel.h"
44 
45 
46 extern "C" {
47 
48 #include <smumps_c.h>
49 #undef F_INT
50 #undef F_DOUBLE
51 #undef F_DOUBLE2
52 #include <dmumps_c.h>
53 #undef F_INT
54 #undef F_DOUBLE
55 #undef F_DOUBLE2
56 #include <cmumps_c.h>
57 #undef F_INT
58 #undef F_DOUBLE
59 #undef F_DOUBLE2
60 #include <zmumps_c.h>
61 #undef F_INT
62 #undef F_DOUBLE
63 #undef F_DOUBLE2
64 
65 }
66 
67 namespace gmm {
68 
69 #define ICNTL(I) icntl[(I)-1]
70 #define INFO(I) info[(I)-1]
71 #define INFOG(I) infog[(I)-1]
72 #define RINFOG(I) rinfog[(I)-1]
73 
74  template <typename T> struct ij_sparse_matrix {
75  std::vector<int> irn;
76  std::vector<int> jcn;
77  std::vector<T> a;
78  bool sym;
79 
80  template <typename L> void store(const L& l, size_type i) {
81  typename linalg_traits<L>::const_iterator it = vect_const_begin(l),
82  ite = vect_const_end(l);
83  for (; it != ite; ++it) {
84  int ir = (int)i + 1, jc = (int)it.index() + 1;
85  if (*it != T(0) && (!sym || ir >= jc))
86  { irn.push_back(ir); jcn.push_back(jc); a.push_back(*it); }
87  }
88  }
89 
90  template <typename L> void build_from(const L& l, row_major) {
91  for (size_type i = 0; i < mat_nrows(l); ++i)
92  store(mat_const_row(l, i), i);
93  }
94 
95  template <typename L> void build_from(const L& l, col_major) {
96  for (size_type i = 0; i < mat_ncols(l); ++i)
97  store(mat_const_col(l, i), i);
98  irn.swap(jcn);
99  }
100 
101  template <typename L> ij_sparse_matrix(const L& A, bool sym_) {
102  size_type nz = nnz(A);
103  sym = sym_;
104  irn.reserve(nz); jcn.reserve(nz); a.reserve(nz);
105  build_from(A, typename principal_orientation_type<typename
106  linalg_traits<L>::sub_orientation>::potype());
107  }
108  };
109 
110  /* ********************************************************************* */
111  /* MUMPS solve interface */
112  /* ********************************************************************* */
113 
114  template <typename T> struct mumps_interf {};
115 
116  template <> struct mumps_interf<float> {
117  typedef SMUMPS_STRUC_C MUMPS_STRUC_C;
118  typedef float value_type;
119 
120  static void mumps_c(MUMPS_STRUC_C &id) { smumps_c(&id); }
121  };
122 
123  template <> struct mumps_interf<double> {
124  typedef DMUMPS_STRUC_C MUMPS_STRUC_C;
125  typedef double value_type;
126  static void mumps_c(MUMPS_STRUC_C &id) { dmumps_c(&id); }
127  };
128 
129  template <> struct mumps_interf<std::complex<float> > {
130  typedef CMUMPS_STRUC_C MUMPS_STRUC_C;
131  typedef mumps_complex value_type;
132  static void mumps_c(MUMPS_STRUC_C &id) { cmumps_c(&id); }
133  };
134 
135  template <> struct mumps_interf<std::complex<double> > {
136  typedef ZMUMPS_STRUC_C MUMPS_STRUC_C;
137  typedef mumps_double_complex value_type;
138  static void mumps_c(MUMPS_STRUC_C &id) { zmumps_c(&id); }
139  };
140 
141 
142  template <typename MUMPS_STRUCT>
143  static inline bool mumps_error_check(MUMPS_STRUCT &id) {
144  if (id.INFO(1) < 0) {
145  switch (id.INFO(1)) {
146  case -2:
147  GMM_ASSERT1(false, "Solve with MUMPS failed: NZ = " << id.INFO(2)
148  << " is out of range");
149  case -6 : case -10 :
150  GMM_WARNING1("Solve with MUMPS failed: matrix is singular");
151  return false;
152  case -9:
153  GMM_ASSERT1(false, "Solve with MUMPS failed: error "
154  << id.INFO(1) << ", increase ICNTL(14)");
155  case -13 :
156  GMM_ASSERT1(false, "Solve with MUMPS failed: not enough memory");
157  default :
158  GMM_ASSERT1(false, "Solve with MUMPS failed with error "
159  << id.INFO(1));
160  }
161  }
162  return true;
163  }
164 
165 
166  /** MUMPS solve interface
167  * Works only with sparse or skyline matrices
168  */
169  template <typename MAT, typename VECTX, typename VECTB>
170  bool MUMPS_solve(const MAT &A, const VECTX &X_, const VECTB &B,
171  bool sym = false, bool distributed = false) {
172  VECTX &X = const_cast<VECTX &>(X_);
173 
174  typedef typename linalg_traits<MAT>::value_type T;
175  typedef typename mumps_interf<T>::value_type MUMPS_T;
176  GMM_ASSERT2(gmm::mat_nrows(A) == gmm::mat_ncols(A), "Non-square matrix");
177 
178  std::vector<T> rhs(gmm::vect_size(B)); gmm::copy(B, rhs);
179 
180  ij_sparse_matrix<T> AA(A, sym);
181 
182  const int JOB_INIT = -1;
183  const int JOB_END = -2;
184  const int USE_COMM_WORLD = -987654;
185 
186  typename mumps_interf<T>::MUMPS_STRUC_C id;
187 
188  int rank(0);
189 #ifdef GMM_USES_MPI
190  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
191 #endif
192 
193  id.job = JOB_INIT;
194  id.par = 1;
195  id.sym = sym ? 2 : 0;
196  id.comm_fortran = USE_COMM_WORLD;
197  mumps_interf<T>::mumps_c(id);
198 
199  if (rank == 0 || distributed) {
200  id.n = int(gmm::mat_nrows(A));
201  if (distributed) {
202  id.nz_loc = int(AA.irn.size());
203  id.irn_loc = &(AA.irn[0]);
204  id.jcn_loc = &(AA.jcn[0]);
205  id.a_loc = (MUMPS_T*)(&(AA.a[0]));
206  } else {
207  id.nz = int(AA.irn.size());
208  id.irn = &(AA.irn[0]);
209  id.jcn = &(AA.jcn[0]);
210  id.a = (MUMPS_T*)(&(AA.a[0]));
211  }
212  if (rank == 0)
213  id.rhs = (MUMPS_T*)(&(rhs[0]));
214  }
215 
216  id.ICNTL(1) = -1; // output stream for error messages
217  id.ICNTL(2) = -1; // output stream for other messages
218  id.ICNTL(3) = -1; // output stream for global information
219  id.ICNTL(4) = 0; // verbosity level
220 
221  if (distributed)
222  id.ICNTL(5) = 0; // assembled input matrix (default)
223 
224  id.ICNTL(14) += 80; /* small boost to the workspace size as we have encountered some problem
225  who did not fit in the default settings of mumps..
226  by default, ICNTL(14) = 15 or 20
227  */
228  //cout << "ICNTL(14): " << id.ICNTL(14) << "\n";
229 
230  if (distributed)
231  id.ICNTL(18) = 3; // strategy for distributed input matrix
232 
233  // id.ICNTL(22) = 1; /* enables out-of-core support */
234 
235  id.job = 6;
236  mumps_interf<T>::mumps_c(id);
237  bool ok = mumps_error_check(id);
238 
239  id.job = JOB_END;
240  mumps_interf<T>::mumps_c(id);
241 
242 #ifdef GMM_USES_MPI
243  MPI_Bcast(&(rhs[0]),id.n,gmm::mpi_type(T()),0,MPI_COMM_WORLD);
244 #endif
245 
246  gmm::copy(rhs, X);
247 
248  return ok;
249 
250  }
251 
252 
253 
254  /** MUMPS solve interface for distributed matrices
255  * Works only with sparse or skyline matrices
256  */
257  template <typename MAT, typename VECTX, typename VECTB>
258  bool MUMPS_distributed_matrix_solve(const MAT &A, const VECTX &X_,
259  const VECTB &B, bool sym = false) {
260  return MUMPS_solve(A, X_, B, sym, true);
261  }
262 
263 
264 
265  template<typename T>
266  inline T real_or_complex(std::complex<T> a) { return a.real(); }
267  template<typename T>
268  inline T real_or_complex(T &a) { return a; }
269 
270 
271  /** Evaluate matrix determinant with MUMPS
272  * Works only with sparse or skyline matrices
273  */
274  template <typename MAT, typename T = typename linalg_traits<MAT>::value_type>
275  T MUMPS_determinant(const MAT &A, int &exponent,
276  bool sym = false, bool distributed = false) {
277  exponent = 0;
278  typedef typename mumps_interf<T>::value_type MUMPS_T;
279  typedef typename number_traits<T>::magnitude_type R;
280  GMM_ASSERT2(gmm::mat_nrows(A) == gmm::mat_ncols(A), "Non-square matrix");
281 
282  ij_sparse_matrix<T> AA(A, sym);
283 
284  const int JOB_INIT = -1;
285  const int JOB_END = -2;
286  const int USE_COMM_WORLD = -987654;
287 
288  typename mumps_interf<T>::MUMPS_STRUC_C id;
289 
290  int rank(0);
291 #ifdef GMM_USES_MPI
292  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
293 #endif
294 
295  id.job = JOB_INIT;
296  id.par = 1;
297  id.sym = sym ? 2 : 0;
298  id.comm_fortran = USE_COMM_WORLD;
299  mumps_interf<T>::mumps_c(id);
300 
301  if (rank == 0 || distributed) {
302  id.n = int(gmm::mat_nrows(A));
303  if (distributed) {
304  id.nz_loc = int(AA.irn.size());
305  id.irn_loc = &(AA.irn[0]);
306  id.jcn_loc = &(AA.jcn[0]);
307  id.a_loc = (MUMPS_T*)(&(AA.a[0]));
308  } else {
309  id.nz = int(AA.irn.size());
310  id.irn = &(AA.irn[0]);
311  id.jcn = &(AA.jcn[0]);
312  id.a = (MUMPS_T*)(&(AA.a[0]));
313  }
314  }
315 
316  id.ICNTL(1) = -1; // output stream for error messages
317  id.ICNTL(2) = -1; // output stream for other messages
318  id.ICNTL(3) = -1; // output stream for global information
319  id.ICNTL(4) = 0; // verbosity level
320 
321  if (distributed)
322  id.ICNTL(5) = 0; // assembled input matrix (default)
323 
324 // id.ICNTL(14) += 80; // small boost to the workspace size
325 
326  if (distributed)
327  id.ICNTL(18) = 3; // strategy for distributed input matrix
328 
329  id.ICNTL(31) = 1; // only factorization, no solution to follow
330  id.ICNTL(33) = 1; // request determinant calculation
331 
332  id.job = 4; // abalysis (job=1) + factorization (job=2)
333  mumps_interf<T>::mumps_c(id);
334  mumps_error_check(id);
335 
336  T det = real_or_complex(std::complex<R>(id.RINFOG(12),id.RINFOG(13)));
337  exponent = id.INFOG(34);
338 
339  id.job = JOB_END;
340  mumps_interf<T>::mumps_c(id);
341 
342  return det;
343  }
344 
345 #undef ICNTL
346 #undef INFO
347 #undef INFOG
348 #undef RINFOG
349 
350 }
351 
352 
353 #endif // GMM_MUMPS_INTERFACE_H
354 
355 #endif // GMM_USES_MUMPS
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
gmm::nnz
size_type nnz(const L &l)
count the number of non-zero entries of a vector or matrix.
Definition: gmm_blas.h:68
gmm_kernel.h
Include the base gmm files.