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