Path Tracer
SuperLUSupport.h
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>
106 struct SluMatrixMapHelper;
107 
115 struct SluMatrix : SuperMatrix
116 {
117  SluMatrix()
118  {
119  Store = &storage;
120  }
121 
122  SluMatrix(const SluMatrix& other)
123  : SuperMatrix(other)
124  {
125  Store = &storage;
126  storage = other.storage;
127  }
128 
129  SluMatrix& operator=(const SluMatrix& other)
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>
158  void setScalarType()
159  {
160  if (internal::is_same<Scalar,float>::value)
161  Dtype = SLU_S;
162  else if (internal::is_same<Scalar,double>::value)
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>
175  static SluMatrix Map(MatrixBase<MatrixType>& _mat)
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>
193  static SluMatrix Map(SparseMatrixBase<MatrixType>& a_mat)
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 (MatrixType::Flags & Upper)
221  res.Mtype = SLU_TRU;
222  if (MatrixType::Flags & Lower)
223  res.Mtype = SLU_TRL;
224 
225  eigen_assert(((MatrixType::Flags & 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 {
234  typedef Matrix<Scalar,Rows,Cols,Options,MRows,MCols> MatrixType;
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>
251 struct SluMatrixMapHelper<SparseMatrixBase<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>
291 SluMatrix asSluMatrix(MatrixType& mat)
292 {
293  return SluMatrix::Map(mat);
294 }
295 
297 template<typename Scalar, int Flags, typename Index>
298 MappedSparseMatrix<Scalar,Flags,Index> map_superlu(SluMatrix& sluMat)
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 
305  return MappedSparseMatrix<Scalar,Flags,Index>(
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:
320  typedef SparseSolverBase<Derived> Base;
321  using Base::derived;
322  using Base::m_isInitialized;
323  public:
324  typedef _MatrixType MatrixType;
325  typedef typename MatrixType::Scalar Scalar;
326  typedef typename MatrixType::RealScalar RealScalar;
327  typedef typename MatrixType::StorageIndex StorageIndex;
328  typedef Matrix<Scalar,Dynamic,1> Vector;
329  typedef Matrix<int, 1, MatrixType::ColsAtCompileTime> IntRowVectorType;
330  typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType;
331  typedef Map<PermutationMatrix<Dynamic,Dynamic,int> > PermutationMap;
332  typedef SparseMatrix<Scalar> LUMatrixType;
333  enum {
334  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
335  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
336  };
337 
338  public:
339 
340  SuperLUBase() {}
341 
342  ~SuperLUBase()
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 
358  ComputationInfo info() const
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 
430  void clearFactors()
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;
457  mutable Matrix<RealScalar,Dynamic,1> m_sluRscale, m_sluCscale;
458  mutable Matrix<RealScalar,Dynamic,1> m_sluFerr, m_sluBerr;
459  mutable char m_sluEqued;
460 
461  mutable ComputationInfo m_info;
462  int m_factorizationIsOk;
463  int m_analysisIsOk;
464  mutable bool m_extractedDataAreDirty;
465 
466  private:
467  SuperLUBase(SuperLUBase& ) { }
468 };
469 
470 
487 template<typename _MatrixType>
488 class SuperLU : public SuperLUBase<_MatrixType,SuperLU<_MatrixType> >
489 {
490  public:
491  typedef SuperLUBase<_MatrixType,SuperLU> Base;
492  typedef _MatrixType MatrixType;
493  typedef typename Base::Scalar Scalar;
494  typedef typename Base::RealScalar RealScalar;
495  typedef typename Base::StorageIndex StorageIndex;
496  typedef typename Base::IntRowVectorType IntRowVectorType;
497  typedef typename Base::IntColVectorType IntColVectorType;
498  typedef typename Base::PermutationMap PermutationMap;
499  typedef typename Base::LUMatrixType LUMatrixType;
500  typedef TriangularView<LUMatrixType, Lower|UnitDiag> LMatrixType;
501  typedef TriangularView<LUMatrixType, Upper> UMatrixType;
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 
514  ~SuperLU()
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:
609  SuperLU(SuperLU& ) { }
610 };
611 
612 template<typename MatrixType>
613 void SuperLU<MatrixType>::factorize(const MatrixType& a)
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>
649 void SuperLU<MatrixType>::_solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest>& x) const
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 size = m_matrix.rows();
654  const Index rhsCols = b.cols();
655  eigen_assert(size==b.rows());
656 
657  m_sluOptions.Trans = NOTRANS;
658  m_sluOptions.Fact = FACTORED;
659  m_sluOptions.IterRefine = NOREFINE;
660 
661 
662  m_sluFerr.resize(rhsCols);
663  m_sluBerr.resize(rhsCols);
664 
665  Ref<const Matrix<typename Rhs::Scalar,Dynamic,Dynamic,ColMajor> > b_ref(b);
666  Ref<const Matrix<typename Dest::Scalar,Dynamic,Dynamic,ColMajor> > x_ref(x);
667 
668  m_sluB = SluMatrix::Map(b_ref.const_cast_derived());
669  m_sluX = SluMatrix::Map(x_ref.const_cast_derived());
670 
671  typename Rhs::PlainObject b_cpy;
672  if(m_sluEqued!='N')
673  {
674  b_cpy = b;
675  m_sluB = SluMatrix::Map(b_cpy.const_cast_derived());
676  }
677 
678  StatInit(&m_sluStat);
679  int info = 0;
680  RealScalar recip_pivot_growth, rcond;
681  SuperLU_gssvx(&m_sluOptions, &m_sluA,
682  m_q.data(), m_p.data(),
683  &m_sluEtree[0], &m_sluEqued,
684  &m_sluRscale[0], &m_sluCscale[0],
685  &m_sluL, &m_sluU,
686  NULL, 0,
687  &m_sluB, &m_sluX,
688  &recip_pivot_growth, &rcond,
689  &m_sluFerr[0], &m_sluBerr[0],
690  &m_sluStat, &info, Scalar());
691  StatFree(&m_sluStat);
692 
693  if(x.derived().data() != x_ref.data())
694  x = x_ref;
695 
696  m_info = info==0 ? Success : NumericalIssue;
697 }
698 
699 // the code of this extractData() function has been adapted from the SuperLU's Matlab support code,
700 //
701 // Copyright (c) 1994 by Xerox Corporation. All rights reserved.
702 //
703 // THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
704 // EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK.
705 //
706 template<typename MatrixType, typename Derived>
707 void SuperLUBase<MatrixType,Derived>::extractData() const
708 {
709  eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for extracting factors, you must first call either compute() or analyzePattern()/factorize()");
710  if (m_extractedDataAreDirty)
711  {
712  int upper;
713  int fsupc, istart, nsupr;
714  int lastl = 0, lastu = 0;
715  SCformat *Lstore = static_cast<SCformat*>(m_sluL.Store);
716  NCformat *Ustore = static_cast<NCformat*>(m_sluU.Store);
717  Scalar *SNptr;
718 
719  const Index size = m_matrix.rows();
720  m_l.resize(size,size);
721  m_l.resizeNonZeros(Lstore->nnz);
722  m_u.resize(size,size);
723  m_u.resizeNonZeros(Ustore->nnz);
724 
725  int* Lcol = m_l.outerIndexPtr();
726  int* Lrow = m_l.innerIndexPtr();
727  Scalar* Lval = m_l.valuePtr();
728 
729  int* Ucol = m_u.outerIndexPtr();
730  int* Urow = m_u.innerIndexPtr();
731  Scalar* Uval = m_u.valuePtr();
732 
733  Ucol[0] = 0;
734  Ucol[0] = 0;
735 
736  /* for each supernode */
737  for (int k = 0; k <= Lstore->nsuper; ++k)
738  {
739  fsupc = L_FST_SUPC(k);
740  istart = L_SUB_START(fsupc);
741  nsupr = L_SUB_START(fsupc+1) - istart;
742  upper = 1;
743 
744  /* for each column in the supernode */
745  for (int j = fsupc; j < L_FST_SUPC(k+1); ++j)
746  {
747  SNptr = &((Scalar*)Lstore->nzval)[L_NZ_START(j)];
748 
749  /* Extract U */
750  for (int i = U_NZ_START(j); i < U_NZ_START(j+1); ++i)
751  {
752  Uval[lastu] = ((Scalar*)Ustore->nzval)[i];
753  /* Matlab doesn't like explicit zero. */
754  if (Uval[lastu] != 0.0)
755  Urow[lastu++] = U_SUB(i);
756  }
757  for (int i = 0; i < upper; ++i)
758  {
759  /* upper triangle in the supernode */
760  Uval[lastu] = SNptr[i];
761  /* Matlab doesn't like explicit zero. */
762  if (Uval[lastu] != 0.0)
763  Urow[lastu++] = L_SUB(istart+i);
764  }
765  Ucol[j+1] = lastu;
766 
767  /* Extract L */
768  Lval[lastl] = 1.0; /* unit diagonal */
769  Lrow[lastl++] = L_SUB(istart + upper - 1);
770  for (int i = upper; i < nsupr; ++i)
771  {
772  Lval[lastl] = SNptr[i];
773  /* Matlab doesn't like explicit zero. */
774  if (Lval[lastl] != 0.0)
775  Lrow[lastl++] = L_SUB(istart+i);
776  }
777  Lcol[j+1] = lastl;
778 
779  ++upper;
780  } /* for j ... */
781 
782  } /* for k ... */
783 
784  // squeeze the matrices :
785  m_l.resizeNonZeros(lastl);
786  m_u.resizeNonZeros(lastu);
787 
788  m_extractedDataAreDirty = false;
789  }
790 }
791 
792 template<typename MatrixType>
793 typename SuperLU<MatrixType>::Scalar SuperLU<MatrixType>::determinant() const
794 {
795  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()");
796 
797  if (m_extractedDataAreDirty)
798  this->extractData();
799 
800  Scalar det = Scalar(1);
801  for (int j=0; j<m_u.cols(); ++j)
802  {
803  if (m_u.outerIndexPtr()[j+1]-m_u.outerIndexPtr()[j] > 0)
804  {
805  int lastId = m_u.outerIndexPtr()[j+1]-1;
806  eigen_assert(m_u.innerIndexPtr()[lastId]<=j);
807  if (m_u.innerIndexPtr()[lastId]==j)
808  det *= m_u.valuePtr()[lastId];
809  }
810  }
811  if(PermutationMap(m_p.data(),m_p.size()).determinant()*PermutationMap(m_q.data(),m_q.size()).determinant()<0)
812  det = -det;
813  if(m_sluEqued!='N')
814  return det/m_sluRscale.prod()/m_sluCscale.prod();
815  else
816  return det;
817 }
818 
819 #ifdef EIGEN_PARSED_BY_DOXYGEN
820 #define EIGEN_SUPERLU_HAS_ILU
821 #endif
822 
823 #ifdef EIGEN_SUPERLU_HAS_ILU
824 
841 template<typename _MatrixType>
842 class SuperILU : public SuperLUBase<_MatrixType,SuperILU<_MatrixType> >
843 {
844  public:
845  typedef SuperLUBase<_MatrixType,SuperILU> Base;
846  typedef _MatrixType MatrixType;
847  typedef typename Base::Scalar Scalar;
848  typedef typename Base::RealScalar RealScalar;
849 
850  public:
851  using Base::_solve_impl;
852 
853  SuperILU() : Base() { init(); }
854 
855  SuperILU(const MatrixType& matrix) : Base()
856  {
857  init();
858  Base::compute(matrix);
859  }
860 
861  ~SuperILU()
862  {
863  }
864 
871  void analyzePattern(const MatrixType& matrix)
872  {
873  Base::analyzePattern(matrix);
874  }
875 
882  void factorize(const MatrixType& matrix);
883 
884  #ifndef EIGEN_PARSED_BY_DOXYGEN
885 
886  template<typename Rhs,typename Dest>
887  void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const;
888  #endif // EIGEN_PARSED_BY_DOXYGEN
889 
890  protected:
891 
892  using Base::m_matrix;
893  using Base::m_sluOptions;
894  using Base::m_sluA;
895  using Base::m_sluB;
896  using Base::m_sluX;
897  using Base::m_p;
898  using Base::m_q;
899  using Base::m_sluEtree;
900  using Base::m_sluEqued;
901  using Base::m_sluRscale;
902  using Base::m_sluCscale;
903  using Base::m_sluL;
904  using Base::m_sluU;
905  using Base::m_sluStat;
906  using Base::m_sluFerr;
907  using Base::m_sluBerr;
908  using Base::m_l;
909  using Base::m_u;
910 
911  using Base::m_analysisIsOk;
912  using Base::m_factorizationIsOk;
913  using Base::m_extractedDataAreDirty;
914  using Base::m_isInitialized;
915  using Base::m_info;
916 
917  void init()
918  {
919  Base::init();
920 
921  ilu_set_default_options(&m_sluOptions);
922  m_sluOptions.PrintStat = NO;
923  m_sluOptions.ConditionNumber = NO;
924  m_sluOptions.Trans = NOTRANS;
925  m_sluOptions.ColPerm = MMD_AT_PLUS_A;
926 
927  // no attempt to preserve column sum
928  m_sluOptions.ILU_MILU = SILU;
929  // only basic ILU(k) support -- no direct control over memory consumption
930  // better to use ILU_DropRule = DROP_BASIC | DROP_AREA
931  // and set ILU_FillFactor to max memory growth
932  m_sluOptions.ILU_DropRule = DROP_BASIC;
933  m_sluOptions.ILU_DropTol = NumTraits<Scalar>::dummy_precision()*10;
934  }
935 
936  private:
937  SuperILU(SuperILU& ) { }
938 };
939 
940 template<typename MatrixType>
941 void SuperILU<MatrixType>::factorize(const MatrixType& a)
942 {
943  eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
944  if(!m_analysisIsOk)
945  {
946  m_info = InvalidInput;
947  return;
948  }
949 
950  this->initFactorization(a);
951 
952  int info = 0;
953  RealScalar recip_pivot_growth, rcond;
954 
955  StatInit(&m_sluStat);
956  SuperLU_gsisx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0],
957  &m_sluEqued, &m_sluRscale[0], &m_sluCscale[0],
958  &m_sluL, &m_sluU,
959  NULL, 0,
960  &m_sluB, &m_sluX,
961  &recip_pivot_growth, &rcond,
962  &m_sluStat, &info, Scalar());
963  StatFree(&m_sluStat);
964 
965  // FIXME how to better check for errors ???
966  m_info = info == 0 ? Success : NumericalIssue;
967  m_factorizationIsOk = true;
968 }
969 
970 #ifndef EIGEN_PARSED_BY_DOXYGEN
971 template<typename MatrixType>
972 template<typename Rhs,typename Dest>
973 void SuperILU<MatrixType>::_solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest>& x) const
974 {
975  eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or analyzePattern()/factorize()");
976 
977  const int size = m_matrix.rows();
978  const int rhsCols = b.cols();
979  eigen_assert(size==b.rows());
980 
981  m_sluOptions.Trans = NOTRANS;
982  m_sluOptions.Fact = FACTORED;
983  m_sluOptions.IterRefine = NOREFINE;
984 
985  m_sluFerr.resize(rhsCols);
986  m_sluBerr.resize(rhsCols);
987 
988  Ref<const Matrix<typename Rhs::Scalar,Dynamic,Dynamic,ColMajor> > b_ref(b);
989  Ref<const Matrix<typename Dest::Scalar,Dynamic,Dynamic,ColMajor> > x_ref(x);
990 
991  m_sluB = SluMatrix::Map(b_ref.const_cast_derived());
992  m_sluX = SluMatrix::Map(x_ref.const_cast_derived());
993 
994  typename Rhs::PlainObject b_cpy;
995  if(m_sluEqued!='N')
996  {
997  b_cpy = b;
998  m_sluB = SluMatrix::Map(b_cpy.const_cast_derived());
999  }
1000 
1001  int info = 0;
1002  RealScalar recip_pivot_growth, rcond;
1003 
1004  StatInit(&m_sluStat);
1005  SuperLU_gsisx(&m_sluOptions, &m_sluA,
1006  m_q.data(), m_p.data(),
1007  &m_sluEtree[0], &m_sluEqued,
1008  &m_sluRscale[0], &m_sluCscale[0],
1009  &m_sluL, &m_sluU,
1010  NULL, 0,
1011  &m_sluB, &m_sluX,
1012  &recip_pivot_growth, &rcond,
1013  &m_sluStat, &info, Scalar());
1014  StatFree(&m_sluStat);
1015 
1016  if(x.derived().data() != x_ref.data())
1017  x = x_ref;
1018 
1019  m_info = info==0 ? Success : NumericalIssue;
1020 }
1021 #endif
1022 
1023 #endif
1024 
1025 } // end namespace Eigen
1026 
1027 #endif // EIGEN_SUPERLUSUPPORT_H
Eigen::SparseMatrix::cols
Index cols() const
Definition: SparseMatrix.h:140
Eigen::NumericalIssue
@ NumericalIssue
Definition: Constants.h:443
Eigen
Namespace containing all symbols from the Eigen library.
Definition: LDLT.h:16
Eigen::SparseMatrix< Scalar >
Eigen::RowMajorBit
const unsigned int RowMajorBit
Definition: Constants.h:65
Eigen::SuperLU::analyzePattern
void analyzePattern(const MatrixType &matrix)
Definition: SuperLUSupport.h:524
Eigen::Upper
@ Upper
Definition: Constants.h:210
Eigen::Success
@ Success
Definition: Constants.h:441
Eigen::RowMajor
@ RowMajor
Definition: Constants.h:320
Eigen::PlainObjectBase::resize
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void resize(Index rows, Index cols)
Definition: PlainObjectBase.h:279
Eigen::SuperLUBase::compute
void compute(const MatrixType &matrix)
Definition: SuperLUSupport.h:365
Eigen::SuperLUBase
The base class for the direct and incomplete LU factorization of SuperLU.
Definition: SuperLUSupport.h:318
Eigen::SuperLU
A sparse direct LU factorization and solver based on the SuperLU library.
Definition: SuperLUSupport.h:489
Eigen::SluMatrix
Definition: SuperLUSupport.h:116
Eigen::SuperLUBase::options
superlu_options_t & options()
Definition: SuperLUSupport.h:351
Eigen::SuperLUBase::analyzePattern
void analyzePattern(const MatrixType &)
Definition: SuperLUSupport.h:377
Eigen::SuperLU::factorize
void factorize(const MatrixType &matrix)
Definition: SuperLUSupport.h:613
Eigen::SluMatrixMapHelper
Definition: SuperLUSupport.h:106
Eigen::Lower
@ Lower
Definition: Constants.h:208
Eigen::Map
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:96
Eigen::SparseSolverBase
A base class for sparse solvers.
Definition: SparseSolverBase.h:68
Eigen::Ref
A matrix or vector expression mapping an existing expression.
Definition: Ref.h:197
Eigen::SparseMatrixBase
Base class of any sparse matrices or sparse expressions.
Definition: SparseMatrixBase.h:28
Eigen::internal::is_same
Definition: Meta.h:115
Eigen::Matrix
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:180
Eigen::InvalidInput
@ InvalidInput
Definition: Constants.h:448
Eigen::MatrixBase
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:50
Eigen::ColMajor
@ ColMajor
Definition: Constants.h:318
Eigen::ComputationInfo
ComputationInfo
Definition: Constants.h:439
Eigen::SparseMatrix::rows
Index rows() const
Definition: SparseMatrix.h:138
Eigen::TriangularView
Expression of a triangular part in a matrix.
Definition: TriangularMatrix.h:189
Eigen::SelfAdjoint
@ SelfAdjoint
Definition: Constants.h:224
Eigen::PlainObjectBase::data
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar * data() const
Definition: PlainObjectBase.h:255
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:42
Eigen::SuperLUBase::info
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: SuperLUSupport.h:358
Generated by   doxygen 1.8.19