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-2011 Gael Guennebaud <gael.guennebaud@inria.fr>
5 //
6 // Eigen is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 3 of the License, or (at your option) any later version.
10 //
11 // Alternatively, you can redistribute it and/or
12 // modify it under the terms of the GNU General Public License as
13 // published by the Free Software Foundation; either version 2 of
14 // the License, or (at your option) any later version.
15 //
16 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
17 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
19 // GNU General Public License for more details.
20 //
21 // You should have received a copy of the GNU Lesser General Public
22 // License and a copy of the GNU General Public License along with
23 // Eigen. If not, see <http://www.gnu.org/licenses/>.
24 
25 #ifndef EIGEN_SUPERLUSUPPORT_H
26 #define EIGEN_SUPERLUSUPPORT_H
27 
28 namespace Eigen {
29 
30 #define DECL_GSSVX(PREFIX,FLOATTYPE,KEYTYPE) \
31  extern "C" { \
32  typedef struct { FLOATTYPE for_lu; FLOATTYPE total_needed; int expansions; } PREFIX##mem_usage_t; \
33  extern void PREFIX##gssvx(superlu_options_t *, SuperMatrix *, int *, int *, int *, \
34  char *, FLOATTYPE *, FLOATTYPE *, SuperMatrix *, SuperMatrix *, \
35  void *, int, SuperMatrix *, SuperMatrix *, \
36  FLOATTYPE *, FLOATTYPE *, FLOATTYPE *, FLOATTYPE *, \
37  PREFIX##mem_usage_t *, SuperLUStat_t *, int *); \
38  } \
39  inline float SuperLU_gssvx(superlu_options_t *options, SuperMatrix *A, \
40  int *perm_c, int *perm_r, int *etree, char *equed, \
41  FLOATTYPE *R, FLOATTYPE *C, SuperMatrix *L, \
42  SuperMatrix *U, void *work, int lwork, \
43  SuperMatrix *B, SuperMatrix *X, \
44  FLOATTYPE *recip_pivot_growth, \
45  FLOATTYPE *rcond, FLOATTYPE *ferr, FLOATTYPE *berr, \
46  SuperLUStat_t *stats, int *info, KEYTYPE) { \
47  PREFIX##mem_usage_t mem_usage; \
48  PREFIX##gssvx(options, A, perm_c, perm_r, etree, equed, R, C, L, \
49  U, work, lwork, B, X, recip_pivot_growth, rcond, \
50  ferr, berr, &mem_usage, stats, info); \
51  return mem_usage.for_lu; /* bytes used by the factor storage */ \
52  }
53 
54 DECL_GSSVX(s,float,float)
55 DECL_GSSVX(c,float,std::complex<float>)
56 DECL_GSSVX(d,double,double)
57 DECL_GSSVX(z,double,std::complex<double>)
58 
59 #ifdef MILU_ALPHA
60 #define EIGEN_SUPERLU_HAS_ILU
61 #endif
62 
63 #ifdef EIGEN_SUPERLU_HAS_ILU
64 
65 // similarly for the incomplete factorization using gsisx
66 #define DECL_GSISX(PREFIX,FLOATTYPE,KEYTYPE) \
67  extern "C" { \
68  extern void PREFIX##gsisx(superlu_options_t *, SuperMatrix *, int *, int *, int *, \
69  char *, FLOATTYPE *, FLOATTYPE *, SuperMatrix *, SuperMatrix *, \
70  void *, int, SuperMatrix *, SuperMatrix *, FLOATTYPE *, FLOATTYPE *, \
71  PREFIX##mem_usage_t *, SuperLUStat_t *, int *); \
72  } \
73  inline float SuperLU_gsisx(superlu_options_t *options, SuperMatrix *A, \
74  int *perm_c, int *perm_r, int *etree, char *equed, \
75  FLOATTYPE *R, FLOATTYPE *C, SuperMatrix *L, \
76  SuperMatrix *U, void *work, int lwork, \
77  SuperMatrix *B, SuperMatrix *X, \
78  FLOATTYPE *recip_pivot_growth, \
79  FLOATTYPE *rcond, \
80  SuperLUStat_t *stats, int *info, KEYTYPE) { \
81  PREFIX##mem_usage_t mem_usage; \
82  PREFIX##gsisx(options, A, perm_c, perm_r, etree, equed, R, C, L, \
83  U, work, lwork, B, X, recip_pivot_growth, rcond, \
84  &mem_usage, stats, info); \
85  return mem_usage.for_lu; /* bytes used by the factor storage */ \
86  }
87 
88 DECL_GSISX(s,float,float)
89 DECL_GSISX(c,float,std::complex<float>)
90 DECL_GSISX(d,double,double)
91 DECL_GSISX(z,double,std::complex<double>)
92 
93 #endif
94 
95 template<typename MatrixType>
96 struct SluMatrixMapHelper;
97 
106 {
108  {
109  Store = &storage;
110  }
111 
112  SluMatrix(const SluMatrix& other)
113  : SuperMatrix(other)
114  {
115  Store = &storage;
116  storage = other.storage;
117  }
118 
120  {
121  SuperMatrix::operator=(static_cast<const SuperMatrix&>(other));
122  Store = &storage;
123  storage = other.storage;
124  return *this;
125  }
126 
127  struct
128  {
129  union {int nnz;int lda;};
130  void *values;
131  int *innerInd;
132  int *outerInd;
133  } storage;
134 
135  void setStorageType(Stype_t t)
136  {
137  Stype = t;
138  if (t==SLU_NC || t==SLU_NR || t==SLU_DN)
139  Store = &storage;
140  else
141  {
142  eigen_assert(false && "storage type not supported");
143  Store = 0;
144  }
145  }
146 
147  template<typename Scalar>
149  {
150  if (internal::is_same<Scalar,float>::value)
151  Dtype = SLU_S;
152  else if (internal::is_same<Scalar,double>::value)
153  Dtype = SLU_D;
154  else if (internal::is_same<Scalar,std::complex<float> >::value)
155  Dtype = SLU_C;
156  else if (internal::is_same<Scalar,std::complex<double> >::value)
157  Dtype = SLU_Z;
158  else
159  {
160  eigen_assert(false && "Scalar type not supported by SuperLU");
161  }
162  }
163 
164  template<typename MatrixType>
166  {
167  MatrixType& mat(_mat.derived());
168  eigen_assert( ((MatrixType::Flags&RowMajorBit)!=RowMajorBit) && "row-major dense matrices are not supported by SuperLU");
169  SluMatrix res;
170  res.setStorageType(SLU_DN);
171  res.setScalarType<typename MatrixType::Scalar>();
172  res.Mtype = SLU_GE;
173 
174  res.nrow = mat.rows();
175  res.ncol = mat.cols();
176 
177  res.storage.lda = MatrixType::IsVectorAtCompileTime ? mat.size() : mat.outerStride();
178  res.storage.values = mat.data();
179  return res;
180  }
181 
182  template<typename MatrixType>
184  {
185  SluMatrix res;
186  if ((MatrixType::Flags&RowMajorBit)==RowMajorBit)
187  {
188  res.setStorageType(SLU_NR);
189  res.nrow = mat.cols();
190  res.ncol = mat.rows();
191  }
192  else
193  {
194  res.setStorageType(SLU_NC);
195  res.nrow = mat.rows();
196  res.ncol = mat.cols();
197  }
198 
199  res.Mtype = SLU_GE;
200 
201  res.storage.nnz = mat.nonZeros();
202  res.storage.values = mat.derived().valuePtr();
203  res.storage.innerInd = mat.derived().innerIndexPtr();
204  res.storage.outerInd = mat.derived().outerIndexPtr();
205 
206  res.setScalarType<typename MatrixType::Scalar>();
207 
208  // FIXME the following is not very accurate
209  if (MatrixType::Flags & Upper)
210  res.Mtype = SLU_TRU;
211  if (MatrixType::Flags & Lower)
212  res.Mtype = SLU_TRL;
213 
214  eigen_assert(((MatrixType::Flags & SelfAdjoint)==0) && "SelfAdjoint matrix shape not supported by SuperLU");
215 
216  return res;
217  }
218 };
219 
220 template<typename Scalar, int Rows, int Cols, int Options, int MRows, int MCols>
221 struct SluMatrixMapHelper<Matrix<Scalar,Rows,Cols,Options,MRows,MCols> >
222 {
224  static void run(MatrixType& mat, SluMatrix& res)
225  {
226  eigen_assert( ((Options&RowMajor)!=RowMajor) && "row-major dense matrices is not supported by SuperLU");
227  res.setStorageType(SLU_DN);
228  res.setScalarType<Scalar>();
229  res.Mtype = SLU_GE;
230 
231  res.nrow = mat.rows();
232  res.ncol = mat.cols();
233 
234  res.storage.lda = mat.outerStride();
235  res.storage.values = mat.data();
236  }
237 };
238 
239 template<typename Derived>
240 struct SluMatrixMapHelper<SparseMatrixBase<Derived> >
241 {
242  typedef Derived MatrixType;
243  static void run(MatrixType& mat, SluMatrix& res)
244  {
245  if ((MatrixType::Flags&RowMajorBit)==RowMajorBit)
246  {
247  res.setStorageType(SLU_NR);
248  res.nrow = mat.cols();
249  res.ncol = mat.rows();
250  }
251  else
252  {
253  res.setStorageType(SLU_NC);
254  res.nrow = mat.rows();
255  res.ncol = mat.cols();
256  }
257 
258  res.Mtype = SLU_GE;
259 
260  res.storage.nnz = mat.nonZeros();
261  res.storage.values = mat.valuePtr();
262  res.storage.innerInd = mat.innerIndexPtr();
263  res.storage.outerInd = mat.outerIndexPtr();
264 
265  res.setScalarType<typename MatrixType::Scalar>();
266 
267  // FIXME the following is not very accurate
268  if (MatrixType::Flags & Upper)
269  res.Mtype = SLU_TRU;
270  if (MatrixType::Flags & Lower)
271  res.Mtype = SLU_TRL;
272 
273  eigen_assert(((MatrixType::Flags & SelfAdjoint)==0) && "SelfAdjoint matrix shape not supported by SuperLU");
274  }
275 };
276 
277 namespace internal {
278 
279 template<typename MatrixType>
280 SluMatrix asSluMatrix(MatrixType& mat)
281 {
282  return SluMatrix::Map(mat);
283 }
284 
286 template<typename Scalar, int Flags, typename Index>
288 {
289  eigen_assert((Flags&RowMajor)==RowMajor && sluMat.Stype == SLU_NR
290  || (Flags&ColMajor)==ColMajor && sluMat.Stype == SLU_NC);
291 
292  Index outerSize = (Flags&RowMajor)==RowMajor ? sluMat.ncol : sluMat.nrow;
293 
295  sluMat.nrow, sluMat.ncol, sluMat.storage.outerInd[outerSize],
296  sluMat.storage.outerInd, sluMat.storage.innerInd, reinterpret_cast<Scalar*>(sluMat.storage.values) );
297 }
298 
299 } // end namespace internal
300 
305 template<typename _MatrixType, typename Derived>
306 class SuperLUBase : internal::noncopyable
307 {
308  public:
309  typedef _MatrixType MatrixType;
310  typedef typename MatrixType::Scalar Scalar;
311  typedef typename MatrixType::RealScalar RealScalar;
312  typedef typename MatrixType::Index Index;
317 
318  public:
319 
321 
323  {
324  clearFactors();
325  }
326 
327  Derived& derived() { return *static_cast<Derived*>(this); }
328  const Derived& derived() const { return *static_cast<const Derived*>(this); }
329 
330  inline Index rows() const { return m_matrix.rows(); }
331  inline Index cols() const { return m_matrix.cols(); }
332 
334  inline superlu_options_t& options() { return m_sluOptions; }
335 
342  {
343  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
344  return m_info;
345  }
346 
348  void compute(const MatrixType& matrix)
349  {
350  derived().analyzePattern(matrix);
351  derived().factorize(matrix);
352  }
353 
358  template<typename Rhs>
359  inline const internal::solve_retval<SuperLUBase, Rhs> solve(const MatrixBase<Rhs>& b) const
360  {
361  eigen_assert(m_isInitialized && "SuperLU is not initialized.");
362  eigen_assert(rows()==b.rows()
363  && "SuperLU::solve(): invalid number of rows of the right hand side matrix b");
364  return internal::solve_retval<SuperLUBase, Rhs>(*this, b.derived());
365  }
366 
371 // template<typename Rhs>
372 // inline const internal::sparse_solve_retval<SuperLU, Rhs> solve(const SparseMatrixBase<Rhs>& b) const
373 // {
374 // eigen_assert(m_isInitialized && "SuperLU is not initialized.");
375 // eigen_assert(rows()==b.rows()
376 // && "SuperLU::solve(): invalid number of rows of the right hand side matrix b");
377 // return internal::sparse_solve_retval<SuperLU, Rhs>(*this, b.derived());
378 // }
379 
386  void analyzePattern(const MatrixType& /*matrix*/)
387  {
388  m_isInitialized = true;
389  m_info = Success;
390  m_analysisIsOk = true;
391  m_factorizationIsOk = false;
392  }
393 
394  template<typename Stream>
395  void dumpMemory(Stream& s)
396  {}
397 
398  protected:
399 
401  {
402  set_default_options(&this->m_sluOptions);
403 
404  const int size = a.rows();
405  m_matrix = a;
406 
408  clearFactors();
409 
410  m_p.resize(size);
411  m_q.resize(size);
412  m_sluRscale.resize(size);
413  m_sluCscale.resize(size);
414  m_sluEtree.resize(size);
415 
416  // set empty B and X
417  m_sluB.setStorageType(SLU_DN);
419  m_sluB.Mtype = SLU_GE;
420  m_sluB.storage.values = 0;
421  m_sluB.nrow = 0;
422  m_sluB.ncol = 0;
423  m_sluB.storage.lda = size;
424  m_sluX = m_sluB;
425 
427  }
428 
429  void init()
430  {
432  m_isInitialized = false;
433  m_sluL.Store = 0;
434  m_sluU.Store = 0;
435  }
436 
437  void extractData() const;
438 
440  {
441  if(m_sluL.Store)
442  Destroy_SuperNode_Matrix(&m_sluL);
443  if(m_sluU.Store)
444  Destroy_CompCol_Matrix(&m_sluU);
445 
446  m_sluL.Store = 0;
447  m_sluU.Store = 0;
448 
449  memset(&m_sluL,0,sizeof m_sluL);
450  memset(&m_sluU,0,sizeof m_sluU);
451  }
452 
453  // cached data to reduce reallocation, etc.
454  mutable LUMatrixType m_l;
455  mutable LUMatrixType m_u;
458 
459  mutable LUMatrixType m_matrix; // copy of the factorized matrix
460  mutable SluMatrix m_sluA;
465  mutable std::vector<int> m_sluEtree;
468  mutable char m_sluEqued;
469 
475 
476  private:
477  SuperLUBase(SuperLUBase& ) { }
478 };
479 
480 
493 template<typename _MatrixType>
494 class SuperLU : public SuperLUBase<_MatrixType,SuperLU<_MatrixType> >
495 {
496  public:
498  typedef _MatrixType MatrixType;
499  typedef typename Base::Scalar Scalar;
500  typedef typename Base::RealScalar RealScalar;
501  typedef typename Base::Index Index;
507 
508  public:
509 
510  SuperLU() : Base() { init(); }
511 
512  SuperLU(const MatrixType& matrix) : Base()
513  {
514  Base::init();
515  compute(matrix);
516  }
517 
519  {
520  }
521 
528  void analyzePattern(const MatrixType& matrix)
529  {
531  m_isInitialized = false;
532  Base::analyzePattern(matrix);
533  }
534 
541  void factorize(const MatrixType& matrix);
542 
543  #ifndef EIGEN_PARSED_BY_DOXYGEN
544 
545  template<typename Rhs,typename Dest>
546  void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const;
547  #endif // EIGEN_PARSED_BY_DOXYGEN
548 
549  inline const LMatrixType& matrixL() const
550  {
552  return m_l;
553  }
554 
555  inline const UMatrixType& matrixU() const
556  {
558  return m_u;
559  }
560 
561  inline const IntColVectorType& permutationP() const
562  {
564  return m_p;
565  }
566 
567  inline const IntRowVectorType& permutationQ() const
568  {
570  return m_q;
571  }
572 
573  Scalar determinant() const;
574 
575  protected:
576 
577  using Base::m_matrix;
578  using Base::m_sluOptions;
579  using Base::m_sluA;
580  using Base::m_sluB;
581  using Base::m_sluX;
582  using Base::m_p;
583  using Base::m_q;
584  using Base::m_sluEtree;
585  using Base::m_sluEqued;
586  using Base::m_sluRscale;
587  using Base::m_sluCscale;
588  using Base::m_sluL;
589  using Base::m_sluU;
590  using Base::m_sluStat;
591  using Base::m_sluFerr;
592  using Base::m_sluBerr;
593  using Base::m_l;
594  using Base::m_u;
595 
596  using Base::m_analysisIsOk;
599  using Base::m_isInitialized;
600  using Base::m_info;
601 
602  void init()
603  {
604  Base::init();
605 
606  set_default_options(&this->m_sluOptions);
607  m_sluOptions.PrintStat = NO;
608  m_sluOptions.ConditionNumber = NO;
609  m_sluOptions.Trans = NOTRANS;
610  m_sluOptions.ColPerm = COLAMD;
611  }
612 
613 
614  private:
615  SuperLU(SuperLU& ) { }
616 };
617 
618 template<typename MatrixType>
620 {
621  eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
622  if(!m_analysisIsOk)
623  {
624  m_info = InvalidInput;
625  return;
626  }
627 
628  this->initFactorization(a);
629 
630  int info = 0;
631  RealScalar recip_pivot_growth, rcond;
632  RealScalar ferr, berr;
633 
634  StatInit(&m_sluStat);
635  SuperLU_gssvx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0],
636  &m_sluEqued, &m_sluRscale[0], &m_sluCscale[0],
637  &m_sluL, &m_sluU,
638  NULL, 0,
639  &m_sluB, &m_sluX,
640  &recip_pivot_growth, &rcond,
641  &ferr, &berr,
642  &m_sluStat, &info, Scalar());
643  StatFree(&m_sluStat);
644 
645  m_extractedDataAreDirty = true;
646 
647  // FIXME how to better check for errors ???
648  m_info = info == 0 ? Success : NumericalIssue;
649  m_factorizationIsOk = true;
650 }
651 
652 template<typename MatrixType>
653 template<typename Rhs,typename Dest>
655 {
656  eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or analyzePattern()/factorize()");
657 
658  const int size = m_matrix.rows();
659  const int rhsCols = b.cols();
660  eigen_assert(size==b.rows());
661 
662  m_sluOptions.Trans = NOTRANS;
663  m_sluOptions.Fact = FACTORED;
664  m_sluOptions.IterRefine = NOREFINE;
665 
666 
667  m_sluFerr.resize(rhsCols);
668  m_sluBerr.resize(rhsCols);
669  m_sluB = SluMatrix::Map(b.const_cast_derived());
670  m_sluX = SluMatrix::Map(x.derived());
671 
672  typename Rhs::PlainObject b_cpy;
673  if(m_sluEqued!='N')
674  {
675  b_cpy = b;
676  m_sluB = SluMatrix::Map(b_cpy.const_cast_derived());
677  }
678 
679  StatInit(&m_sluStat);
680  int info = 0;
681  RealScalar recip_pivot_growth, rcond;
682  SuperLU_gssvx(&m_sluOptions, &m_sluA,
683  m_q.data(), m_p.data(),
684  &m_sluEtree[0], &m_sluEqued,
685  &m_sluRscale[0], &m_sluCscale[0],
686  &m_sluL, &m_sluU,
687  NULL, 0,
688  &m_sluB, &m_sluX,
689  &recip_pivot_growth, &rcond,
690  &m_sluFerr[0], &m_sluBerr[0],
691  &m_sluStat, &info, Scalar());
692  StatFree(&m_sluStat);
693  m_info = info==0 ? Success : NumericalIssue;
694 }
695 
696 // the code of this extractData() function has been adapted from the SuperLU's Matlab support code,
697 //
698 // Copyright (c) 1994 by Xerox Corporation. All rights reserved.
699 //
700 // THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
701 // EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK.
702 //
703 template<typename MatrixType, typename Derived>
705 {
706  eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for extracting factors, you must first call either compute() or analyzePattern()/factorize()");
707  if (m_extractedDataAreDirty)
708  {
709  int upper;
710  int fsupc, istart, nsupr;
711  int lastl = 0, lastu = 0;
712  SCformat *Lstore = static_cast<SCformat*>(m_sluL.Store);
713  NCformat *Ustore = static_cast<NCformat*>(m_sluU.Store);
714  Scalar *SNptr;
715 
716  const int size = m_matrix.rows();
717  m_l.resize(size,size);
718  m_l.resizeNonZeros(Lstore->nnz);
719  m_u.resize(size,size);
720  m_u.resizeNonZeros(Ustore->nnz);
721 
722  int* Lcol = m_l.outerIndexPtr();
723  int* Lrow = m_l.innerIndexPtr();
724  Scalar* Lval = m_l.valuePtr();
725 
726  int* Ucol = m_u.outerIndexPtr();
727  int* Urow = m_u.innerIndexPtr();
728  Scalar* Uval = m_u.valuePtr();
729 
730  Ucol[0] = 0;
731  Ucol[0] = 0;
732 
733  /* for each supernode */
734  for (int k = 0; k <= Lstore->nsuper; ++k)
735  {
736  fsupc = L_FST_SUPC(k);
737  istart = L_SUB_START(fsupc);
738  nsupr = L_SUB_START(fsupc+1) - istart;
739  upper = 1;
740 
741  /* for each column in the supernode */
742  for (int j = fsupc; j < L_FST_SUPC(k+1); ++j)
743  {
744  SNptr = &((Scalar*)Lstore->nzval)[L_NZ_START(j)];
745 
746  /* Extract U */
747  for (int i = U_NZ_START(j); i < U_NZ_START(j+1); ++i)
748  {
749  Uval[lastu] = ((Scalar*)Ustore->nzval)[i];
750  /* Matlab doesn't like explicit zero. */
751  if (Uval[lastu] != 0.0)
752  Urow[lastu++] = U_SUB(i);
753  }
754  for (int i = 0; i < upper; ++i)
755  {
756  /* upper triangle in the supernode */
757  Uval[lastu] = SNptr[i];
758  /* Matlab doesn't like explicit zero. */
759  if (Uval[lastu] != 0.0)
760  Urow[lastu++] = L_SUB(istart+i);
761  }
762  Ucol[j+1] = lastu;
763 
764  /* Extract L */
765  Lval[lastl] = 1.0; /* unit diagonal */
766  Lrow[lastl++] = L_SUB(istart + upper - 1);
767  for (int i = upper; i < nsupr; ++i)
768  {
769  Lval[lastl] = SNptr[i];
770  /* Matlab doesn't like explicit zero. */
771  if (Lval[lastl] != 0.0)
772  Lrow[lastl++] = L_SUB(istart+i);
773  }
774  Lcol[j+1] = lastl;
775 
776  ++upper;
777  } /* for j ... */
778 
779  } /* for k ... */
780 
781  // squeeze the matrices :
782  m_l.resizeNonZeros(lastl);
783  m_u.resizeNonZeros(lastu);
784 
785  m_extractedDataAreDirty = false;
786  }
787 }
788 
789 template<typename MatrixType>
791 {
792  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()");
793 
794  if (m_extractedDataAreDirty)
795  this->extractData();
796 
797  Scalar det = Scalar(1);
798  for (int j=0; j<m_u.cols(); ++j)
799  {
800  if (m_u.outerIndexPtr()[j+1]-m_u.outerIndexPtr()[j] > 0)
801  {
802  int lastId = m_u.outerIndexPtr()[j+1]-1;
803  eigen_assert(m_u.innerIndexPtr()[lastId]<=j);
804  if (m_u.innerIndexPtr()[lastId]==j)
805  det *= m_u.valuePtr()[lastId];
806  }
807  }
808  if(m_sluEqued!='N')
809  return det/m_sluRscale.prod()/m_sluCscale.prod();
810  else
811  return det;
812 }
813 
814 #ifdef EIGEN_PARSED_BY_DOXYGEN
815 #define EIGEN_SUPERLU_HAS_ILU
816 #endif
817 
818 #ifdef EIGEN_SUPERLU_HAS_ILU
819 
834 template<typename _MatrixType>
835 class SuperILU : public SuperLUBase<_MatrixType,SuperILU<_MatrixType> >
836 {
837  public:
839  typedef _MatrixType MatrixType;
840  typedef typename Base::Scalar Scalar;
841  typedef typename Base::RealScalar RealScalar;
842  typedef typename Base::Index Index;
843 
844  public:
845 
846  SuperILU() : Base() { init(); }
847 
848  SuperILU(const MatrixType& matrix) : Base()
849  {
850  init();
851  compute(matrix);
852  }
853 
855  {
856  }
857 
864  void analyzePattern(const MatrixType& matrix)
865  {
866  Base::analyzePattern(matrix);
867  }
868 
875  void factorize(const MatrixType& matrix);
876 
877  #ifndef EIGEN_PARSED_BY_DOXYGEN
878 
879  template<typename Rhs,typename Dest>
880  void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const;
881  #endif // EIGEN_PARSED_BY_DOXYGEN
882 
883  protected:
884 
885  using Base::m_matrix;
886  using Base::m_sluOptions;
887  using Base::m_sluA;
888  using Base::m_sluB;
889  using Base::m_sluX;
890  using Base::m_p;
891  using Base::m_q;
892  using Base::m_sluEtree;
893  using Base::m_sluEqued;
894  using Base::m_sluRscale;
895  using Base::m_sluCscale;
896  using Base::m_sluL;
897  using Base::m_sluU;
898  using Base::m_sluStat;
899  using Base::m_sluFerr;
900  using Base::m_sluBerr;
901  using Base::m_l;
902  using Base::m_u;
903 
904  using Base::m_analysisIsOk;
907  using Base::m_isInitialized;
908  using Base::m_info;
909 
910  void init()
911  {
912  Base::init();
913 
914  ilu_set_default_options(&m_sluOptions);
915  m_sluOptions.PrintStat = NO;
916  m_sluOptions.ConditionNumber = NO;
917  m_sluOptions.Trans = NOTRANS;
918  m_sluOptions.ColPerm = MMD_AT_PLUS_A;
919 
920  // no attempt to preserve column sum
921  m_sluOptions.ILU_MILU = SILU;
922  // only basic ILU(k) support -- no direct control over memory consumption
923  // better to use ILU_DropRule = DROP_BASIC | DROP_AREA
924  // and set ILU_FillFactor to max memory growth
925  m_sluOptions.ILU_DropRule = DROP_BASIC;
927  }
928 
929  private:
930  SuperILU(SuperILU& ) { }
931 };
932 
933 template<typename MatrixType>
935 {
936  eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
937  if(!m_analysisIsOk)
938  {
939  m_info = InvalidInput;
940  return;
941  }
942 
943  this->initFactorization(a);
944 
945  int info = 0;
946  RealScalar recip_pivot_growth, rcond;
947 
948  StatInit(&m_sluStat);
949  SuperLU_gsisx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0],
950  &m_sluEqued, &m_sluRscale[0], &m_sluCscale[0],
951  &m_sluL, &m_sluU,
952  NULL, 0,
953  &m_sluB, &m_sluX,
954  &recip_pivot_growth, &rcond,
955  &m_sluStat, &info, Scalar());
956  StatFree(&m_sluStat);
957 
958  // FIXME how to better check for errors ???
959  m_info = info == 0 ? Success : NumericalIssue;
960  m_factorizationIsOk = true;
961 }
962 
963 template<typename MatrixType>
964 template<typename Rhs,typename Dest>
966 {
967  eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or analyzePattern()/factorize()");
968 
969  const int size = m_matrix.rows();
970  const int rhsCols = b.cols();
971  eigen_assert(size==b.rows());
972 
973  m_sluOptions.Trans = NOTRANS;
974  m_sluOptions.Fact = FACTORED;
975  m_sluOptions.IterRefine = NOREFINE;
976 
977  m_sluFerr.resize(rhsCols);
978  m_sluBerr.resize(rhsCols);
979  m_sluB = SluMatrix::Map(b.const_cast_derived());
980  m_sluX = SluMatrix::Map(x.derived());
981 
982  typename Rhs::PlainObject b_cpy;
983  if(m_sluEqued!='N')
984  {
985  b_cpy = b;
986  m_sluB = SluMatrix::Map(b_cpy.const_cast_derived());
987  }
988 
989  int info = 0;
990  RealScalar recip_pivot_growth, rcond;
991 
992  StatInit(&m_sluStat);
993  SuperLU_gsisx(&m_sluOptions, &m_sluA,
994  m_q.data(), m_p.data(),
995  &m_sluEtree[0], &m_sluEqued,
996  &m_sluRscale[0], &m_sluCscale[0],
997  &m_sluL, &m_sluU,
998  NULL, 0,
999  &m_sluB, &m_sluX,
1000  &recip_pivot_growth, &rcond,
1001  &m_sluStat, &info, Scalar());
1002  StatFree(&m_sluStat);
1003 
1004  m_info = info==0 ? Success : NumericalIssue;
1005 }
1006 #endif
1007 
1008 namespace internal {
1009 
1010 template<typename _MatrixType, typename Derived, typename Rhs>
1011 struct solve_retval<SuperLUBase<_MatrixType,Derived>, Rhs>
1012  : solve_retval_base<SuperLUBase<_MatrixType,Derived>, Rhs>
1013 {
1014  typedef SuperLUBase<_MatrixType,Derived> Dec;
1015  EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
1016 
1017  template<typename Dest> void evalTo(Dest& dst) const
1018  {
1019  dec().derived()._solve(rhs(),dst);
1020  }
1021 };
1022 
1023 template<typename _MatrixType, typename Derived, typename Rhs>
1024 struct sparse_solve_retval<SuperLUBase<_MatrixType,Derived>, Rhs>
1025  : sparse_solve_retval_base<SuperLUBase<_MatrixType,Derived>, Rhs>
1026 {
1027  typedef SuperLUBase<_MatrixType,Derived> Dec;
1029 
1030  template<typename Dest> void evalTo(Dest& dst) const
1031  {
1032  dec().derived()._solve(rhs(),dst);
1033  }
1034 };
1035 
1036 } // end namespace internal
1037 
1038 } // end namespace Eigen
1039 
1040 #endif // EIGEN_SUPERLUSUPPORT_H