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