Welcome to the NetCologne GmbH open source mirroring service!

This machine mirrors various open-source projects. 20 Gbit/s uplink.

If there are any issues or you want another project mirrored, please contact mirror-service -=AT=- netcologne DOT de !

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