GetFEM  5.4.2
getfem_superlu.cc
1 /*===========================================================================
2 
3  Copyright (C) 2004-2020 Julien Pommier
4 
5  This file is a part of GetFEM
6 
7  GetFEM is free software; you can redistribute it and/or modify it
8  under the terms of the GNU Lesser General Public License as published
9  by the Free Software Foundation; either version 3 of the License, or
10  (at your option) any later version along with the GCC Runtime Library
11  Exception either version 3.1 or (at your option) any later version.
12  This program is distributed in the hope that it will be useful, but
13  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15  License and GCC Runtime Library Exception for more details.
16  You should have received a copy of the GNU Lesser General Public License
17  along with this program; if not, write to the Free Software Foundation,
18  Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
19 
20 ===========================================================================*/
21 
22 #include "getfem/getfem_superlu.h"
23 
24 typedef int int_t;
25 
26 /* because SRC/util.h defines TRUE and FALSE ... */
27 #ifdef TRUE
28 # undef TRUE
29 #endif
30 #ifdef FALSE
31 # undef FALSE
32 #endif
33 
34 #include "superlu/slu_Cnames.h"
35 #include "superlu/supermatrix.h"
36 #include "superlu/slu_util.h"
37 
38 namespace SuperLU_S {
39 #include "superlu/slu_sdefs.h"
40 }
41 namespace SuperLU_D {
42 #include "superlu/slu_ddefs.h"
43 }
44 namespace SuperLU_C {
45 #include "superlu/slu_cdefs.h"
46 }
47 namespace SuperLU_Z {
48 #include "superlu/slu_zdefs.h"
49 }
50 
51 
52 
53 namespace gmm {
54 
55  /* interface for Create_CompCol_Matrix */
56 
57  inline void Create_CompCol_Matrix(SuperMatrix *A, int m, int n, int nnz,
58  float *a, int *ir, int *jc) {
59  SuperLU_S::sCreate_CompCol_Matrix(A, m, n, nnz, a, ir, jc,
60  SLU_NC, SLU_S, SLU_GE);
61  }
62 
63  inline void Create_CompCol_Matrix(SuperMatrix *A, int m, int n, int nnz,
64  double *a, int *ir, int *jc) {
65  SuperLU_D::dCreate_CompCol_Matrix(A, m, n, nnz, a, ir, jc,
66  SLU_NC, SLU_D, SLU_GE);
67  }
68 
69  inline void Create_CompCol_Matrix(SuperMatrix *A, int m, int n, int nnz,
70  std::complex<float> *a, int *ir, int *jc) {
71  SuperLU_C::cCreate_CompCol_Matrix(A, m, n, nnz, (SuperLU_C::complex *)(a),
72  ir, jc, SLU_NC, SLU_C, SLU_GE);
73  }
74 
75  inline void Create_CompCol_Matrix(SuperMatrix *A, int m, int n, int nnz,
76  std::complex<double> *a, int *ir, int *jc) {
77  SuperLU_Z::zCreate_CompCol_Matrix(A, m, n, nnz,
78  (SuperLU_Z::doublecomplex *)(a), ir, jc,
79  SLU_NC, SLU_Z, SLU_GE);
80  }
81 
82  /* interface for Create_Dense_Matrix */
83 
84  inline void Create_Dense_Matrix(SuperMatrix *A, int m, int n, float *a, int k)
85  { SuperLU_S::sCreate_Dense_Matrix(A, m, n, a, k, SLU_DN, SLU_S, SLU_GE); }
86  inline void Create_Dense_Matrix(SuperMatrix *A, int m, int n, double *a, int k)
87  { SuperLU_D::dCreate_Dense_Matrix(A, m, n, a, k, SLU_DN, SLU_D, SLU_GE); }
88  inline void Create_Dense_Matrix(SuperMatrix *A, int m, int n,
89  std::complex<float> *a, int k) {
90  SuperLU_C::cCreate_Dense_Matrix(A, m, n, (SuperLU_C::complex *)(a),
91  k, SLU_DN, SLU_C, SLU_GE);
92  }
93  inline void Create_Dense_Matrix(SuperMatrix *A, int m, int n,
94  std::complex<double> *a, int k) {
95  SuperLU_Z::zCreate_Dense_Matrix(A, m, n, (SuperLU_Z::doublecomplex *)(a),
96  k, SLU_DN, SLU_Z, SLU_GE);
97  }
98 
99  /* interface for gssv */
100 
101 #define DECL_GSSV(NAMESPACE,FNAME,FLOATTYPE,KEYTYPE) \
102  inline void SuperLU_gssv(superlu_options_t *options, SuperMatrix *A, int *p, \
103  int *q, SuperMatrix *L, SuperMatrix *U, SuperMatrix *B, \
104  SuperLUStat_t *stats, int *info, KEYTYPE) { \
105  NAMESPACE::FNAME(options, A, p, q, L, U, B, stats, info); \
106  }
107 
108  DECL_GSSV(SuperLU_S,sgssv,float,float)
109  DECL_GSSV(SuperLU_C,cgssv,float,std::complex<float>)
110  DECL_GSSV(SuperLU_D,dgssv,double,double)
111  DECL_GSSV(SuperLU_Z,zgssv,double,std::complex<double>)
112 
113  /* interface for gssvx */
114 
115 #define DECL_GSSVX(NAMESPACE,FNAME,FLOATTYPE,KEYTYPE) \
116  inline float SuperLU_gssvx(superlu_options_t *options, SuperMatrix *A, \
117  int *perm_c, int *perm_r, int *etree, char *equed, \
118  FLOATTYPE *R, FLOATTYPE *C, SuperMatrix *L, \
119  SuperMatrix *U, void *work, int lwork, \
120  SuperMatrix *B, SuperMatrix *X, \
121  FLOATTYPE *recip_pivot_growth, \
122  FLOATTYPE *rcond, FLOATTYPE *ferr, FLOATTYPE *berr, \
123  SuperLUStat_t *stats, int *info, KEYTYPE) { \
124  NAMESPACE::mem_usage_t mem_usage; \
125  NAMESPACE::FNAME(options, A, perm_c, perm_r, etree, equed, R, C, L, \
126  U, work, lwork, B, X, recip_pivot_growth, rcond, \
127  ferr, berr, &mem_usage, stats, info); \
128  return mem_usage.for_lu; /* bytes used by the factor storage */ \
129  }
130 
131  DECL_GSSVX(SuperLU_S,sgssvx,float,float)
132  DECL_GSSVX(SuperLU_C,cgssvx,float,std::complex<float>)
133  DECL_GSSVX(SuperLU_D,dgssvx,double,double)
134  DECL_GSSVX(SuperLU_Z,zgssvx,double,std::complex<double>)
135 
136  /* ********************************************************************* */
137  /* SuperLU solve interface */
138  /* ********************************************************************* */
139 
140  template<typename T>
141  int SuperLU_solve(const gmm::csc_matrix<T> &csc_A, T *sol, T *rhs,
142  double& rcond_, int permc_spec) {
143  /*
144  * Get column permutation vector perm_c[], according to permc_spec:
145  * permc_spec = 0: use the natural ordering
146  * permc_spec = 1: use minimum degree ordering on structure of A'*A
147  * permc_spec = 2: use minimum degree ordering on structure of A'+A
148  * permc_spec = 3: use approximate minimum degree column ordering
149  */
150  typedef typename gmm::number_traits<T>::magnitude_type R;
151 
152  int m = int(mat_nrows(csc_A)), n = int(mat_ncols(csc_A));
153  int nrhs = 1, info = 0, nz = int(nnz(csc_A));
154 
155  GMM_ASSERT1(nz != 0, "Cannot factor a matrix full of zeros!");
156  GMM_ASSERT1(n == m, "Cannot factor a non-square matrix");
157 
158  if ((2 * nz / n) >= m)
159  GMM_WARNING2("CAUTION : it seems that SuperLU has a problem"
160  " for nearly dense sparse matrices");
161 
162  superlu_options_t options;
163  set_default_options(&options);
164  options.ColPerm = NATURAL;
165  options.PrintStat = NO;
166  options.ConditionNumber = YES;
167  switch (permc_spec) {
168  case 1 : options.ColPerm = MMD_ATA; break;
169  case 2 : options.ColPerm = MMD_AT_PLUS_A; break;
170  case 3 : options.ColPerm = COLAMD; break;
171  }
172  SuperLUStat_t stat;
173  StatInit(&stat);
174 
175  SuperMatrix SA, SL, SU, SB, SX; // SuperLU format.
176  Create_CompCol_Matrix(&SA, m, n, nz, const_cast<T*>(&csc_A.pr[0]),
177  const_cast<int *>((const int *)(&csc_A.ir[0])),
178  const_cast<int *>((const int *)(&csc_A.jc[0])));
179  Create_Dense_Matrix(&SB, m, nrhs, &rhs[0], m);
180  Create_Dense_Matrix(&SX, m, nrhs, &sol[0], m);
181 
182  memset(&SL,0,sizeof SL);
183  memset(&SU,0,sizeof SU);
184 
185  std::vector<int> etree(n);
186  char equed[] = "B";
187  std::vector<R> Rscale(m),Cscale(n); // row scale factors
188  std::vector<R> ferr(nrhs), berr(nrhs);
189  R recip_pivot_gross, rcond;
190  std::vector<int> perm_r(m), perm_c(n);
191 
192  SuperLU_gssvx(&options, &SA, &perm_c[0], &perm_r[0],
193  &etree[0] /* output */, equed /* output */,
194  &Rscale[0] /* row scale factors (output) */,
195  &Cscale[0] /* col scale factors (output) */,
196  &SL /* fact L (output)*/, &SU /* fact U (output)*/,
197  NULL /* work */,
198  0 /* lwork: superlu auto allocates (input) */,
199  &SB /* rhs */, &SX /* solution */,
200  &recip_pivot_gross /* reciprocal pivot growth */
201  /* factor max_j( norm(A_j)/norm(U_j) ). */,
202  &rcond /*estimate of the reciprocal condition */
203  /* number of the matrix A after equilibration */,
204  &ferr[0] /* estimated forward error */,
205  &berr[0] /* relative backward error */,
206  &stat, &info, T());
207  rcond_ = rcond;
208  if (SB.Store) Destroy_SuperMatrix_Store(&SB);
209  if (SX.Store) Destroy_SuperMatrix_Store(&SX);
210  if (SA.Store) Destroy_SuperMatrix_Store(&SA);
211  if (SL.Store) Destroy_SuperNode_Matrix(&SL);
212  if (SU.Store) Destroy_CompCol_Matrix(&SU);
213  StatFree(&stat);
214  GMM_ASSERT1(info != -333333333, "SuperLU was cancelled."); // user interruption (for matlab interface)
215 
216  GMM_ASSERT1(info >= 0, "SuperLU solve failed: info =" << info);
217  if (info > 0) GMM_WARNING1("SuperLU solve failed: info =" << info);
218  return info;
219  }
220 
221  template int SuperLU_solve(const gmm::csc_matrix<float> &csc_A, float *sol, float *rhs, double& rcond_, int permc_spec);
222  template int SuperLU_solve(const gmm::csc_matrix<double> &csc_A, double *sol, double *rhs, double& rcond_, int permc_spec);
223  template int SuperLU_solve(const gmm::csc_matrix<std::complex<float> > &csc_A, std::complex<float> *sol, std::complex<float> *rhs, double& rcond_, int permc_spec);
224  template int SuperLU_solve(const gmm::csc_matrix<std::complex<double> > &csc_A, std::complex<double> *sol, std::complex<double> *rhs, double& rcond_, int permc_spec);
225 
226  struct SuperLU_factor_impl_common {
227  mutable SuperMatrix SA, SL, SB, SU, SX;
228  mutable SuperLUStat_t stat;
229  mutable superlu_options_t options;
230  float memory_used;
231  mutable bool is_init;
232  mutable char equed;
233  void free_supermatrix() {
234  if (is_init) {
235  if (SB.Store) Destroy_SuperMatrix_Store(&SB);
236  if (SX.Store) Destroy_SuperMatrix_Store(&SX);
237  if (SA.Store) Destroy_SuperMatrix_Store(&SA);
238  if (SL.Store) Destroy_SuperNode_Matrix(&SL);
239  if (SU.Store) Destroy_CompCol_Matrix(&SU);
240  }
241  }
242  SuperLU_factor_impl_common() : is_init(false) {}
243  virtual ~SuperLU_factor_impl_common() { free_supermatrix(); }
244  };
245 
246  template <typename T> struct SuperLU_factor_impl : public SuperLU_factor_impl_common {
247  typedef typename gmm::number_traits<T>::magnitude_type R;
248 
249  std::vector<int> etree, perm_r, perm_c;
250  std::vector<R> Rscale, Cscale;
251  std::vector<R> ferr, berr;
252  std::vector<T> rhs;
253  std::vector<T> sol;
254  void build_with(const gmm::csc_matrix<T> &A, int permc_spec);
255  void solve(int transp);
256  };
257 
258  template <typename T>
259  void SuperLU_factor_impl<T>::build_with(const gmm::csc_matrix<T> &A, int permc_spec) {
260  /*
261  * Get column permutation vector perm_c[], according to permc_spec:
262  * permc_spec = 0: use the natural ordering
263  * permc_spec = 1: use minimum degree ordering on structure of A'*A
264  * permc_spec = 2: use minimum degree ordering on structure of A'+A
265  * permc_spec = 3: use approximate minimum degree column ordering
266  */
267  free_supermatrix();
268  int n = int(mat_nrows(A)), m = int(mat_ncols(A)), info = 0;
269 
270  rhs.resize(m); sol.resize(m);
271  gmm::clear(rhs);
272  int nz = int(nnz(A));
273 
274  GMM_ASSERT1(nz != 0, "Cannot factor a matrix full of zeros!");
275  GMM_ASSERT1(n == m, "Cannot factor a non-square matrix");
276 
277  set_default_options(&options);
278  options.ColPerm = NATURAL;
279  options.PrintStat = NO;
280  options.ConditionNumber = NO;
281  switch (permc_spec) {
282  case 1 : options.ColPerm = MMD_ATA; break;
283  case 2 : options.ColPerm = MMD_AT_PLUS_A; break;
284  case 3 : options.ColPerm = COLAMD; break;
285  }
286  StatInit(&stat);
287 
288  Create_CompCol_Matrix(&SA, m, n, nz, const_cast<T*>(&A.pr[0]),
289  const_cast<int *>((const int *)(&A.ir[0])),
290  const_cast<int *>((const int *)(&A.jc[0])));
291  Create_Dense_Matrix(&SB, m, 0, &rhs[0], m);
292  Create_Dense_Matrix(&SX, m, 0, &sol[0], m);
293  memset(&SL,0,sizeof SL);
294  memset(&SU,0,sizeof SU);
295 
296 
297  equed = 'B';
298  Rscale.resize(m); Cscale.resize(n); etree.resize(n);
299  ferr.resize(1); berr.resize(1);
300  R recip_pivot_gross, rcond;
301  perm_r.resize(m); perm_c.resize(n);
302  memory_used = SuperLU_gssvx(&options, &SA, &perm_c[0], &perm_r[0],
303  &etree[0] /* output */, &equed /* output */,
304  &Rscale[0] /* row scale factors (output) */,
305  &Cscale[0] /* col scale factors (output) */,
306  &SL /* fact L (output)*/, &SU /* fact U (output)*/,
307  NULL /* work */,
308  0 /* lwork: superlu auto allocates (input) */,
309  &SB /* rhs */, &SX /* solution */,
310  &recip_pivot_gross /* reciprocal pivot growth */
311  /* factor max_j( norm(A_j)/norm(U_j) ). */,
312  &rcond /*estimate of the reciprocal condition */
313  /* number of the matrix A after equilibration */,
314  &ferr[0] /* estimated forward error */,
315  &berr[0] /* relative backward error */,
316  &stat, &info, T());
317 
318  Destroy_SuperMatrix_Store(&SB);
319  Destroy_SuperMatrix_Store(&SX);
320  Create_Dense_Matrix(&SB, m, 1, &rhs[0], m);
321  Create_Dense_Matrix(&SX, m, 1, &sol[0], m);
322  StatFree(&stat);
323 
324  GMM_ASSERT1(info != -333333333, "SuperLU was cancelled.");
325  GMM_ASSERT1(info == 0, "SuperLU solve failed: info=" << info);
326  is_init = true;
327  }
328 
329  template <typename T>
330  void SuperLU_factor_impl<T>::solve(int transp) {
331  options.Fact = FACTORED;
332  options.IterRefine = NOREFINE;
333  switch (transp) {
334  case SuperLU_factor<T>::LU_NOTRANSP: options.Trans = NOTRANS; break;
335  case SuperLU_factor<T>::LU_TRANSP: options.Trans = TRANS; break;
336  case SuperLU_factor<T>::LU_CONJUGATED: options.Trans = CONJ; break;
337  default: GMM_ASSERT1(false, "invalid value for transposition option");
338  }
339  StatInit(&stat);
340  int info = 0;
341  R recip_pivot_gross, rcond;
342  SuperLU_gssvx(&options, &SA, &perm_c[0], &perm_r[0],
343  &etree[0] /* output */, &equed /* output */,
344  &Rscale[0] /* row scale factors (output) */,
345  &Cscale[0] /* col scale factors (output) */,
346  &SL /* fact L (output)*/, &SU /* fact U (output)*/,
347  NULL /* work */,
348  0 /* lwork: superlu auto allocates (input) */,
349  &SB /* rhs */, &SX /* solution */,
350  &recip_pivot_gross /* reciprocal pivot growth */
351  /* factor max_j( norm(A_j)/norm(U_j) ). */,
352  &rcond /*estimate of the reciprocal condition */
353  /* number of the matrix A after equilibration */,
354  &ferr[0] /* estimated forward error */,
355  &berr[0] /* relative backward error */,
356  &stat, &info, T());
357  StatFree(&stat);
358  GMM_ASSERT1(info == 0, "SuperLU solve failed: info=" << info);
359  }
360 
361  template<typename T> void
362  SuperLU_factor<T>::build_with(const gmm::csc_matrix<T> &A, int permc_spec) {
363  ((SuperLU_factor_impl<T>*)impl.get())->build_with(A,permc_spec);
364  }
365 
366  template<typename T> void
367  SuperLU_factor<T>::solve(int transp) const {
368  ((SuperLU_factor_impl<T>*)impl.get())->solve(transp);
369  }
370 
371  template<typename T> std::vector<T> &
372  SuperLU_factor<T>::sol() const {
373  return ((SuperLU_factor_impl<T>*)impl.get())->sol;
374  }
375 
376  template<typename T> std::vector<T> &
377  SuperLU_factor<T>::rhs() const {
378  return ((SuperLU_factor_impl<T>*)impl.get())->rhs;
379  }
380 
381  template<typename T>
382  SuperLU_factor<T>::SuperLU_factor() {
383  impl = std::make_shared<SuperLU_factor_impl<T>>();
384  }
385 
386  template<typename T>
387  SuperLU_factor<T>::SuperLU_factor(const SuperLU_factor& other) {
388  impl = std::make_shared<SuperLU_factor_impl<T>>();
389  GMM_ASSERT1(!(other.impl->is_init),
390  "copy of initialized SuperLU_factor is forbidden");
391  other.impl->is_init = false;
392  }
393 
394  template<typename T> SuperLU_factor<T>&
395  SuperLU_factor<T>::operator=(const SuperLU_factor& other) {
396  GMM_ASSERT1(!(other.impl->is_init) && !(impl->is_init),
397  "assignment of initialized SuperLU_factor is forbidden");
398  return *this;
399  }
400 
401  template<typename T> float
402  SuperLU_factor<T>::memsize() const {
403  return impl->memory_used;
404  }
405 
406  /* void force_instantiation() {
407  SuperLU_factor<float> a;
408  SuperLU_factor<double> b;
409  SuperLU_factor<std::complex<float> > c;
410  SuperLU_factor<std::complex<double> > d;
411  //a = 0; b = 0; c = 0; d = 0;
412  }
413  */
414 }
415 
416 template class gmm::SuperLU_factor<float>;
417 template class gmm::SuperLU_factor<double>;
418 template class gmm::SuperLU_factor<std::complex<float> >;
419 template class gmm::SuperLU_factor<std::complex<double> >;
420 
421 static int (*superlu_callback)();
422 
423 /* this one is called from superlu (see dcolumn_bmod) */
424 extern "C" int handle_getfem_callback() {
425  if (superlu_callback) return superlu_callback();
426  else return 0;
427 }
428 
429 extern "C" void set_superlu_callback(int (*cb)()) {
430  superlu_callback = cb;
431 }
gmm::clear
void clear(L &l)
clear (fill with zeros) a vector or matrix.
Definition: gmm_blas.h:59
getfem_superlu.h
SuperLU interface for getfem.
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::SuperLU_factor::build_with
void build_with(const MAT &A, int permc_spec=3)
Do the factorization of the supplied sparse matrix.
Definition: getfem_superlu.h:90
gmm::SuperLU_factor
Factorization of a sparse matrix with SuperLU.
Definition: getfem_superlu.h:84
gmm::SuperLU_factor::solve
void solve(const VECTX &X, const VECTB &B, int transp=LU_NOTRANSP) const
After factorization, do the triangular solves.
Definition: getfem_superlu.h:103