SuperLUSupport.h
Go to the documentation of this file.
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2008-2015 Gael Guennebaud <gael.guennebaud@inria.fr>
5 //
6 // This Source Code Form is subject to the terms of the Mozilla
7 // Public License v. 2.0. If a copy of the MPL was not distributed
8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9 
10 #ifndef EIGEN_SUPERLUSUPPORT_H
11 #define EIGEN_SUPERLUSUPPORT_H
12 
13 namespace Eigen {
14 
15 #if defined(SUPERLU_MAJOR_VERSION) && (SUPERLU_MAJOR_VERSION >= 5)
16 #define DECL_GSSVX(PREFIX,FLOATTYPE,KEYTYPE) \
17  extern "C" { \
18  extern void PREFIX##gssvx(superlu_options_t *, SuperMatrix *, int *, int *, int *, \
19  char *, FLOATTYPE *, FLOATTYPE *, SuperMatrix *, SuperMatrix *, \
20  void *, int, SuperMatrix *, SuperMatrix *, \
21  FLOATTYPE *, FLOATTYPE *, FLOATTYPE *, FLOATTYPE *, \
22  GlobalLU_t *, mem_usage_t *, SuperLUStat_t *, int *); \
23  } \
24  inline float SuperLU_gssvx(superlu_options_t *options, SuperMatrix *A, \
25  int *perm_c, int *perm_r, int *etree, char *equed, \
26  FLOATTYPE *R, FLOATTYPE *C, SuperMatrix *L, \
27  SuperMatrix *U, void *work, int lwork, \
28  SuperMatrix *B, SuperMatrix *X, \
29  FLOATTYPE *recip_pivot_growth, \
30  FLOATTYPE *rcond, FLOATTYPE *ferr, FLOATTYPE *berr, \
31  SuperLUStat_t *stats, int *info, KEYTYPE) { \
32  mem_usage_t mem_usage; \
33  GlobalLU_t gLU; \
34  PREFIX##gssvx(options, A, perm_c, perm_r, etree, equed, R, C, L, \
35  U, work, lwork, B, X, recip_pivot_growth, rcond, \
36  ferr, berr, &gLU, &mem_usage, stats, info); \
37  return mem_usage.for_lu; /* bytes used by the factor storage */ \
38  }
39 #else // version < 5.0
40 #define DECL_GSSVX(PREFIX,FLOATTYPE,KEYTYPE) \
41  extern "C" { \
42  extern void PREFIX##gssvx(superlu_options_t *, SuperMatrix *, int *, int *, int *, \
43  char *, FLOATTYPE *, FLOATTYPE *, SuperMatrix *, SuperMatrix *, \
44  void *, int, SuperMatrix *, SuperMatrix *, \
45  FLOATTYPE *, FLOATTYPE *, FLOATTYPE *, FLOATTYPE *, \
46  mem_usage_t *, SuperLUStat_t *, int *); \
47  } \
48  inline float SuperLU_gssvx(superlu_options_t *options, SuperMatrix *A, \
49  int *perm_c, int *perm_r, int *etree, char *equed, \
50  FLOATTYPE *R, FLOATTYPE *C, SuperMatrix *L, \
51  SuperMatrix *U, void *work, int lwork, \
52  SuperMatrix *B, SuperMatrix *X, \
53  FLOATTYPE *recip_pivot_growth, \
54  FLOATTYPE *rcond, FLOATTYPE *ferr, FLOATTYPE *berr, \
55  SuperLUStat_t *stats, int *info, KEYTYPE) { \
56  mem_usage_t mem_usage; \
57  PREFIX##gssvx(options, A, perm_c, perm_r, etree, equed, R, C, L, \
58  U, work, lwork, B, X, recip_pivot_growth, rcond, \
59  ferr, berr, &mem_usage, stats, info); \
60  return mem_usage.for_lu; /* bytes used by the factor storage */ \
61  }
62 #endif
63 
64 DECL_GSSVX(s,float,float)
65 DECL_GSSVX(c,float,std::complex<float>)
66 DECL_GSSVX(d,double,double)
67 DECL_GSSVX(z,double,std::complex<double>)
68 
69 #ifdef MILU_ALPHA
70 #define EIGEN_SUPERLU_HAS_ILU
71 #endif
72 
73 #ifdef EIGEN_SUPERLU_HAS_ILU
74 
75 // similarly for the incomplete factorization using gsisx
76 #define DECL_GSISX(PREFIX,FLOATTYPE,KEYTYPE) \
77  extern "C" { \
78  extern void PREFIX##gsisx(superlu_options_t *, SuperMatrix *, int *, int *, int *, \
79  char *, FLOATTYPE *, FLOATTYPE *, SuperMatrix *, SuperMatrix *, \
80  void *, int, SuperMatrix *, SuperMatrix *, FLOATTYPE *, FLOATTYPE *, \
81  mem_usage_t *, SuperLUStat_t *, int *); \
82  } \
83  inline float SuperLU_gsisx(superlu_options_t *options, SuperMatrix *A, \
84  int *perm_c, int *perm_r, int *etree, char *equed, \
85  FLOATTYPE *R, FLOATTYPE *C, SuperMatrix *L, \
86  SuperMatrix *U, void *work, int lwork, \
87  SuperMatrix *B, SuperMatrix *X, \
88  FLOATTYPE *recip_pivot_growth, \
89  FLOATTYPE *rcond, \
90  SuperLUStat_t *stats, int *info, KEYTYPE) { \
91  mem_usage_t mem_usage; \
92  PREFIX##gsisx(options, A, perm_c, perm_r, etree, equed, R, C, L, \
93  U, work, lwork, B, X, recip_pivot_growth, rcond, \
94  &mem_usage, stats, info); \
95  return mem_usage.for_lu; /* bytes used by the factor storage */ \
96  }
97 
98 DECL_GSISX(s,float,float)
99 DECL_GSISX(c,float,std::complex<float>)
100 DECL_GSISX(d,double,double)
101 DECL_GSISX(z,double,std::complex<double>)
102 
103 #endif
104 
105 template<typename MatrixType>
107 
115 struct SluMatrix : SuperMatrix
116 {
118  {
119  Store = &storage;
120  }
121 
123  : SuperMatrix(other)
124  {
125  Store = &storage;
126  storage = other.storage;
127  }
128 
130  {
131  SuperMatrix::operator=(static_cast<const SuperMatrix&>(other));
132  Store = &storage;
133  storage = other.storage;
134  return *this;
135  }
136 
137  struct
138  {
139  union {int nnz;int lda;};
140  void *values;
141  int *innerInd;
142  int *outerInd;
143  } storage;
144 
145  void setStorageType(Stype_t t)
146  {
147  Stype = t;
148  if (t==SLU_NC || t==SLU_NR || t==SLU_DN)
149  Store = &storage;
150  else
151  {
152  eigen_assert(false && "storage type not supported");
153  Store = 0;
154  }
155  }
156 
157  template<typename Scalar>
159  {
161  Dtype = SLU_S;
163  Dtype = SLU_D;
164  else if (internal::is_same<Scalar,std::complex<float> >::value)
165  Dtype = SLU_C;
166  else if (internal::is_same<Scalar,std::complex<double> >::value)
167  Dtype = SLU_Z;
168  else
169  {
170  eigen_assert(false && "Scalar type not supported by SuperLU");
171  }
172  }
173 
174  template<typename MatrixType>
176  {
177  MatrixType& mat(_mat.derived());
178  eigen_assert( ((MatrixType::Flags&RowMajorBit)!=RowMajorBit) && "row-major dense matrices are not supported by SuperLU");
179  SluMatrix res;
180  res.setStorageType(SLU_DN);
181  res.setScalarType<typename MatrixType::Scalar>();
182  res.Mtype = SLU_GE;
183 
184  res.nrow = internal::convert_index<int>(mat.rows());
185  res.ncol = internal::convert_index<int>(mat.cols());
186 
187  res.storage.lda = internal::convert_index<int>(MatrixType::IsVectorAtCompileTime ? mat.size() : mat.outerStride());
188  res.storage.values = (void*)(mat.data());
189  return res;
190  }
191 
192  template<typename MatrixType>
194  {
195  MatrixType &mat(a_mat.derived());
196  SluMatrix res;
197  if ((MatrixType::Flags&RowMajorBit)==RowMajorBit)
198  {
199  res.setStorageType(SLU_NR);
200  res.nrow = internal::convert_index<int>(mat.cols());
201  res.ncol = internal::convert_index<int>(mat.rows());
202  }
203  else
204  {
205  res.setStorageType(SLU_NC);
206  res.nrow = internal::convert_index<int>(mat.rows());
207  res.ncol = internal::convert_index<int>(mat.cols());
208  }
209 
210  res.Mtype = SLU_GE;
211 
212  res.storage.nnz = internal::convert_index<int>(mat.nonZeros());
213  res.storage.values = mat.valuePtr();
214  res.storage.innerInd = mat.innerIndexPtr();
215  res.storage.outerInd = mat.outerIndexPtr();
216 
217  res.setScalarType<typename MatrixType::Scalar>();
218 
219  // FIXME the following is not very accurate
220  if (int(MatrixType::Flags) & int(Upper))
221  res.Mtype = SLU_TRU;
222  if (int(MatrixType::Flags) & int(Lower))
223  res.Mtype = SLU_TRL;
224 
225  eigen_assert(((int(MatrixType::Flags) & int(SelfAdjoint))==0) && "SelfAdjoint matrix shape not supported by SuperLU");
226 
227  return res;
228  }
229 };
230 
231 template<typename Scalar, int Rows, int Cols, int Options, int MRows, int MCols>
232 struct SluMatrixMapHelper<Matrix<Scalar,Rows,Cols,Options,MRows,MCols> >
233 {
235  static void run(MatrixType& mat, SluMatrix& res)
236  {
237  eigen_assert( ((Options&RowMajor)!=RowMajor) && "row-major dense matrices is not supported by SuperLU");
238  res.setStorageType(SLU_DN);
239  res.setScalarType<Scalar>();
240  res.Mtype = SLU_GE;
241 
242  res.nrow = mat.rows();
243  res.ncol = mat.cols();
244 
245  res.storage.lda = mat.outerStride();
246  res.storage.values = mat.data();
247  }
248 };
249 
250 template<typename Derived>
252 {
253  typedef Derived MatrixType;
254  static void run(MatrixType& mat, SluMatrix& res)
255  {
256  if ((MatrixType::Flags&RowMajorBit)==RowMajorBit)
257  {
258  res.setStorageType(SLU_NR);
259  res.nrow = mat.cols();
260  res.ncol = mat.rows();
261  }
262  else
263  {
264  res.setStorageType(SLU_NC);
265  res.nrow = mat.rows();
266  res.ncol = mat.cols();
267  }
268 
269  res.Mtype = SLU_GE;
270 
271  res.storage.nnz = mat.nonZeros();
272  res.storage.values = mat.valuePtr();
273  res.storage.innerInd = mat.innerIndexPtr();
274  res.storage.outerInd = mat.outerIndexPtr();
275 
276  res.setScalarType<typename MatrixType::Scalar>();
277 
278  // FIXME the following is not very accurate
279  if (MatrixType::Flags & Upper)
280  res.Mtype = SLU_TRU;
281  if (MatrixType::Flags & Lower)
282  res.Mtype = SLU_TRL;
283 
284  eigen_assert(((MatrixType::Flags & SelfAdjoint)==0) && "SelfAdjoint matrix shape not supported by SuperLU");
285  }
286 };
287 
288 namespace internal {
289 
290 template<typename MatrixType>
292 {
293  return SluMatrix::Map(mat);
294 }
295 
297 template<typename Scalar, int Flags, typename Index>
299 {
300  eigen_assert(((Flags&RowMajor)==RowMajor && sluMat.Stype == SLU_NR)
301  || ((Flags&ColMajor)==ColMajor && sluMat.Stype == SLU_NC));
302 
303  Index outerSize = (Flags&RowMajor)==RowMajor ? sluMat.ncol : sluMat.nrow;
304 
306  sluMat.nrow, sluMat.ncol, sluMat.storage.outerInd[outerSize],
307  sluMat.storage.outerInd, sluMat.storage.innerInd, reinterpret_cast<Scalar*>(sluMat.storage.values) );
308 }
309 
310 } // end namespace internal
311 
316 template<typename _MatrixType, typename Derived>
317 class SuperLUBase : public SparseSolverBase<Derived>
318 {
319  protected:
321  using Base::derived;
322  using Base::m_isInitialized;
323  public:
324  typedef _MatrixType MatrixType;
325  typedef typename MatrixType::Scalar Scalar;
327  typedef typename MatrixType::StorageIndex StorageIndex;
333  enum {
334  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
335  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
336  };
337 
338  public:
339 
341 
343  {
344  clearFactors();
345  }
346 
347  inline Index rows() const { return m_matrix.rows(); }
348  inline Index cols() const { return m_matrix.cols(); }
349 
351  inline superlu_options_t& options() { return m_sluOptions; }
352 
359  {
360  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
361  return m_info;
362  }
363 
365  void compute(const MatrixType& matrix)
366  {
367  derived().analyzePattern(matrix);
368  derived().factorize(matrix);
369  }
370 
377  void analyzePattern(const MatrixType& /*matrix*/)
378  {
379  m_isInitialized = true;
380  m_info = Success;
381  m_analysisIsOk = true;
382  m_factorizationIsOk = false;
383  }
384 
385  template<typename Stream>
386  void dumpMemory(Stream& /*s*/)
387  {}
388 
389  protected:
390 
391  void initFactorization(const MatrixType& a)
392  {
393  set_default_options(&this->m_sluOptions);
394 
395  const Index size = a.rows();
396  m_matrix = a;
397 
398  m_sluA = internal::asSluMatrix(m_matrix);
399  clearFactors();
400 
401  m_p.resize(size);
402  m_q.resize(size);
403  m_sluRscale.resize(size);
404  m_sluCscale.resize(size);
405  m_sluEtree.resize(size);
406 
407  // set empty B and X
408  m_sluB.setStorageType(SLU_DN);
409  m_sluB.setScalarType<Scalar>();
410  m_sluB.Mtype = SLU_GE;
411  m_sluB.storage.values = 0;
412  m_sluB.nrow = 0;
413  m_sluB.ncol = 0;
414  m_sluB.storage.lda = internal::convert_index<int>(size);
415  m_sluX = m_sluB;
416 
417  m_extractedDataAreDirty = true;
418  }
419 
420  void init()
421  {
422  m_info = InvalidInput;
423  m_isInitialized = false;
424  m_sluL.Store = 0;
425  m_sluU.Store = 0;
426  }
427 
428  void extractData() const;
429 
431  {
432  if(m_sluL.Store)
433  Destroy_SuperNode_Matrix(&m_sluL);
434  if(m_sluU.Store)
435  Destroy_CompCol_Matrix(&m_sluU);
436 
437  m_sluL.Store = 0;
438  m_sluU.Store = 0;
439 
440  memset(&m_sluL,0,sizeof m_sluL);
441  memset(&m_sluU,0,sizeof m_sluU);
442  }
443 
444  // cached data to reduce reallocation, etc.
445  mutable LUMatrixType m_l;
446  mutable LUMatrixType m_u;
447  mutable IntColVectorType m_p;
448  mutable IntRowVectorType m_q;
449 
450  mutable LUMatrixType m_matrix; // copy of the factorized matrix
451  mutable SluMatrix m_sluA;
452  mutable SuperMatrix m_sluL, m_sluU;
453  mutable SluMatrix m_sluB, m_sluX;
454  mutable SuperLUStat_t m_sluStat;
455  mutable superlu_options_t m_sluOptions;
456  mutable std::vector<int> m_sluEtree;
459  mutable char m_sluEqued;
460 
465 
466  private:
468 };
469 
470 
487 template<typename _MatrixType>
488 class SuperLU : public SuperLUBase<_MatrixType,SuperLU<_MatrixType> >
489 {
490  public:
492  typedef _MatrixType MatrixType;
493  typedef typename Base::Scalar Scalar;
494  typedef typename Base::RealScalar RealScalar;
502 
503  public:
504  using Base::_solve_impl;
505 
506  SuperLU() : Base() { init(); }
507 
508  explicit SuperLU(const MatrixType& matrix) : Base()
509  {
510  init();
511  Base::compute(matrix);
512  }
513 
515  {
516  }
517 
524  void analyzePattern(const MatrixType& matrix)
525  {
526  m_info = InvalidInput;
527  m_isInitialized = false;
528  Base::analyzePattern(matrix);
529  }
530 
537  void factorize(const MatrixType& matrix);
538 
540  template<typename Rhs,typename Dest>
541  void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const;
542 
543  inline const LMatrixType& matrixL() const
544  {
545  if (m_extractedDataAreDirty) this->extractData();
546  return m_l;
547  }
548 
549  inline const UMatrixType& matrixU() const
550  {
551  if (m_extractedDataAreDirty) this->extractData();
552  return m_u;
553  }
554 
555  inline const IntColVectorType& permutationP() const
556  {
557  if (m_extractedDataAreDirty) this->extractData();
558  return m_p;
559  }
560 
561  inline const IntRowVectorType& permutationQ() const
562  {
563  if (m_extractedDataAreDirty) this->extractData();
564  return m_q;
565  }
566 
567  Scalar determinant() const;
568 
569  protected:
570 
571  using Base::m_matrix;
572  using Base::m_sluOptions;
573  using Base::m_sluA;
574  using Base::m_sluB;
575  using Base::m_sluX;
576  using Base::m_p;
577  using Base::m_q;
578  using Base::m_sluEtree;
579  using Base::m_sluEqued;
580  using Base::m_sluRscale;
581  using Base::m_sluCscale;
582  using Base::m_sluL;
583  using Base::m_sluU;
584  using Base::m_sluStat;
585  using Base::m_sluFerr;
586  using Base::m_sluBerr;
587  using Base::m_l;
588  using Base::m_u;
589 
590  using Base::m_analysisIsOk;
591  using Base::m_factorizationIsOk;
592  using Base::m_extractedDataAreDirty;
593  using Base::m_isInitialized;
594  using Base::m_info;
595 
596  void init()
597  {
598  Base::init();
599 
600  set_default_options(&this->m_sluOptions);
601  m_sluOptions.PrintStat = NO;
602  m_sluOptions.ConditionNumber = NO;
603  m_sluOptions.Trans = NOTRANS;
604  m_sluOptions.ColPerm = COLAMD;
605  }
606 
607 
608  private:
610 };
611 
612 template<typename MatrixType>
614 {
615  eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
616  if(!m_analysisIsOk)
617  {
618  m_info = InvalidInput;
619  return;
620  }
621 
622  this->initFactorization(a);
623 
624  m_sluOptions.ColPerm = COLAMD;
625  int info = 0;
626  RealScalar recip_pivot_growth, rcond;
627  RealScalar ferr, berr;
628 
629  StatInit(&m_sluStat);
630  SuperLU_gssvx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0],
631  &m_sluEqued, &m_sluRscale[0], &m_sluCscale[0],
632  &m_sluL, &m_sluU,
633  NULL, 0,
634  &m_sluB, &m_sluX,
635  &recip_pivot_growth, &rcond,
636  &ferr, &berr,
637  &m_sluStat, &info, Scalar());
638  StatFree(&m_sluStat);
639 
640  m_extractedDataAreDirty = true;
641 
642  // FIXME how to better check for errors ???
643  m_info = info == 0 ? Success : NumericalIssue;
644  m_factorizationIsOk = true;
645 }
646 
647 template<typename MatrixType>
648 template<typename Rhs,typename Dest>
650 {
651  eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or analyzePattern()/factorize()");
652 
653  const Index rhsCols = b.cols();
654  eigen_assert(m_matrix.rows()==b.rows());
655 
656  m_sluOptions.Trans = NOTRANS;
657  m_sluOptions.Fact = FACTORED;
658  m_sluOptions.IterRefine = NOREFINE;
659 
660 
661  m_sluFerr.resize(rhsCols);
662  m_sluBerr.resize(rhsCols);
663 
666 
667  m_sluB = SluMatrix::Map(b_ref.const_cast_derived());
668  m_sluX = SluMatrix::Map(x_ref.const_cast_derived());
669 
670  typename Rhs::PlainObject b_cpy;
671  if(m_sluEqued!='N')
672  {
673  b_cpy = b;
674  m_sluB = SluMatrix::Map(b_cpy.const_cast_derived());
675  }
676 
677  StatInit(&m_sluStat);
678  int info = 0;
679  RealScalar recip_pivot_growth, rcond;
680  SuperLU_gssvx(&m_sluOptions, &m_sluA,
681  m_q.data(), m_p.data(),
682  &m_sluEtree[0], &m_sluEqued,
683  &m_sluRscale[0], &m_sluCscale[0],
684  &m_sluL, &m_sluU,
685  NULL, 0,
686  &m_sluB, &m_sluX,
687  &recip_pivot_growth, &rcond,
688  &m_sluFerr[0], &m_sluBerr[0],
689  &m_sluStat, &info, Scalar());
690  StatFree(&m_sluStat);
691 
692  if(x.derived().data() != x_ref.data())
693  x = x_ref;
694 
695  m_info = info==0 ? Success : NumericalIssue;
696 }
697 
698 // the code of this extractData() function has been adapted from the SuperLU's Matlab support code,
699 //
700 // Copyright (c) 1994 by Xerox Corporation. All rights reserved.
701 //
702 // THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
703 // EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK.
704 //
705 template<typename MatrixType, typename Derived>
707 {
708  eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for extracting factors, you must first call either compute() or analyzePattern()/factorize()");
709  if (m_extractedDataAreDirty)
710  {
711  int upper;
712  int fsupc, istart, nsupr;
713  int lastl = 0, lastu = 0;
714  SCformat *Lstore = static_cast<SCformat*>(m_sluL.Store);
715  NCformat *Ustore = static_cast<NCformat*>(m_sluU.Store);
716  Scalar *SNptr;
717 
718  const Index size = m_matrix.rows();
719  m_l.resize(size,size);
720  m_l.resizeNonZeros(Lstore->nnz);
721  m_u.resize(size,size);
722  m_u.resizeNonZeros(Ustore->nnz);
723 
724  int* Lcol = m_l.outerIndexPtr();
725  int* Lrow = m_l.innerIndexPtr();
726  Scalar* Lval = m_l.valuePtr();
727 
728  int* Ucol = m_u.outerIndexPtr();
729  int* Urow = m_u.innerIndexPtr();
730  Scalar* Uval = m_u.valuePtr();
731 
732  Ucol[0] = 0;
733  Ucol[0] = 0;
734 
735  /* for each supernode */
736  for (int k = 0; k <= Lstore->nsuper; ++k)
737  {
738  fsupc = L_FST_SUPC(k);
739  istart = L_SUB_START(fsupc);
740  nsupr = L_SUB_START(fsupc+1) - istart;
741  upper = 1;
742 
743  /* for each column in the supernode */
744  for (int j = fsupc; j < L_FST_SUPC(k+1); ++j)
745  {
746  SNptr = &((Scalar*)Lstore->nzval)[L_NZ_START(j)];
747 
748  /* Extract U */
749  for (int i = U_NZ_START(j); i < U_NZ_START(j+1); ++i)
750  {
751  Uval[lastu] = ((Scalar*)Ustore->nzval)[i];
752  /* Matlab doesn't like explicit zero. */
753  if (Uval[lastu] != 0.0)
754  Urow[lastu++] = U_SUB(i);
755  }
756  for (int i = 0; i < upper; ++i)
757  {
758  /* upper triangle in the supernode */
759  Uval[lastu] = SNptr[i];
760  /* Matlab doesn't like explicit zero. */
761  if (Uval[lastu] != 0.0)
762  Urow[lastu++] = L_SUB(istart+i);
763  }
764  Ucol[j+1] = lastu;
765 
766  /* Extract L */
767  Lval[lastl] = 1.0; /* unit diagonal */
768  Lrow[lastl++] = L_SUB(istart + upper - 1);
769  for (int i = upper; i < nsupr; ++i)
770  {
771  Lval[lastl] = SNptr[i];
772  /* Matlab doesn't like explicit zero. */
773  if (Lval[lastl] != 0.0)
774  Lrow[lastl++] = L_SUB(istart+i);
775  }
776  Lcol[j+1] = lastl;
777 
778  ++upper;
779  } /* for j ... */
780 
781  } /* for k ... */
782 
783  // squeeze the matrices :
784  m_l.resizeNonZeros(lastl);
785  m_u.resizeNonZeros(lastu);
786 
787  m_extractedDataAreDirty = false;
788  }
789 }
790 
791 template<typename MatrixType>
793 {
794  eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for computing the determinant, you must first call either compute() or analyzePattern()/factorize()");
795 
796  if (m_extractedDataAreDirty)
797  this->extractData();
798 
799  Scalar det = Scalar(1);
800  for (int j=0; j<m_u.cols(); ++j)
801  {
802  if (m_u.outerIndexPtr()[j+1]-m_u.outerIndexPtr()[j] > 0)
803  {
804  int lastId = m_u.outerIndexPtr()[j+1]-1;
805  eigen_assert(m_u.innerIndexPtr()[lastId]<=j);
806  if (m_u.innerIndexPtr()[lastId]==j)
807  det *= m_u.valuePtr()[lastId];
808  }
809  }
810  if(PermutationMap(m_p.data(),m_p.size()).determinant()*PermutationMap(m_q.data(),m_q.size()).determinant()<0)
811  det = -det;
812  if(m_sluEqued!='N')
813  return det/m_sluRscale.prod()/m_sluCscale.prod();
814  else
815  return det;
816 }
817 
818 #ifdef EIGEN_PARSED_BY_DOXYGEN
819 #define EIGEN_SUPERLU_HAS_ILU
820 #endif
821 
822 #ifdef EIGEN_SUPERLU_HAS_ILU
823 
840 template<typename _MatrixType>
841 class SuperILU : public SuperLUBase<_MatrixType,SuperILU<_MatrixType> >
842 {
843  public:
845  typedef _MatrixType MatrixType;
846  typedef typename Base::Scalar Scalar;
847  typedef typename Base::RealScalar RealScalar;
848 
849  public:
850  using Base::_solve_impl;
851 
852  SuperILU() : Base() { init(); }
853 
854  SuperILU(const MatrixType& matrix) : Base()
855  {
856  init();
857  Base::compute(matrix);
858  }
859 
860  ~SuperILU()
861  {
862  }
863 
870  void analyzePattern(const MatrixType& matrix)
871  {
872  Base::analyzePattern(matrix);
873  }
874 
881  void factorize(const MatrixType& matrix);
882 
883  #ifndef EIGEN_PARSED_BY_DOXYGEN
884 
885  template<typename Rhs,typename Dest>
886  void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const;
887  #endif // EIGEN_PARSED_BY_DOXYGEN
888 
889  protected:
890 
891  using Base::m_matrix;
892  using Base::m_sluOptions;
893  using Base::m_sluA;
894  using Base::m_sluB;
895  using Base::m_sluX;
896  using Base::m_p;
897  using Base::m_q;
898  using Base::m_sluEtree;
899  using Base::m_sluEqued;
900  using Base::m_sluRscale;
901  using Base::m_sluCscale;
902  using Base::m_sluL;
903  using Base::m_sluU;
904  using Base::m_sluStat;
905  using Base::m_sluFerr;
906  using Base::m_sluBerr;
907  using Base::m_l;
908  using Base::m_u;
909 
910  using Base::m_analysisIsOk;
911  using Base::m_factorizationIsOk;
912  using Base::m_extractedDataAreDirty;
913  using Base::m_isInitialized;
914  using Base::m_info;
915 
916  void init()
917  {
918  Base::init();
919 
920  ilu_set_default_options(&m_sluOptions);
921  m_sluOptions.PrintStat = NO;
922  m_sluOptions.ConditionNumber = NO;
923  m_sluOptions.Trans = NOTRANS;
924  m_sluOptions.ColPerm = MMD_AT_PLUS_A;
925 
926  // no attempt to preserve column sum
927  m_sluOptions.ILU_MILU = SILU;
928  // only basic ILU(k) support -- no direct control over memory consumption
929  // better to use ILU_DropRule = DROP_BASIC | DROP_AREA
930  // and set ILU_FillFactor to max memory growth
931  m_sluOptions.ILU_DropRule = DROP_BASIC;
932  m_sluOptions.ILU_DropTol = NumTraits<Scalar>::dummy_precision()*10;
933  }
934 
935  private:
936  SuperILU(SuperILU& ) { }
937 };
938 
939 template<typename MatrixType>
940 void SuperILU<MatrixType>::factorize(const MatrixType& a)
941 {
942  eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
943  if(!m_analysisIsOk)
944  {
945  m_info = InvalidInput;
946  return;
947  }
948 
949  this->initFactorization(a);
950 
951  int info = 0;
952  RealScalar recip_pivot_growth, rcond;
953 
954  StatInit(&m_sluStat);
955  SuperLU_gsisx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0],
956  &m_sluEqued, &m_sluRscale[0], &m_sluCscale[0],
957  &m_sluL, &m_sluU,
958  NULL, 0,
959  &m_sluB, &m_sluX,
960  &recip_pivot_growth, &rcond,
961  &m_sluStat, &info, Scalar());
962  StatFree(&m_sluStat);
963 
964  // FIXME how to better check for errors ???
965  m_info = info == 0 ? Success : NumericalIssue;
966  m_factorizationIsOk = true;
967 }
968 
969 #ifndef EIGEN_PARSED_BY_DOXYGEN
970 template<typename MatrixType>
971 template<typename Rhs,typename Dest>
972 void SuperILU<MatrixType>::_solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest>& x) const
973 {
974  eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or analyzePattern()/factorize()");
975 
976  const int rhsCols = b.cols();
977  eigen_assert(m_matrix.rows()==b.rows());
978 
979  m_sluOptions.Trans = NOTRANS;
980  m_sluOptions.Fact = FACTORED;
981  m_sluOptions.IterRefine = NOREFINE;
982 
983  m_sluFerr.resize(rhsCols);
984  m_sluBerr.resize(rhsCols);
985 
988 
989  m_sluB = SluMatrix::Map(b_ref.const_cast_derived());
990  m_sluX = SluMatrix::Map(x_ref.const_cast_derived());
991 
992  typename Rhs::PlainObject b_cpy;
993  if(m_sluEqued!='N')
994  {
995  b_cpy = b;
996  m_sluB = SluMatrix::Map(b_cpy.const_cast_derived());
997  }
998 
999  int info = 0;
1000  RealScalar recip_pivot_growth, rcond;
1001 
1002  StatInit(&m_sluStat);
1003  SuperLU_gsisx(&m_sluOptions, &m_sluA,
1004  m_q.data(), m_p.data(),
1005  &m_sluEtree[0], &m_sluEqued,
1006  &m_sluRscale[0], &m_sluCscale[0],
1007  &m_sluL, &m_sluU,
1008  NULL, 0,
1009  &m_sluB, &m_sluX,
1010  &recip_pivot_growth, &rcond,
1011  &m_sluStat, &info, Scalar());
1012  StatFree(&m_sluStat);
1013 
1014  if(x.derived().data() != x_ref.data())
1015  x = x_ref;
1016 
1017  m_info = info==0 ? Success : NumericalIssue;
1018 }
1019 #endif
1020 
1021 #endif
1022 
1023 } // end namespace Eigen
1024 
1025 #endif // EIGEN_SUPERLUSUPPORT_H
MappedSparseMatrix< Scalar, Flags, Index > map_superlu(SluMatrix &sluMat)
SCALAR Scalar
Definition: bench_gemm.cpp:46
const UMatrixType & matrixU() const
SluMatrix asSluMatrix(MatrixType &mat)
SparseSolverBase< Derived > Base
void compute(const MatrixType &matrix)
Matrix< int, 1, MatrixType::ColsAtCompileTime > IntRowVectorType
Scalar * b
Definition: benchVecAdd.cpp:17
const IntRowVectorType & permutationQ() const
Base::PermutationMap PermutationMap
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar * data() const
A sparse direct LU factorization and solver based on the SuperLU library.
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:94
MatrixType::Scalar Scalar
void determinant(const MatrixType &m)
Definition: determinant.cpp:14
static SluMatrix Map(SparseMatrixBase< MatrixType > &a_mat)
A base class for sparse solvers.
Scalar Scalar * c
Definition: benchVecAdd.cpp:17
#define DECL_GSSVX(PREFIX, FLOATTYPE, KEYTYPE)
Index rows() const
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
void _solve_impl(const MatrixBase< Rhs > &b, MatrixBase< Dest > &dest) const
MatrixXf MatrixType
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:232
if((m *x).isApprox(y))
struct Eigen::SluMatrix::@895 storage
EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR Index outerStride() const EIGEN_NOEXCEPT
MatrixType::StorageIndex StorageIndex
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
const unsigned int RowMajorBit
Definition: Constants.h:66
Base::Scalar Scalar
ComputationInfo m_info
void analyzePattern(const MatrixType &)
LUMatrixType m_matrix
Base::IntColVectorType IntColVectorType
else if n * info
cout<< "Here is the matrix m:"<< endl<< m<< endl;Matrix< ptrdiff_t, 3, 1 > res
Scalar Scalar int size
Definition: benchVecAdd.cpp:17
SuperLU(const MatrixType &matrix)
Index cols() const
Base class of any sparse matrices or sparse expressions.
Base::RealScalar RealScalar
SluMatrix & operator=(const SluMatrix &other)
superlu_options_t m_sluOptions
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
const IntColVectorType & permutationP() const
#define eigen_assert(x)
Definition: Macros.h:1037
void dumpMemory(Stream &)
The base class for the direct and incomplete LU factorization of SuperLU.
SuperLUBase(SuperLUBase &)
RealScalar s
Matrix< RealScalar, Dynamic, 1 > m_sluRscale
Scalar determinant() const
#define NULL
Definition: ccolamd.c:609
SuperLU(SuperLU &)
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:47
std::vector< int > m_sluEtree
TriangularView< LUMatrixType, Lower|UnitDiag > LMatrixType
void initFactorization(const MatrixType &a)
Map< PermutationMatrix< Dynamic, Dynamic, int > > PermutationMap
static const int Cols
Base::LUMatrixType LUMatrixType
A matrix or vector expression mapping an existing expression.
Definition: Ref.h:281
IntRowVectorType m_q
TriangularView< LUMatrixType, Upper > UMatrixType
SuperLUStat_t m_sluStat
void extractData() const
Matrix< int, MatrixType::RowsAtCompileTime, 1 > IntColVectorType
detail::initimpl::constructor< Args... > init()
Binds an existing constructor taking arguments Args...
Definition: pybind11.h:1882
Matrix< RealScalar, Dynamic, 1 > m_sluFerr
SuperLUBase< _MatrixType, SuperLU > Base
Expression of a triangular part in a matrix.
SparseMatrix< Scalar > LUMatrixType
_MatrixType MatrixType
MatrixType::RealScalar RealScalar
const Derived & derived() const
void setStorageType(Stype_t t)
_MatrixType MatrixType
static SluMatrix Map(MatrixBase< MatrixType > &_mat)
void factorize(const MatrixType &matrix)
EIGEN_DONT_INLINE void compute(Solver &solver, const MatrixType &A)
superlu_options_t & options()
Map< Matrix< T, Dynamic, Dynamic, ColMajor >, 0, OuterStride<> > matrix(T *data, int rows, int cols, int stride)
The matrix class, also used for vectors and row-vectors.
Base::StorageIndex StorageIndex
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy x
SluMatrix(const SluMatrix &other)
ComputationInfo
Definition: Constants.h:440
Matrix< Scalar, Dynamic, 1 > Vector
IntColVectorType m_p
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
const LMatrixType & matrixL() const
Base::IntRowVectorType IntRowVectorType
void analyzePattern(const MatrixType &matrix)
std::ptrdiff_t j
Point2 t(10, 10)
ComputationInfo info() const
Reports whether previous computation was successful.
static void run(MatrixType &mat, SluMatrix &res)


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:36:28