GetFEM  5.4.3
gmm_superlu_interface.h
Go to the documentation of this file.
1 /* -*- c++ -*- (enables emacs c++ mode) */
2 /*===========================================================================
3 
4  Copyright (C) 2003-2024 Yves Renard
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_superlu_interface.h
33  @author Yves Renard <[email protected]>
34  @date October 17, 2003.
35  @brief Interface with SuperLU (LU direct solver for sparse matrices).
36 */
37 #if defined(GMM_USES_SUPERLU)
38 
39 #ifndef GMM_SUPERLU_INTERFACE_H
40 #define GMM_SUPERLU_INTERFACE_H
41 
42 #include "gmm_kernel.h"
43 
44 
45 typedef int int_t;
46 
47 /* because slu_util.h defines TRUE, FALSE, EMPTY ... */
48 #ifdef TRUE
49 # undef TRUE
50 #endif
51 #ifdef FALSE
52 # undef FALSE
53 #endif
54 #ifdef EMPTY
55 # undef EMPTY
56 #endif
57 
58 #if defined(GMM_NO_SUPERLU_INCLUDE_SUBDIR)
59 
60  #include "slu_Cnames.h"
61  #include "supermatrix.h"
62  #include "slu_util.h"
63 
64  namespace SuperLU_S {
65  #include "slu_sdefs.h"
66  }
67  namespace SuperLU_D {
68  #include "slu_ddefs.h"
69  }
70  namespace SuperLU_C {
71  #include "slu_cdefs.h"
72  }
73  namespace SuperLU_Z {
74  #include "slu_zdefs.h"
75  }
76 
77 #else
78 
79  #include "superlu/slu_Cnames.h"
80  #include "superlu/supermatrix.h"
81  #include "superlu/slu_util.h"
82 
83  namespace SuperLU_S {
84  #include "superlu/slu_sdefs.h"
85  }
86  namespace SuperLU_D {
87  #include "superlu/slu_ddefs.h"
88  }
89  namespace SuperLU_C {
90  #include "superlu/slu_cdefs.h"
91  }
92  namespace SuperLU_Z {
93  #include "superlu/slu_zdefs.h"
94  }
95 
96 #endif
97 
98 #if (SUPERLU_MAJOR_VERSION > 6)
99 # define SuperLuComplexFloat SuperLU_C::singlecomplex
100 #else
101 # define SuperLuComplexFloat SuperLU_C::complex
102 #endif
103 
104 namespace gmm {
105 
106  /* interface for Create_CompCol_Matrix */
107 
108  inline void Create_CompCol_Matrix(SuperMatrix *A, int m, int n, int nnz,
109  float *a, int *ir, int *jc) {
110  SuperLU_S::sCreate_CompCol_Matrix(A, m, n, nnz, a, ir, jc,
111  SLU_NC, SLU_S, SLU_GE);
112  }
113 
114  inline void Create_CompCol_Matrix(SuperMatrix *A, int m, int n, int nnz,
115  double *a, int *ir, int *jc) {
116  SuperLU_D::dCreate_CompCol_Matrix(A, m, n, nnz, a, ir, jc,
117  SLU_NC, SLU_D, SLU_GE);
118  }
119 
120  inline void Create_CompCol_Matrix(SuperMatrix *A, int m, int n, int nnz,
121  std::complex<float> *a, int *ir, int *jc) {
122  SuperLU_C::cCreate_CompCol_Matrix(A, m, n, nnz, (SuperLuComplexFloat *)(a),
123  ir, jc, SLU_NC, SLU_C, SLU_GE);
124  }
125 
126  inline void Create_CompCol_Matrix(SuperMatrix *A, int m, int n, int nnz,
127  std::complex<double> *a, int *ir, int *jc) {
128  SuperLU_Z::zCreate_CompCol_Matrix(A, m, n, nnz,
129  (SuperLU_Z::doublecomplex *)(a), ir, jc,
130  SLU_NC, SLU_Z, SLU_GE);
131  }
132 
133  /* interface for Create_Dense_Matrix */
134 
135  inline void Create_Dense_Matrix(SuperMatrix *A, int m, int n, float *a, int k)
136  { SuperLU_S::sCreate_Dense_Matrix(A, m, n, a, k, SLU_DN, SLU_S, SLU_GE); }
137  inline void Create_Dense_Matrix(SuperMatrix *A, int m, int n, double *a,
138  int k)
139  { SuperLU_D::dCreate_Dense_Matrix(A, m, n, a, k, SLU_DN, SLU_D, SLU_GE); }
140  inline void Create_Dense_Matrix(SuperMatrix *A, int m, int n,
141  std::complex<float> *a, int k) {
142  SuperLU_C::cCreate_Dense_Matrix(A, m, n, (SuperLuComplexFloat *)(a),
143  k, SLU_DN, SLU_C, SLU_GE);
144  }
145  inline void Create_Dense_Matrix(SuperMatrix *A, int m, int n,
146  std::complex<double> *a, int k) {
147  SuperLU_Z::zCreate_Dense_Matrix(A, m, n, (SuperLU_Z::doublecomplex *)(a),
148  k, SLU_DN, SLU_Z, SLU_GE);
149  }
150 
151  /* interface for gssv */
152 
153 #define DECL_GSSV(NAMESPACE,FNAME,KEYTYPE) \
154  inline void SuperLU_gssv(superlu_options_t *options, SuperMatrix *A, int *p, \
155  int *q, SuperMatrix *L, SuperMatrix *U, SuperMatrix *B, \
156  SuperLUStat_t *stats, int *info, KEYTYPE) { \
157  NAMESPACE::FNAME(options, A, p, q, L, U, B, stats, info); \
158  }
159 
160  DECL_GSSV(SuperLU_S,sgssv, float)
161  DECL_GSSV(SuperLU_C,cgssv, std::complex<float>)
162  DECL_GSSV(SuperLU_D,dgssv, double)
163  DECL_GSSV(SuperLU_Z,zgssv, std::complex<double>)
164 
165  /* interface for gssvx */
166 
167 #define DECL_GSSVX(NAMESPACE,FNAME,FLOATTYPE,KEYTYPE) \
168  inline float SuperLU_gssvx(superlu_options_t *options, SuperMatrix *A, \
169  int *perm_c, int *perm_r, int *etree, \
170  char *equed, \
171  FLOATTYPE *R, FLOATTYPE *C, SuperMatrix *L, \
172  SuperMatrix *U, void *work, int lwork, \
173  SuperMatrix *B, SuperMatrix *X, \
174  FLOATTYPE *recip_pivot_growth, \
175  FLOATTYPE *rcond, FLOATTYPE *ferr, \
176  FLOATTYPE *berr, \
177  SuperLUStat_t *stats, int *info, KEYTYPE) { \
178  mem_usage_t mem_usage; \
179  GlobalLU_t Glu; \
180  NAMESPACE::FNAME(options, A, perm_c, perm_r, etree, equed, R, C, L, \
181  U, work, lwork, B, X, recip_pivot_growth, rcond, \
182  ferr, berr, &Glu, &mem_usage, stats, info); \
183  return mem_usage.for_lu; /* bytes used by the factor storage */ \
184  }
185 
186  DECL_GSSVX(SuperLU_S, sgssvx, float, float)
187  DECL_GSSVX(SuperLU_C, cgssvx, float, std::complex<float>)
188  DECL_GSSVX(SuperLU_D, dgssvx, double, double)
189  DECL_GSSVX(SuperLU_Z, zgssvx, double, std::complex<double>)
190 
191  /* ********************************************************************* */
192  /* SuperLU solve interface */
193  /* ********************************************************************* */
194 
195  template <typename MAT, typename VECTX, typename VECTB>
196  int SuperLU_solve(const MAT &A, const VECTX &X, const VECTB &B,
197  double& rcond_, int permc_spec = 3) {
198  /*
199  * Get column permutation vector perm_c[], according to permc_spec:
200  * permc_spec = 0: use the natural ordering
201  * permc_spec = 1: use minimum degree ordering on structure of A'*A
202  * permc_spec = 2: use minimum degree ordering on structure of A'+A
203  * permc_spec = 3: use approximate minimum degree column ordering
204  */
205  typedef typename linalg_traits<MAT>::value_type T;
206  typedef typename number_traits<T>::magnitude_type R;
207 
208  int m = int(mat_nrows(A)), n = int(mat_ncols(A)), nrhs = 1, info = 0;
209 
210  csc_matrix<T> csc_A(m, n);
211  gmm::copy(A, csc_A);
212  std::vector<T> rhs(m), sol(m);
213  gmm::copy(B, rhs);
214 
215  int nz = int(nnz(csc_A));
216  if ((2 * nz / n) >= m)
217  GMM_WARNING2("CAUTION : it seems that SuperLU has a problem"
218  " for nearly dense sparse matrices");
219 
220  superlu_options_t options;
221  set_default_options(&options);
222  options.ColPerm = NATURAL;
223  options.PrintStat = NO;
224  options.ConditionNumber = YES;
225  switch (permc_spec) {
226  case 1 : options.ColPerm = MMD_ATA; break;
227  case 2 : options.ColPerm = MMD_AT_PLUS_A; break;
228  case 3 : options.ColPerm = COLAMD; break;
229  }
230  SuperLUStat_t stat;
231  StatInit(&stat);
232 
233  SuperMatrix SA, SL, SU, SB, SX; // SuperLU format.
234  Create_CompCol_Matrix(&SA, m, n, nz, (T *)(&csc_A.pr[0]),
235  (int *)(&csc_A.ir[0]),
236  (int *)(&csc_A.jc[0]));
237  Create_Dense_Matrix(&SB, m, nrhs, &rhs[0], m);
238  Create_Dense_Matrix(&SX, m, nrhs, &sol[0], m);
239  memset(&SL,0,sizeof SL);
240  memset(&SU,0,sizeof SU);
241 
242  std::vector<int> etree(n);
243  char equed[] = "B";
244  std::vector<R> Rscale(m),Cscale(n); // row scale factors
245  std::vector<R> ferr(nrhs), berr(nrhs);
246  R recip_pivot_gross, rcond;
247  std::vector<int> perm_r(m), perm_c(n);
248 
249  SuperLU_gssvx(&options, &SA, &perm_c[0], &perm_r[0],
250  &etree[0] /* output */, equed /* output */,
251  &Rscale[0] /* row scale factors (output) */,
252  &Cscale[0] /* col scale factors (output) */,
253  &SL /* fact L (output)*/, &SU /* fact U (output)*/,
254  NULL /* work */,
255  0 /* lwork: superlu auto allocates (input) */,
256  &SB /* rhs */, &SX /* solution */,
257  &recip_pivot_gross /* reciprocal pivot growth */
258  /* factor max_j( norm(A_j)/norm(U_j) ). */,
259  &rcond /*estimate of the reciprocal condition */
260  /* number of the matrix A after equilibration */,
261  &ferr[0] /* estimated forward error */,
262  &berr[0] /* relative backward error */,
263  &stat, &info, T());
264  rcond_ = rcond;
265  if (SB.Store) Destroy_SuperMatrix_Store(&SB);
266  if (SX.Store) Destroy_SuperMatrix_Store(&SX);
267  if (SA.Store) Destroy_SuperMatrix_Store(&SA);
268  if (SL.Store) Destroy_SuperNode_Matrix(&SL);
269  if (SU.Store) Destroy_CompCol_Matrix(&SU);
270  StatFree(&stat);
271  GMM_ASSERT1(info != -333333333, "SuperLU was cancelled."); // user interruption (for matlab interface)
272 
273  GMM_ASSERT1(info >= 0, "SuperLU solve failed: info =" << info);
274  if (info > 0) GMM_WARNING1("SuperLU solve failed: info =" << info);
275  gmm::copy(sol, const_cast<VECTX &>(X));
276  return info;
277  }
278 
279  template <class T>
280  class SuperLU_factor {
281  typedef typename number_traits<T>::magnitude_type R;
282 
283  csc_matrix<T> csc_A;
284  mutable SuperMatrix SA, SL, SB, SU, SX;
285  mutable SuperLUStat_t stat;
286  mutable superlu_options_t options;
287  float memory_used;
288  mutable std::vector<int> etree, perm_r, perm_c;
289  mutable std::vector<R> Rscale, Cscale;
290  mutable std::vector<R> ferr, berr;
291  mutable std::vector<T> rhs;
292  mutable std::vector<T> sol;
293  mutable bool is_init;
294  mutable char equed;
295 
296  public :
297  enum { LU_NOTRANSP, LU_TRANSP, LU_CONJUGATED };
298  void free_supermatrix() {
299  if (is_init) {
300  if (SB.Store) Destroy_SuperMatrix_Store(&SB);
301  if (SX.Store) Destroy_SuperMatrix_Store(&SX);
302  if (SA.Store) Destroy_SuperMatrix_Store(&SA);
303  if (SL.Store) Destroy_SuperNode_Matrix(&SL);
304  if (SU.Store) Destroy_CompCol_Matrix(&SU);
305  }
306  }
307  template <class MAT> void build_with(const MAT &A, int permc_spec = 3);
308  template <typename VECTX, typename VECTB>
309  /* transp = LU_NOTRANSP -> solves Ax = B
310  transp = LU_TRANSP -> solves A'x = B
311  transp = LU_CONJUGATED -> solves conj(A)X = B */
312  void solve(const VECTX &X_, const VECTB &B, int transp=LU_NOTRANSP) const;
313  SuperLU_factor() { is_init = false; }
314  SuperLU_factor(const SuperLU_factor& other) {
315  GMM_ASSERT2(!(other.is_init),
316  "copy of initialized SuperLU_factor is forbidden");
317  is_init = false;
318  }
319  SuperLU_factor& operator=(const SuperLU_factor& other) {
320  GMM_ASSERT2(!(other.is_init) && !is_init,
321  "assignment of initialized SuperLU_factor is forbidden");
322  return *this;
323  }
324  ~SuperLU_factor() { free_supermatrix(); }
325  float memsize() { return memory_used; }
326  };
327 
328 
329  template <class T> template <class MAT>
330  void SuperLU_factor<T>::build_with(const MAT &A, int permc_spec) {
331  /*
332  * Get column permutation vector perm_c[], according to permc_spec:
333  * permc_spec = 0: use the natural ordering
334  * permc_spec = 1: use minimum degree ordering on structure of A'*A
335  * permc_spec = 2: use minimum degree ordering on structure of A'+A
336  * permc_spec = 3: use approximate minimum degree column ordering
337  */
338  free_supermatrix();
339  int n = int(mat_nrows(A)), m = int(mat_ncols(A)), info = 0;
340  csc_A.init_with(A);
341 
342  rhs.resize(m); sol.resize(m);
343  gmm::clear(rhs);
344  int nz = int(nnz(csc_A));
345 
346  set_default_options(&options);
347  options.ColPerm = NATURAL;
348  options.PrintStat = NO;
349  options.ConditionNumber = NO;
350  switch (permc_spec) {
351  case 1 : options.ColPerm = MMD_ATA; break;
352  case 2 : options.ColPerm = MMD_AT_PLUS_A; break;
353  case 3 : options.ColPerm = COLAMD; break;
354  }
355  StatInit(&stat);
356 
357  Create_CompCol_Matrix(&SA, m, n, nz, (T *)(&csc_A.pr[0]),
358  (int *)(&csc_A.ir[0]),
359  (int *)(&csc_A.jc[0]));
360  Create_Dense_Matrix(&SB, m, 0, &rhs[0], m);
361  Create_Dense_Matrix(&SX, m, 0, &sol[0], m);
362  memset(&SL,0,sizeof SL);
363  memset(&SU,0,sizeof SU);
364  equed = 'B';
365  Rscale.resize(m); Cscale.resize(n); etree.resize(n);
366  ferr.resize(1); berr.resize(1);
367  R recip_pivot_gross, rcond;
368  perm_r.resize(m); perm_c.resize(n);
369  memory_used = SuperLU_gssvx(&options, &SA, &perm_c[0], &perm_r[0],
370  &etree[0] /* output */, &equed /* output */,
371  &Rscale[0] /* row scale factors (output) */,
372  &Cscale[0] /* col scale factors (output) */,
373  &SL /* fact L (output) */,
374  &SU /* fact U (output) */,
375  NULL /* work */,
376  0 /* lwork: superlu auto allocates (input) */,
377  &SB /* rhs */, &SX /* solution */,
378  &recip_pivot_gross /* reciprocal pivot growth*/
379  /* factor max_j(norm(A_j)/norm(U_j)). */,
380  &rcond /*estimate of the reciprocal condition*/
381  /* number of the matrix A after equilibration*/,
382  &ferr[0] /* estimated forward error */,
383  &berr[0] /* relative backward error */,
384  &stat, &info, T());
385 
386  Destroy_SuperMatrix_Store(&SB);
387  Destroy_SuperMatrix_Store(&SX);
388  Create_Dense_Matrix(&SB, m, 1, &rhs[0], m);
389  Create_Dense_Matrix(&SX, m, 1, &sol[0], m);
390  StatFree(&stat);
391 
392  GMM_ASSERT1(info != -333333333, "SuperLU was cancelled.");
393  GMM_ASSERT1(info == 0, "SuperLU solve failed: info=" << info);
394  is_init = true;
395  }
396 
397  template <class T> template <typename VECTX, typename VECTB>
398  void SuperLU_factor<T>::solve(const VECTX &X, const VECTB &B,
399  int transp) const {
400  gmm::copy(B, rhs);
401  options.Fact = FACTORED;
402  options.IterRefine = NOREFINE;
403  switch (transp) {
404  case LU_NOTRANSP: options.Trans = NOTRANS; break;
405  case LU_TRANSP: options.Trans = TRANS; break;
406  case LU_CONJUGATED: options.Trans = CONJ; break;
407  default: GMM_ASSERT1(false, "invalid value for transposition option");
408  }
409  StatInit(&stat);
410  int info = 0;
411  R recip_pivot_gross, rcond;
412  SuperLU_gssvx(&options, &SA, &perm_c[0], &perm_r[0],
413  &etree[0] /* output */, &equed /* output */,
414  &Rscale[0] /* row scale factors (output) */,
415  &Cscale[0] /* col scale factors (output) */,
416  &SL /* fact L (output)*/, &SU /* fact U (output)*/,
417  NULL /* work */,
418  0 /* lwork: superlu auto allocates (input) */,
419  &SB /* rhs */, &SX /* solution */,
420  &recip_pivot_gross /* reciprocal pivot growth */
421  /* factor max_j( norm(A_j)/norm(U_j) ). */,
422  &rcond /*estimate of the reciprocal condition */
423  /* number of the matrix A after equilibration */,
424  &ferr[0] /* estimated forward error */,
425  &berr[0] /* relative backward error */,
426  &stat, &info, T());
427  StatFree(&stat);
428  GMM_ASSERT1(info == 0, "SuperLU solve failed: info=" << info);
429  gmm::copy(sol, const_cast<VECTX &>(X));
430  }
431 
432  template <typename T, typename V1, typename V2> inline
433  void mult(const SuperLU_factor<T>& P, const V1 &v1, const V2 &v2) {
434  P.solve(v2,v1);
435  }
436 
437  template <typename T, typename V1, typename V2> inline
438  void transposed_mult(const SuperLU_factor<T>& P,const V1 &v1,const V2 &v2) {
439  P.solve(v2, v1, SuperLU_factor<T>::LU_TRANSP);
440  }
441 
442 }
443 
444 #ifdef TRUE
445 # undef TRUE
446 #endif
447 #ifdef FALSE
448 # undef FALSE
449 #endif
450 #ifdef EMPTY
451 # undef EMPTY
452 #endif
453 
454 #endif // GMM_SUPERLU_INTERFACE_H
455 
456 #endif // GMM_USES_SUPERLU
gmm::nnz
size_type nnz(const L &l)
count the number of non-zero entries of a vector or matrix.
Definition: gmm_blas.h:69
gmm::copy
void copy(const L1 &l1, L2 &l2)
*‍/
Definition: gmm_blas.h:978
gmm::clear
void clear(L &l)
clear (fill with zeros) a vector or matrix.
Definition: gmm_blas.h:59
gmm::mult
void mult(const L1 &l1, const L2 &l2, L3 &l3)
*‍/
Definition: gmm_blas.h:1664
gmm_kernel.h
Include the base gmm files.