10 #ifndef EIGEN_SUPERLUSUPPORT_H
11 #define EIGEN_SUPERLUSUPPORT_H
15 #define DECL_GSSVX(PREFIX,FLOATTYPE,KEYTYPE) \
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 *); \
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; \
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>)
45 #define EIGEN_SUPERLU_HAS_ILU
48 #ifdef EIGEN_SUPERLU_HAS_ILU
51 #define DECL_GSISX(PREFIX,FLOATTYPE,KEYTYPE) \
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 *); \
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, \
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; \
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>)
80 template<
typename MatrixType>
81 struct SluMatrixMapHelper;
90 struct SluMatrix : SuperMatrix
97 SluMatrix(
const SluMatrix& other)
101 storage = other.storage;
104 SluMatrix& operator=(
const SluMatrix& other)
106 SuperMatrix::operator=(static_cast<const SuperMatrix&>(other));
108 storage = other.storage;
114 union {
int nnz;
int lda;};
120 void setStorageType(Stype_t t)
123 if (t==SLU_NC || t==SLU_NR || t==SLU_DN)
127 eigen_assert(
false &&
"storage type not supported");
132 template<
typename Scalar>
135 if (internal::is_same<Scalar,float>::value)
137 else if (internal::is_same<Scalar,double>::value)
139 else if (internal::is_same<Scalar,std::complex<float> >::value)
141 else if (internal::is_same<Scalar,std::complex<double> >::value)
145 eigen_assert(
false &&
"Scalar type not supported by SuperLU");
149 template<
typename MatrixType>
150 static SluMatrix Map(MatrixBase<MatrixType>& _mat)
152 MatrixType& mat(_mat.derived());
153 eigen_assert( ((MatrixType::Flags&
RowMajorBit)!=RowMajorBit) &&
"row-major dense matrices are not supported by SuperLU");
155 res.setStorageType(SLU_DN);
156 res.setScalarType<
typename MatrixType::Scalar>();
159 res.nrow = mat.rows();
160 res.ncol = mat.cols();
162 res.storage.lda = MatrixType::IsVectorAtCompileTime ? mat.size() : mat.outerStride();
163 res.storage.values = mat.data();
167 template<
typename MatrixType>
168 static SluMatrix Map(SparseMatrixBase<MatrixType>& mat)
173 res.setStorageType(SLU_NR);
174 res.nrow = mat.cols();
175 res.ncol = mat.rows();
179 res.setStorageType(SLU_NC);
180 res.nrow = mat.rows();
181 res.ncol = mat.cols();
186 res.storage.nnz = mat.nonZeros();
187 res.storage.values = mat.derived().valuePtr();
188 res.storage.innerInd = mat.derived().innerIndexPtr();
189 res.storage.outerInd = mat.derived().outerIndexPtr();
191 res.setScalarType<
typename MatrixType::Scalar>();
194 if (MatrixType::Flags &
Upper)
196 if (MatrixType::Flags &
Lower)
199 eigen_assert(((MatrixType::Flags &
SelfAdjoint)==0) &&
"SelfAdjoint matrix shape not supported by SuperLU");
205 template<
typename Scalar,
int Rows,
int Cols,
int Options,
int MRows,
int MCols>
206 struct SluMatrixMapHelper<Matrix<Scalar,Rows,Cols,Options,MRows,MCols> >
208 typedef Matrix<Scalar,Rows,Cols,Options,MRows,MCols> MatrixType;
209 static void run(MatrixType& mat, SluMatrix& res)
211 eigen_assert( ((Options&
RowMajor)!=RowMajor) &&
"row-major dense matrices is not supported by SuperLU");
212 res.setStorageType(SLU_DN);
213 res.setScalarType<Scalar>();
216 res.nrow = mat.rows();
217 res.ncol = mat.cols();
219 res.storage.lda = mat.outerStride();
220 res.storage.values = mat.data();
224 template<
typename Derived>
225 struct SluMatrixMapHelper<SparseMatrixBase<Derived> >
227 typedef Derived MatrixType;
228 static void run(MatrixType& mat, SluMatrix& res)
230 if ((MatrixType::Flags&RowMajorBit)==RowMajorBit)
232 res.setStorageType(SLU_NR);
233 res.nrow = mat.cols();
234 res.ncol = mat.rows();
238 res.setStorageType(SLU_NC);
239 res.nrow = mat.rows();
240 res.ncol = mat.cols();
245 res.storage.nnz = mat.nonZeros();
246 res.storage.values = mat.valuePtr();
247 res.storage.innerInd = mat.innerIndexPtr();
248 res.storage.outerInd = mat.outerIndexPtr();
250 res.setScalarType<
typename MatrixType::Scalar>();
253 if (MatrixType::Flags & Upper)
255 if (MatrixType::Flags & Lower)
258 eigen_assert(((MatrixType::Flags &
SelfAdjoint)==0) &&
"SelfAdjoint matrix shape not supported by SuperLU");
264 template<
typename MatrixType>
265 SluMatrix asSluMatrix(MatrixType& mat)
267 return SluMatrix::Map(mat);
271 template<
typename Scalar,
int Flags,
typename Index>
274 eigen_assert((Flags&
RowMajor)==RowMajor && sluMat.Stype == SLU_NR
275 || (Flags&
ColMajor)==ColMajor && sluMat.Stype == SLU_NC);
277 Index outerSize = (Flags&
RowMajor)==RowMajor ? sluMat.ncol : sluMat.nrow;
280 sluMat.nrow, sluMat.ncol, sluMat.storage.outerInd[outerSize],
281 sluMat.storage.outerInd, sluMat.storage.innerInd, reinterpret_cast<Scalar*>(sluMat.storage.values) );
290 template<
typename _MatrixType,
typename Derived>
294 typedef _MatrixType MatrixType;
295 typedef typename MatrixType::Scalar Scalar;
296 typedef typename MatrixType::RealScalar RealScalar;
297 typedef typename MatrixType::Index Index;
312 Derived& derived() {
return *
static_cast<Derived*
>(
this); }
313 const Derived& derived()
const {
return *
static_cast<const Derived*
>(
this); }
315 inline Index rows()
const {
return m_matrix.
rows(); }
316 inline Index cols()
const {
return m_matrix.
cols(); }
319 inline superlu_options_t&
options() {
return m_sluOptions; }
328 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
335 derived().analyzePattern(matrix);
336 derived().factorize(matrix);
343 template<
typename Rhs>
346 eigen_assert(m_isInitialized &&
"SuperLU is not initialized.");
347 eigen_assert(rows()==b.rows()
348 &&
"SuperLU::solve(): invalid number of rows of the right hand side matrix b");
349 return internal::solve_retval<SuperLUBase, Rhs>(*
this, b.derived());
373 m_isInitialized =
true;
375 m_analysisIsOk =
true;
376 m_factorizationIsOk =
false;
379 template<
typename Stream>
380 void dumpMemory(Stream& s)
385 void initFactorization(
const MatrixType& a)
387 set_default_options(&this->m_sluOptions);
389 const int size = a.rows();
392 m_sluA = internal::asSluMatrix(m_matrix);
399 m_sluEtree.resize(size);
402 m_sluB.setStorageType(SLU_DN);
403 m_sluB.setScalarType<Scalar>();
404 m_sluB.Mtype = SLU_GE;
405 m_sluB.storage.values = 0;
408 m_sluB.storage.lda = size;
411 m_extractedDataAreDirty =
true;
417 m_isInitialized =
false;
422 void extractData()
const;
427 Destroy_SuperNode_Matrix(&m_sluL);
429 Destroy_CompCol_Matrix(&m_sluU);
434 memset(&m_sluL,0,
sizeof m_sluL);
435 memset(&m_sluU,0,
sizeof m_sluU);
439 mutable LUMatrixType m_l;
440 mutable LUMatrixType m_u;
441 mutable IntColVectorType m_p;
442 mutable IntRowVectorType m_q;
444 mutable LUMatrixType m_matrix;
445 mutable SluMatrix m_sluA;
446 mutable SuperMatrix m_sluL, m_sluU;
447 mutable SluMatrix m_sluB, m_sluX;
448 mutable SuperLUStat_t m_sluStat;
449 mutable superlu_options_t m_sluOptions;
450 mutable std::vector<int> m_sluEtree;
451 mutable Matrix<RealScalar,Dynamic,1> m_sluRscale, m_sluCscale;
452 mutable Matrix<RealScalar,Dynamic,1> m_sluFerr, m_sluBerr;
453 mutable char m_sluEqued;
456 bool m_isInitialized;
457 int m_factorizationIsOk;
459 mutable bool m_extractedDataAreDirty;
462 SuperLUBase(SuperLUBase& ) { }
478 template<
typename _MatrixType>
483 typedef _MatrixType MatrixType;
484 typedef typename Base::Scalar Scalar;
485 typedef typename Base::RealScalar RealScalar;
486 typedef typename Base::Index Index;
516 m_isInitialized =
false;
526 void factorize(
const MatrixType& matrix);
528 #ifndef EIGEN_PARSED_BY_DOXYGEN
530 template<
typename Rhs,
typename Dest>
532 #endif // EIGEN_PARSED_BY_DOXYGEN
534 inline const LMatrixType& matrixL()
const
536 if (m_extractedDataAreDirty) this->extractData();
540 inline const UMatrixType& matrixU()
const
542 if (m_extractedDataAreDirty) this->extractData();
546 inline const IntColVectorType& permutationP()
const
548 if (m_extractedDataAreDirty) this->extractData();
552 inline const IntRowVectorType& permutationQ()
const
554 if (m_extractedDataAreDirty) this->extractData();
558 Scalar determinant()
const;
562 using Base::m_matrix;
563 using Base::m_sluOptions;
569 using Base::m_sluEtree;
570 using Base::m_sluEqued;
571 using Base::m_sluRscale;
572 using Base::m_sluCscale;
575 using Base::m_sluStat;
576 using Base::m_sluFerr;
577 using Base::m_sluBerr;
581 using Base::m_analysisIsOk;
582 using Base::m_factorizationIsOk;
583 using Base::m_extractedDataAreDirty;
584 using Base::m_isInitialized;
591 set_default_options(&this->m_sluOptions);
592 m_sluOptions.PrintStat = NO;
593 m_sluOptions.ConditionNumber = NO;
594 m_sluOptions.Trans = NOTRANS;
595 m_sluOptions.ColPerm = COLAMD;
600 SuperLU(SuperLU& ) { }
603 template<
typename MatrixType>
606 eigen_assert(m_analysisIsOk &&
"You must first call analyzePattern()");
613 this->initFactorization(a);
616 RealScalar recip_pivot_growth, rcond;
617 RealScalar ferr, berr;
619 StatInit(&m_sluStat);
620 SuperLU_gssvx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0],
621 &m_sluEqued, &m_sluRscale[0], &m_sluCscale[0],
625 &recip_pivot_growth, &rcond,
627 &m_sluStat, &info, Scalar());
628 StatFree(&m_sluStat);
630 m_extractedDataAreDirty =
true;
634 m_factorizationIsOk =
true;
637 template<
typename MatrixType>
638 template<
typename Rhs,
typename Dest>
641 eigen_assert(m_factorizationIsOk &&
"The decomposition is not in a valid state for solving, you must first call either compute() or analyzePattern()/factorize()");
643 const int size = m_matrix.rows();
644 const int rhsCols = b.cols();
645 eigen_assert(size==b.rows());
647 m_sluOptions.Trans = NOTRANS;
648 m_sluOptions.Fact = FACTORED;
649 m_sluOptions.IterRefine = NOREFINE;
652 m_sluFerr.resize(rhsCols);
653 m_sluBerr.resize(rhsCols);
654 m_sluB = SluMatrix::Map(b.const_cast_derived());
655 m_sluX = SluMatrix::Map(x.derived());
657 typename Rhs::PlainObject b_cpy;
661 m_sluB = SluMatrix::Map(b_cpy.const_cast_derived());
664 StatInit(&m_sluStat);
666 RealScalar recip_pivot_growth, rcond;
667 SuperLU_gssvx(&m_sluOptions, &m_sluA,
668 m_q.data(), m_p.data(),
669 &m_sluEtree[0], &m_sluEqued,
670 &m_sluRscale[0], &m_sluCscale[0],
674 &recip_pivot_growth, &rcond,
675 &m_sluFerr[0], &m_sluBerr[0],
676 &m_sluStat, &info, Scalar());
677 StatFree(&m_sluStat);
688 template<
typename MatrixType,
typename Derived>
689 void SuperLUBase<MatrixType,Derived>::extractData()
const
691 eigen_assert(m_factorizationIsOk &&
"The decomposition is not in a valid state for extracting factors, you must first call either compute() or analyzePattern()/factorize()");
692 if (m_extractedDataAreDirty)
695 int fsupc, istart, nsupr;
696 int lastl = 0, lastu = 0;
697 SCformat *Lstore =
static_cast<SCformat*
>(m_sluL.Store);
698 NCformat *Ustore =
static_cast<NCformat*
>(m_sluU.Store);
701 const int size = m_matrix.rows();
702 m_l.resize(size,size);
703 m_l.resizeNonZeros(Lstore->nnz);
704 m_u.resize(size,size);
705 m_u.resizeNonZeros(Ustore->nnz);
707 int* Lcol = m_l.outerIndexPtr();
708 int* Lrow = m_l.innerIndexPtr();
709 Scalar* Lval = m_l.valuePtr();
711 int* Ucol = m_u.outerIndexPtr();
712 int* Urow = m_u.innerIndexPtr();
713 Scalar* Uval = m_u.valuePtr();
719 for (
int k = 0; k <= Lstore->nsuper; ++k)
721 fsupc = L_FST_SUPC(k);
722 istart = L_SUB_START(fsupc);
723 nsupr = L_SUB_START(fsupc+1) - istart;
727 for (
int j = fsupc; j < L_FST_SUPC(k+1); ++j)
729 SNptr = &((Scalar*)Lstore->nzval)[L_NZ_START(j)];
732 for (
int i = U_NZ_START(j); i < U_NZ_START(j+1); ++i)
734 Uval[lastu] = ((Scalar*)Ustore->nzval)[i];
736 if (Uval[lastu] != 0.0)
737 Urow[lastu++] = U_SUB(i);
739 for (
int i = 0; i < upper; ++i)
742 Uval[lastu] = SNptr[i];
744 if (Uval[lastu] != 0.0)
745 Urow[lastu++] = L_SUB(istart+i);
751 Lrow[lastl++] = L_SUB(istart + upper - 1);
752 for (
int i = upper; i < nsupr; ++i)
754 Lval[lastl] = SNptr[i];
756 if (Lval[lastl] != 0.0)
757 Lrow[lastl++] = L_SUB(istart+i);
767 m_l.resizeNonZeros(lastl);
768 m_u.resizeNonZeros(lastu);
770 m_extractedDataAreDirty =
false;
774 template<
typename MatrixType>
775 typename SuperLU<MatrixType>::Scalar SuperLU<MatrixType>::determinant()
const
777 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()");
779 if (m_extractedDataAreDirty)
782 Scalar det = Scalar(1);
783 for (
int j=0; j<m_u.cols(); ++j)
785 if (m_u.outerIndexPtr()[j+1]-m_u.outerIndexPtr()[j] > 0)
787 int lastId = m_u.outerIndexPtr()[j+1]-1;
788 eigen_assert(m_u.innerIndexPtr()[lastId]<=j);
789 if (m_u.innerIndexPtr()[lastId]==j)
790 det *= m_u.valuePtr()[lastId];
794 return det/m_sluRscale.prod()/m_sluCscale.prod();
799 #ifdef EIGEN_PARSED_BY_DOXYGEN
800 #define EIGEN_SUPERLU_HAS_ILU
803 #ifdef EIGEN_SUPERLU_HAS_ILU
819 template<
typename _MatrixType>
824 typedef _MatrixType MatrixType;
825 typedef typename Base::Scalar Scalar;
826 typedef typename Base::RealScalar RealScalar;
827 typedef typename Base::Index Index;
860 void factorize(
const MatrixType& matrix);
862 #ifndef EIGEN_PARSED_BY_DOXYGEN
864 template<
typename Rhs,
typename Dest>
866 #endif // EIGEN_PARSED_BY_DOXYGEN
870 using Base::m_matrix;
871 using Base::m_sluOptions;
877 using Base::m_sluEtree;
878 using Base::m_sluEqued;
879 using Base::m_sluRscale;
880 using Base::m_sluCscale;
883 using Base::m_sluStat;
884 using Base::m_sluFerr;
885 using Base::m_sluBerr;
889 using Base::m_analysisIsOk;
890 using Base::m_factorizationIsOk;
891 using Base::m_extractedDataAreDirty;
892 using Base::m_isInitialized;
899 ilu_set_default_options(&m_sluOptions);
900 m_sluOptions.PrintStat = NO;
901 m_sluOptions.ConditionNumber = NO;
902 m_sluOptions.Trans = NOTRANS;
903 m_sluOptions.ColPerm = MMD_AT_PLUS_A;
906 m_sluOptions.ILU_MILU = SILU;
910 m_sluOptions.ILU_DropRule = DROP_BASIC;
915 SuperILU(SuperILU& ) { }
918 template<
typename MatrixType>
921 eigen_assert(m_analysisIsOk &&
"You must first call analyzePattern()");
928 this->initFactorization(a);
931 RealScalar recip_pivot_growth, rcond;
933 StatInit(&m_sluStat);
934 SuperLU_gsisx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0],
935 &m_sluEqued, &m_sluRscale[0], &m_sluCscale[0],
939 &recip_pivot_growth, &rcond,
940 &m_sluStat, &info, Scalar());
941 StatFree(&m_sluStat);
945 m_factorizationIsOk =
true;
948 template<
typename MatrixType>
949 template<
typename Rhs,
typename Dest>
952 eigen_assert(m_factorizationIsOk &&
"The decomposition is not in a valid state for solving, you must first call either compute() or analyzePattern()/factorize()");
954 const int size = m_matrix.rows();
955 const int rhsCols = b.cols();
956 eigen_assert(size==b.rows());
958 m_sluOptions.Trans = NOTRANS;
959 m_sluOptions.Fact = FACTORED;
960 m_sluOptions.IterRefine = NOREFINE;
962 m_sluFerr.resize(rhsCols);
963 m_sluBerr.resize(rhsCols);
964 m_sluB = SluMatrix::Map(b.const_cast_derived());
965 m_sluX = SluMatrix::Map(x.derived());
967 typename Rhs::PlainObject b_cpy;
971 m_sluB = SluMatrix::Map(b_cpy.const_cast_derived());
975 RealScalar recip_pivot_growth, rcond;
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],
985 &recip_pivot_growth, &rcond,
986 &m_sluStat, &info, Scalar());
987 StatFree(&m_sluStat);
995 template<
typename _MatrixType,
typename Derived,
typename Rhs>
996 struct solve_retval<SuperLUBase<_MatrixType,Derived>, Rhs>
997 : solve_retval_base<SuperLUBase<_MatrixType,Derived>, Rhs>
999 typedef SuperLUBase<_MatrixType,Derived> Dec;
1000 EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
1002 template<typename Dest>
void evalTo(Dest& dst)
const
1004 dec().derived()._solve(rhs(),dst);
1008 template<
typename _MatrixType,
typename Derived,
typename Rhs>
1009 struct sparse_solve_retval<SuperLUBase<_MatrixType,Derived>, Rhs>
1010 : sparse_solve_retval_base<SuperLUBase<_MatrixType,Derived>, Rhs>
1012 typedef SuperLUBase<_MatrixType,Derived> Dec;
1013 EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
1015 template<typename Dest>
void evalTo(Dest& dst)
const
1017 dec().derived()._solve(rhs(),dst);
1025 #endif // EIGEN_SUPERLUSUPPORT_H