Eigen  3.2.92
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 #define DECL_GSSVX(PREFIX,FLOATTYPE,KEYTYPE) \
16  extern "C" { \
17  extern void PREFIX##gssvx(superlu_options_t *, SuperMatrix *, int *, int *, int *, \
18  char *, FLOATTYPE *, FLOATTYPE *, SuperMatrix *, SuperMatrix *, \
19  void *, int, SuperMatrix *, SuperMatrix *, \
20  FLOATTYPE *, FLOATTYPE *, FLOATTYPE *, FLOATTYPE *, \
21  mem_usage_t *, SuperLUStat_t *, int *); \
22  } \
23  inline float SuperLU_gssvx(superlu_options_t *options, SuperMatrix *A, \
24  int *perm_c, int *perm_r, int *etree, char *equed, \
25  FLOATTYPE *R, FLOATTYPE *C, SuperMatrix *L, \
26  SuperMatrix *U, void *work, int lwork, \
27  SuperMatrix *B, SuperMatrix *X, \
28  FLOATTYPE *recip_pivot_growth, \
29  FLOATTYPE *rcond, FLOATTYPE *ferr, FLOATTYPE *berr, \
30  SuperLUStat_t *stats, int *info, KEYTYPE) { \
31  mem_usage_t mem_usage; \
32  PREFIX##gssvx(options, A, perm_c, perm_r, etree, equed, R, C, L, \
33  U, work, lwork, B, X, recip_pivot_growth, rcond, \
34  ferr, berr, &mem_usage, stats, info); \
35  return mem_usage.for_lu; /* bytes used by the factor storage */ \
36  }
37 
38 DECL_GSSVX(s,float,float)
39 DECL_GSSVX(c,float,std::complex<float>)
40 DECL_GSSVX(d,double,double)
41 DECL_GSSVX(z,double,std::complex<double>)
42 
43 #ifdef MILU_ALPHA
44 #define EIGEN_SUPERLU_HAS_ILU
45 #endif
46 
47 #ifdef EIGEN_SUPERLU_HAS_ILU
48 
49 // similarly for the incomplete factorization using gsisx
50 #define DECL_GSISX(PREFIX,FLOATTYPE,KEYTYPE) \
51  extern "C" { \
52  extern void PREFIX##gsisx(superlu_options_t *, SuperMatrix *, int *, int *, int *, \
53  char *, FLOATTYPE *, FLOATTYPE *, SuperMatrix *, SuperMatrix *, \
54  void *, int, SuperMatrix *, SuperMatrix *, FLOATTYPE *, FLOATTYPE *, \
55  mem_usage_t *, SuperLUStat_t *, int *); \
56  } \
57  inline float SuperLU_gsisx(superlu_options_t *options, SuperMatrix *A, \
58  int *perm_c, int *perm_r, int *etree, char *equed, \
59  FLOATTYPE *R, FLOATTYPE *C, SuperMatrix *L, \
60  SuperMatrix *U, void *work, int lwork, \
61  SuperMatrix *B, SuperMatrix *X, \
62  FLOATTYPE *recip_pivot_growth, \
63  FLOATTYPE *rcond, \
64  SuperLUStat_t *stats, int *info, KEYTYPE) { \
65  mem_usage_t mem_usage; \
66  PREFIX##gsisx(options, A, perm_c, perm_r, etree, equed, R, C, L, \
67  U, work, lwork, B, X, recip_pivot_growth, rcond, \
68  &mem_usage, stats, info); \
69  return mem_usage.for_lu; /* bytes used by the factor storage */ \
70  }
71 
72 DECL_GSISX(s,float,float)
73 DECL_GSISX(c,float,std::complex<float>)
74 DECL_GSISX(d,double,double)
75 DECL_GSISX(z,double,std::complex<double>)
76 
77 #endif
78 
79 template<typename MatrixType>
80 struct SluMatrixMapHelper;
81 
89 struct SluMatrix : SuperMatrix
90 {
91  SluMatrix()
92  {
93  Store = &storage;
94  }
95 
96  SluMatrix(const SluMatrix& other)
97  : SuperMatrix(other)
98  {
99  Store = &storage;
100  storage = other.storage;
101  }
102 
103  SluMatrix& operator=(const SluMatrix& other)
104  {
105  SuperMatrix::operator=(static_cast<const SuperMatrix&>(other));
106  Store = &storage;
107  storage = other.storage;
108  return *this;
109  }
110 
111  struct
112  {
113  union {int nnz;int lda;};
114  void *values;
115  int *innerInd;
116  int *outerInd;
117  } storage;
118 
119  void setStorageType(Stype_t t)
120  {
121  Stype = t;
122  if (t==SLU_NC || t==SLU_NR || t==SLU_DN)
123  Store = &storage;
124  else
125  {
126  eigen_assert(false && "storage type not supported");
127  Store = 0;
128  }
129  }
130 
131  template<typename Scalar>
132  void setScalarType()
133  {
134  if (internal::is_same<Scalar,float>::value)
135  Dtype = SLU_S;
136  else if (internal::is_same<Scalar,double>::value)
137  Dtype = SLU_D;
138  else if (internal::is_same<Scalar,std::complex<float> >::value)
139  Dtype = SLU_C;
140  else if (internal::is_same<Scalar,std::complex<double> >::value)
141  Dtype = SLU_Z;
142  else
143  {
144  eigen_assert(false && "Scalar type not supported by SuperLU");
145  }
146  }
147 
148  template<typename MatrixType>
149  static SluMatrix Map(MatrixBase<MatrixType>& _mat)
150  {
151  MatrixType& mat(_mat.derived());
152  eigen_assert( ((MatrixType::Flags&RowMajorBit)!=RowMajorBit) && "row-major dense matrices are not supported by SuperLU");
153  SluMatrix res;
154  res.setStorageType(SLU_DN);
155  res.setScalarType<typename MatrixType::Scalar>();
156  res.Mtype = SLU_GE;
157 
158  res.nrow = internal::convert_index<int>(mat.rows());
159  res.ncol = internal::convert_index<int>(mat.cols());
160 
161  res.storage.lda = internal::convert_index<int>(MatrixType::IsVectorAtCompileTime ? mat.size() : mat.outerStride());
162  res.storage.values = (void*)(mat.data());
163  return res;
164  }
165 
166  template<typename MatrixType>
167  static SluMatrix Map(SparseMatrixBase<MatrixType>& a_mat)
168  {
169  MatrixType &mat(a_mat.derived());
170  SluMatrix res;
171  if ((MatrixType::Flags&RowMajorBit)==RowMajorBit)
172  {
173  res.setStorageType(SLU_NR);
174  res.nrow = internal::convert_index<int>(mat.cols());
175  res.ncol = internal::convert_index<int>(mat.rows());
176  }
177  else
178  {
179  res.setStorageType(SLU_NC);
180  res.nrow = internal::convert_index<int>(mat.rows());
181  res.ncol = internal::convert_index<int>(mat.cols());
182  }
183 
184  res.Mtype = SLU_GE;
185 
186  res.storage.nnz = internal::convert_index<int>(mat.nonZeros());
187  res.storage.values = mat.valuePtr();
188  res.storage.innerInd = mat.innerIndexPtr();
189  res.storage.outerInd = mat.outerIndexPtr();
190 
191  res.setScalarType<typename MatrixType::Scalar>();
192 
193  // FIXME the following is not very accurate
194  if (MatrixType::Flags & Upper)
195  res.Mtype = SLU_TRU;
196  if (MatrixType::Flags & Lower)
197  res.Mtype = SLU_TRL;
198 
199  eigen_assert(((MatrixType::Flags & SelfAdjoint)==0) && "SelfAdjoint matrix shape not supported by SuperLU");
200 
201  return res;
202  }
203 };
204 
205 template<typename Scalar, int Rows, int Cols, int Options, int MRows, int MCols>
206 struct SluMatrixMapHelper<Matrix<Scalar,Rows,Cols,Options,MRows,MCols> >
207 {
208  typedef Matrix<Scalar,Rows,Cols,Options,MRows,MCols> MatrixType;
209  static void run(MatrixType& mat, SluMatrix& res)
210  {
211  eigen_assert( ((Options&RowMajor)!=RowMajor) && "row-major dense matrices is not supported by SuperLU");
212  res.setStorageType(SLU_DN);
213  res.setScalarType<Scalar>();
214  res.Mtype = SLU_GE;
215 
216  res.nrow = mat.rows();
217  res.ncol = mat.cols();
218 
219  res.storage.lda = mat.outerStride();
220  res.storage.values = mat.data();
221  }
222 };
223 
224 template<typename Derived>
225 struct SluMatrixMapHelper<SparseMatrixBase<Derived> >
226 {
227  typedef Derived MatrixType;
228  static void run(MatrixType& mat, SluMatrix& res)
229  {
230  if ((MatrixType::Flags&RowMajorBit)==RowMajorBit)
231  {
232  res.setStorageType(SLU_NR);
233  res.nrow = mat.cols();
234  res.ncol = mat.rows();
235  }
236  else
237  {
238  res.setStorageType(SLU_NC);
239  res.nrow = mat.rows();
240  res.ncol = mat.cols();
241  }
242 
243  res.Mtype = SLU_GE;
244 
245  res.storage.nnz = mat.nonZeros();
246  res.storage.values = mat.valuePtr();
247  res.storage.innerInd = mat.innerIndexPtr();
248  res.storage.outerInd = mat.outerIndexPtr();
249 
250  res.setScalarType<typename MatrixType::Scalar>();
251 
252  // FIXME the following is not very accurate
253  if (MatrixType::Flags & Upper)
254  res.Mtype = SLU_TRU;
255  if (MatrixType::Flags & Lower)
256  res.Mtype = SLU_TRL;
257 
258  eigen_assert(((MatrixType::Flags & SelfAdjoint)==0) && "SelfAdjoint matrix shape not supported by SuperLU");
259  }
260 };
261 
262 namespace internal {
263 
264 template<typename MatrixType>
265 SluMatrix asSluMatrix(MatrixType& mat)
266 {
267  return SluMatrix::Map(mat);
268 }
269 
271 template<typename Scalar, int Flags, typename Index>
272 MappedSparseMatrix<Scalar,Flags,Index> map_superlu(SluMatrix& sluMat)
273 {
274  eigen_assert((Flags&RowMajor)==RowMajor && sluMat.Stype == SLU_NR
275  || (Flags&ColMajor)==ColMajor && sluMat.Stype == SLU_NC);
276 
277  Index outerSize = (Flags&RowMajor)==RowMajor ? sluMat.ncol : sluMat.nrow;
278 
279  return MappedSparseMatrix<Scalar,Flags,Index>(
280  sluMat.nrow, sluMat.ncol, sluMat.storage.outerInd[outerSize],
281  sluMat.storage.outerInd, sluMat.storage.innerInd, reinterpret_cast<Scalar*>(sluMat.storage.values) );
282 }
283 
284 } // end namespace internal
285 
290 template<typename _MatrixType, typename Derived>
291 class SuperLUBase : public SparseSolverBase<Derived>
292 {
293  protected:
295  using Base::derived;
296  using Base::m_isInitialized;
297  public:
298  typedef _MatrixType MatrixType;
299  typedef typename MatrixType::Scalar Scalar;
300  typedef typename MatrixType::RealScalar RealScalar;
301  typedef typename MatrixType::StorageIndex StorageIndex;
307  enum {
308  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
309  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
310  };
311 
312  public:
313 
314  SuperLUBase() {}
315 
316  ~SuperLUBase()
317  {
318  clearFactors();
319  }
320 
321  inline Index rows() const { return m_matrix.rows(); }
322  inline Index cols() const { return m_matrix.cols(); }
323 
325  inline superlu_options_t& options() { return m_sluOptions; }
326 
333  {
334  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
335  return m_info;
336  }
337 
339  void compute(const MatrixType& matrix)
340  {
341  derived().analyzePattern(matrix);
342  derived().factorize(matrix);
343  }
344 
351  void analyzePattern(const MatrixType& /*matrix*/)
352  {
353  m_isInitialized = true;
354  m_info = Success;
355  m_analysisIsOk = true;
356  m_factorizationIsOk = false;
357  }
358 
359  template<typename Stream>
360  void dumpMemory(Stream& /*s*/)
361  {}
362 
363  protected:
364 
365  void initFactorization(const MatrixType& a)
366  {
367  set_default_options(&this->m_sluOptions);
368 
369  const Index size = a.rows();
370  m_matrix = a;
371 
372  m_sluA = internal::asSluMatrix(m_matrix);
373  clearFactors();
374 
375  m_p.resize(size);
376  m_q.resize(size);
377  m_sluRscale.resize(size);
378  m_sluCscale.resize(size);
379  m_sluEtree.resize(size);
380 
381  // set empty B and X
382  m_sluB.setStorageType(SLU_DN);
383  m_sluB.setScalarType<Scalar>();
384  m_sluB.Mtype = SLU_GE;
385  m_sluB.storage.values = 0;
386  m_sluB.nrow = 0;
387  m_sluB.ncol = 0;
388  m_sluB.storage.lda = internal::convert_index<int>(size);
389  m_sluX = m_sluB;
390 
391  m_extractedDataAreDirty = true;
392  }
393 
394  void init()
395  {
396  m_info = InvalidInput;
397  m_isInitialized = false;
398  m_sluL.Store = 0;
399  m_sluU.Store = 0;
400  }
401 
402  void extractData() const;
403 
404  void clearFactors()
405  {
406  if(m_sluL.Store)
407  Destroy_SuperNode_Matrix(&m_sluL);
408  if(m_sluU.Store)
409  Destroy_CompCol_Matrix(&m_sluU);
410 
411  m_sluL.Store = 0;
412  m_sluU.Store = 0;
413 
414  memset(&m_sluL,0,sizeof m_sluL);
415  memset(&m_sluU,0,sizeof m_sluU);
416  }
417 
418  // cached data to reduce reallocation, etc.
419  mutable LUMatrixType m_l;
420  mutable LUMatrixType m_u;
421  mutable IntColVectorType m_p;
422  mutable IntRowVectorType m_q;
423 
424  mutable LUMatrixType m_matrix; // copy of the factorized matrix
425  mutable SluMatrix m_sluA;
426  mutable SuperMatrix m_sluL, m_sluU;
427  mutable SluMatrix m_sluB, m_sluX;
428  mutable SuperLUStat_t m_sluStat;
429  mutable superlu_options_t m_sluOptions;
430  mutable std::vector<int> m_sluEtree;
431  mutable Matrix<RealScalar,Dynamic,1> m_sluRscale, m_sluCscale;
432  mutable Matrix<RealScalar,Dynamic,1> m_sluFerr, m_sluBerr;
433  mutable char m_sluEqued;
434 
435  mutable ComputationInfo m_info;
436  int m_factorizationIsOk;
437  int m_analysisIsOk;
438  mutable bool m_extractedDataAreDirty;
439 
440  private:
441  SuperLUBase(SuperLUBase& ) { }
442 };
443 
444 
461 template<typename _MatrixType>
462 class SuperLU : public SuperLUBase<_MatrixType,SuperLU<_MatrixType> >
463 {
464  public:
466  typedef _MatrixType MatrixType;
467  typedef typename Base::Scalar Scalar;
468  typedef typename Base::RealScalar RealScalar;
469  typedef typename Base::StorageIndex StorageIndex;
470  typedef typename Base::IntRowVectorType IntRowVectorType;
471  typedef typename Base::IntColVectorType IntColVectorType;
472  typedef typename Base::PermutationMap PermutationMap;
473  typedef typename Base::LUMatrixType LUMatrixType;
476 
477  public:
478  using Base::_solve_impl;
479 
480  SuperLU() : Base() { init(); }
481 
482  explicit SuperLU(const MatrixType& matrix) : Base()
483  {
484  init();
485  Base::compute(matrix);
486  }
487 
488  ~SuperLU()
489  {
490  }
491 
498  void analyzePattern(const MatrixType& matrix)
499  {
500  m_info = InvalidInput;
501  m_isInitialized = false;
502  Base::analyzePattern(matrix);
503  }
504 
511  void factorize(const MatrixType& matrix);
512 
514  template<typename Rhs,typename Dest>
515  void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const;
516 
517  inline const LMatrixType& matrixL() const
518  {
519  if (m_extractedDataAreDirty) this->extractData();
520  return m_l;
521  }
522 
523  inline const UMatrixType& matrixU() const
524  {
525  if (m_extractedDataAreDirty) this->extractData();
526  return m_u;
527  }
528 
529  inline const IntColVectorType& permutationP() const
530  {
531  if (m_extractedDataAreDirty) this->extractData();
532  return m_p;
533  }
534 
535  inline const IntRowVectorType& permutationQ() const
536  {
537  if (m_extractedDataAreDirty) this->extractData();
538  return m_q;
539  }
540 
541  Scalar determinant() const;
542 
543  protected:
544 
545  using Base::m_matrix;
546  using Base::m_sluOptions;
547  using Base::m_sluA;
548  using Base::m_sluB;
549  using Base::m_sluX;
550  using Base::m_p;
551  using Base::m_q;
552  using Base::m_sluEtree;
553  using Base::m_sluEqued;
554  using Base::m_sluRscale;
555  using Base::m_sluCscale;
556  using Base::m_sluL;
557  using Base::m_sluU;
558  using Base::m_sluStat;
559  using Base::m_sluFerr;
560  using Base::m_sluBerr;
561  using Base::m_l;
562  using Base::m_u;
563 
564  using Base::m_analysisIsOk;
565  using Base::m_factorizationIsOk;
566  using Base::m_extractedDataAreDirty;
567  using Base::m_isInitialized;
568  using Base::m_info;
569 
570  void init()
571  {
572  Base::init();
573 
574  set_default_options(&this->m_sluOptions);
575  m_sluOptions.PrintStat = NO;
576  m_sluOptions.ConditionNumber = NO;
577  m_sluOptions.Trans = NOTRANS;
578  m_sluOptions.ColPerm = COLAMD;
579  }
580 
581 
582  private:
583  SuperLU(SuperLU& ) { }
584 };
585 
586 template<typename MatrixType>
587 void SuperLU<MatrixType>::factorize(const MatrixType& a)
588 {
589  eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
590  if(!m_analysisIsOk)
591  {
592  m_info = InvalidInput;
593  return;
594  }
595 
596  this->initFactorization(a);
597 
598  m_sluOptions.ColPerm = COLAMD;
599  int info = 0;
600  RealScalar recip_pivot_growth, rcond;
601  RealScalar ferr, berr;
602 
603  StatInit(&m_sluStat);
604  SuperLU_gssvx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0],
605  &m_sluEqued, &m_sluRscale[0], &m_sluCscale[0],
606  &m_sluL, &m_sluU,
607  NULL, 0,
608  &m_sluB, &m_sluX,
609  &recip_pivot_growth, &rcond,
610  &ferr, &berr,
611  &m_sluStat, &info, Scalar());
612  StatFree(&m_sluStat);
613 
614  m_extractedDataAreDirty = true;
615 
616  // FIXME how to better check for errors ???
617  m_info = info == 0 ? Success : NumericalIssue;
618  m_factorizationIsOk = true;
619 }
620 
621 template<typename MatrixType>
622 template<typename Rhs,typename Dest>
624 {
625  eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or analyzePattern()/factorize()");
626 
627  const Index size = m_matrix.rows();
628  const Index rhsCols = b.cols();
629  eigen_assert(size==b.rows());
630 
631  m_sluOptions.Trans = NOTRANS;
632  m_sluOptions.Fact = FACTORED;
633  m_sluOptions.IterRefine = NOREFINE;
634 
635 
636  m_sluFerr.resize(rhsCols);
637  m_sluBerr.resize(rhsCols);
638 
641 
642  m_sluB = SluMatrix::Map(b_ref.const_cast_derived());
643  m_sluX = SluMatrix::Map(x_ref.const_cast_derived());
644 
645  typename Rhs::PlainObject b_cpy;
646  if(m_sluEqued!='N')
647  {
648  b_cpy = b;
649  m_sluB = SluMatrix::Map(b_cpy.const_cast_derived());
650  }
651 
652  StatInit(&m_sluStat);
653  int info = 0;
654  RealScalar recip_pivot_growth, rcond;
655  SuperLU_gssvx(&m_sluOptions, &m_sluA,
656  m_q.data(), m_p.data(),
657  &m_sluEtree[0], &m_sluEqued,
658  &m_sluRscale[0], &m_sluCscale[0],
659  &m_sluL, &m_sluU,
660  NULL, 0,
661  &m_sluB, &m_sluX,
662  &recip_pivot_growth, &rcond,
663  &m_sluFerr[0], &m_sluBerr[0],
664  &m_sluStat, &info, Scalar());
665  StatFree(&m_sluStat);
666 
667  if(x.derived().data() != x_ref.data())
668  x = x_ref;
669 
670  m_info = info==0 ? Success : NumericalIssue;
671 }
672 
673 // the code of this extractData() function has been adapted from the SuperLU's Matlab support code,
674 //
675 // Copyright (c) 1994 by Xerox Corporation. All rights reserved.
676 //
677 // THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
678 // EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK.
679 //
680 template<typename MatrixType, typename Derived>
681 void SuperLUBase<MatrixType,Derived>::extractData() const
682 {
683  eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for extracting factors, you must first call either compute() or analyzePattern()/factorize()");
684  if (m_extractedDataAreDirty)
685  {
686  int upper;
687  int fsupc, istart, nsupr;
688  int lastl = 0, lastu = 0;
689  SCformat *Lstore = static_cast<SCformat*>(m_sluL.Store);
690  NCformat *Ustore = static_cast<NCformat*>(m_sluU.Store);
691  Scalar *SNptr;
692 
693  const Index size = m_matrix.rows();
694  m_l.resize(size,size);
695  m_l.resizeNonZeros(Lstore->nnz);
696  m_u.resize(size,size);
697  m_u.resizeNonZeros(Ustore->nnz);
698 
699  int* Lcol = m_l.outerIndexPtr();
700  int* Lrow = m_l.innerIndexPtr();
701  Scalar* Lval = m_l.valuePtr();
702 
703  int* Ucol = m_u.outerIndexPtr();
704  int* Urow = m_u.innerIndexPtr();
705  Scalar* Uval = m_u.valuePtr();
706 
707  Ucol[0] = 0;
708  Ucol[0] = 0;
709 
710  /* for each supernode */
711  for (int k = 0; k <= Lstore->nsuper; ++k)
712  {
713  fsupc = L_FST_SUPC(k);
714  istart = L_SUB_START(fsupc);
715  nsupr = L_SUB_START(fsupc+1) - istart;
716  upper = 1;
717 
718  /* for each column in the supernode */
719  for (int j = fsupc; j < L_FST_SUPC(k+1); ++j)
720  {
721  SNptr = &((Scalar*)Lstore->nzval)[L_NZ_START(j)];
722 
723  /* Extract U */
724  for (int i = U_NZ_START(j); i < U_NZ_START(j+1); ++i)
725  {
726  Uval[lastu] = ((Scalar*)Ustore->nzval)[i];
727  /* Matlab doesn't like explicit zero. */
728  if (Uval[lastu] != 0.0)
729  Urow[lastu++] = U_SUB(i);
730  }
731  for (int i = 0; i < upper; ++i)
732  {
733  /* upper triangle in the supernode */
734  Uval[lastu] = SNptr[i];
735  /* Matlab doesn't like explicit zero. */
736  if (Uval[lastu] != 0.0)
737  Urow[lastu++] = L_SUB(istart+i);
738  }
739  Ucol[j+1] = lastu;
740 
741  /* Extract L */
742  Lval[lastl] = 1.0; /* unit diagonal */
743  Lrow[lastl++] = L_SUB(istart + upper - 1);
744  for (int i = upper; i < nsupr; ++i)
745  {
746  Lval[lastl] = SNptr[i];
747  /* Matlab doesn't like explicit zero. */
748  if (Lval[lastl] != 0.0)
749  Lrow[lastl++] = L_SUB(istart+i);
750  }
751  Lcol[j+1] = lastl;
752 
753  ++upper;
754  } /* for j ... */
755 
756  } /* for k ... */
757 
758  // squeeze the matrices :
759  m_l.resizeNonZeros(lastl);
760  m_u.resizeNonZeros(lastu);
761 
762  m_extractedDataAreDirty = false;
763  }
764 }
765 
766 template<typename MatrixType>
767 typename SuperLU<MatrixType>::Scalar SuperLU<MatrixType>::determinant() const
768 {
769  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()");
770 
771  if (m_extractedDataAreDirty)
772  this->extractData();
773 
774  Scalar det = Scalar(1);
775  for (int j=0; j<m_u.cols(); ++j)
776  {
777  if (m_u.outerIndexPtr()[j+1]-m_u.outerIndexPtr()[j] > 0)
778  {
779  int lastId = m_u.outerIndexPtr()[j+1]-1;
780  eigen_assert(m_u.innerIndexPtr()[lastId]<=j);
781  if (m_u.innerIndexPtr()[lastId]==j)
782  det *= m_u.valuePtr()[lastId];
783  }
784  }
785  if(PermutationMap(m_p.data(),m_p.size()).determinant()*PermutationMap(m_q.data(),m_q.size()).determinant()<0)
786  det = -det;
787  if(m_sluEqued!='N')
788  return det/m_sluRscale.prod()/m_sluCscale.prod();
789  else
790  return det;
791 }
792 
793 #ifdef EIGEN_PARSED_BY_DOXYGEN
794 #define EIGEN_SUPERLU_HAS_ILU
795 #endif
796 
797 #ifdef EIGEN_SUPERLU_HAS_ILU
798 
815 template<typename _MatrixType>
816 class SuperILU : public SuperLUBase<_MatrixType,SuperILU<_MatrixType> >
817 {
818  public:
820  typedef _MatrixType MatrixType;
821  typedef typename Base::Scalar Scalar;
822  typedef typename Base::RealScalar RealScalar;
823 
824  public:
825  using Base::_solve_impl;
826 
827  SuperILU() : Base() { init(); }
828 
829  SuperILU(const MatrixType& matrix) : Base()
830  {
831  init();
832  Base::compute(matrix);
833  }
834 
835  ~SuperILU()
836  {
837  }
838 
845  void analyzePattern(const MatrixType& matrix)
846  {
847  Base::analyzePattern(matrix);
848  }
849 
856  void factorize(const MatrixType& matrix);
857 
858  #ifndef EIGEN_PARSED_BY_DOXYGEN
859 
860  template<typename Rhs,typename Dest>
861  void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const;
862  #endif // EIGEN_PARSED_BY_DOXYGEN
863 
864  protected:
865 
866  using Base::m_matrix;
867  using Base::m_sluOptions;
868  using Base::m_sluA;
869  using Base::m_sluB;
870  using Base::m_sluX;
871  using Base::m_p;
872  using Base::m_q;
873  using Base::m_sluEtree;
874  using Base::m_sluEqued;
875  using Base::m_sluRscale;
876  using Base::m_sluCscale;
877  using Base::m_sluL;
878  using Base::m_sluU;
879  using Base::m_sluStat;
880  using Base::m_sluFerr;
881  using Base::m_sluBerr;
882  using Base::m_l;
883  using Base::m_u;
884 
885  using Base::m_analysisIsOk;
886  using Base::m_factorizationIsOk;
887  using Base::m_extractedDataAreDirty;
888  using Base::m_isInitialized;
889  using Base::m_info;
890 
891  void init()
892  {
893  Base::init();
894 
895  ilu_set_default_options(&m_sluOptions);
896  m_sluOptions.PrintStat = NO;
897  m_sluOptions.ConditionNumber = NO;
898  m_sluOptions.Trans = NOTRANS;
899  m_sluOptions.ColPerm = MMD_AT_PLUS_A;
900 
901  // no attempt to preserve column sum
902  m_sluOptions.ILU_MILU = SILU;
903  // only basic ILU(k) support -- no direct control over memory consumption
904  // better to use ILU_DropRule = DROP_BASIC | DROP_AREA
905  // and set ILU_FillFactor to max memory growth
906  m_sluOptions.ILU_DropRule = DROP_BASIC;
907  m_sluOptions.ILU_DropTol = NumTraits<Scalar>::dummy_precision()*10;
908  }
909 
910  private:
911  SuperILU(SuperILU& ) { }
912 };
913 
914 template<typename MatrixType>
915 void SuperILU<MatrixType>::factorize(const MatrixType& a)
916 {
917  eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
918  if(!m_analysisIsOk)
919  {
920  m_info = InvalidInput;
921  return;
922  }
923 
924  this->initFactorization(a);
925 
926  int info = 0;
927  RealScalar recip_pivot_growth, rcond;
928 
929  StatInit(&m_sluStat);
930  SuperLU_gsisx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0],
931  &m_sluEqued, &m_sluRscale[0], &m_sluCscale[0],
932  &m_sluL, &m_sluU,
933  NULL, 0,
934  &m_sluB, &m_sluX,
935  &recip_pivot_growth, &rcond,
936  &m_sluStat, &info, Scalar());
937  StatFree(&m_sluStat);
938 
939  // FIXME how to better check for errors ???
940  m_info = info == 0 ? Success : NumericalIssue;
941  m_factorizationIsOk = true;
942 }
943 
944 template<typename MatrixType>
945 template<typename Rhs,typename Dest>
947 {
948  eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or analyzePattern()/factorize()");
949 
950  const int size = m_matrix.rows();
951  const int rhsCols = b.cols();
952  eigen_assert(size==b.rows());
953 
954  m_sluOptions.Trans = NOTRANS;
955  m_sluOptions.Fact = FACTORED;
956  m_sluOptions.IterRefine = NOREFINE;
957 
958  m_sluFerr.resize(rhsCols);
959  m_sluBerr.resize(rhsCols);
960 
963 
964  m_sluB = SluMatrix::Map(b_ref.const_cast_derived());
965  m_sluX = SluMatrix::Map(x_ref.const_cast_derived());
966 
967  typename Rhs::PlainObject b_cpy;
968  if(m_sluEqued!='N')
969  {
970  b_cpy = b;
971  m_sluB = SluMatrix::Map(b_cpy.const_cast_derived());
972  }
973 
974  int info = 0;
975  RealScalar recip_pivot_growth, rcond;
976 
977  StatInit(&m_sluStat);
978  SuperLU_gsisx(&m_sluOptions, &m_sluA,
979  m_q.data(), m_p.data(),
980  &m_sluEtree[0], &m_sluEqued,
981  &m_sluRscale[0], &m_sluCscale[0],
982  &m_sluL, &m_sluU,
983  NULL, 0,
984  &m_sluB, &m_sluX,
985  &recip_pivot_growth, &rcond,
986  &m_sluStat, &info, Scalar());
987  StatFree(&m_sluStat);
988 
989  if(&x.coeffRef(0) != x_ref.data())
990  x = x_ref;
991 
992  m_info = info==0 ? Success : NumericalIssue;
993 }
994 #endif
995 
996 } // end namespace Eigen
997 
998 #endif // EIGEN_SUPERLUSUPPORT_H
Index cols() const
Definition: SparseMatrix.h:133
void compute(const MatrixType &matrix)
Definition: SuperLUSupport.h:339
void analyzePattern(const MatrixType &matrix)
Definition: SuperLUSupport.h:845
A sparse direct LU factorization and solver based on the SuperLU library.
Definition: SuperLUSupport.h:462
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:89
A sparse direct incomplete LU factorization and solver based on the SuperLU library.
Definition: SuperLUSupport.h:816
A base class for sparse solvers.
Definition: SparseSolverBase.h:53
Definition: LDLT.h:16
Definition: StdDeque.h:58
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:107
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: SuperLUSupport.h:332
const unsigned int RowMajorBit
Definition: Constants.h:61
void resize(Index rows, Index cols)
Definition: PlainObjectBase.h:252
void analyzePattern(const MatrixType &)
Definition: SuperLUSupport.h:351
Definition: Constants.h:320
Definition: Constants.h:204
Definition: Constants.h:322
void factorize(const MatrixType &matrix)
Definition: SuperLUSupport.h:915
Definition: Constants.h:434
The base class for the direct and incomplete LU factorization of SuperLU.
Definition: SuperLUSupport.h:291
Definition: Constants.h:439
Definition: Constants.h:432
A matrix or vector expression mapping an existing expression.
Definition: Ref.h:186
Definition: Eigen_Colamd.h:54
Expression of a triangular part in a matrix.
Definition: TriangularMatrix.h:186
void factorize(const MatrixType &matrix)
Definition: SuperLUSupport.h:587
Definition: Constants.h:220
superlu_options_t & options()
Definition: SuperLUSupport.h:325
Definition: Constants.h:206
ComputationInfo
Definition: Constants.h:430
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
void analyzePattern(const MatrixType &matrix)
Definition: SuperLUSupport.h:498
Index rows() const
Definition: SparseMatrix.h:131