10 #ifndef EIGEN_SUPERLUSUPPORT_H
11 #define EIGEN_SUPERLUSUPPORT_H
15 #if defined(SUPERLU_MAJOR_VERSION) && (SUPERLU_MAJOR_VERSION >= 5)
16 #define DECL_GSSVX(PREFIX,FLOATTYPE,KEYTYPE) \
18 extern void PREFIX##gssvx(superlu_options_t *, SuperMatrix *, int *, int *, int *, \
19 char *, FLOATTYPE *, FLOATTYPE *, SuperMatrix *, SuperMatrix *, \
20 void *, int, SuperMatrix *, SuperMatrix *, \
21 FLOATTYPE *, FLOATTYPE *, FLOATTYPE *, FLOATTYPE *, \
22 GlobalLU_t *, mem_usage_t *, SuperLUStat_t *, int *); \
24 inline float SuperLU_gssvx(superlu_options_t *options, SuperMatrix *A, \
25 int *perm_c, int *perm_r, int *etree, char *equed, \
26 FLOATTYPE *R, FLOATTYPE *C, SuperMatrix *L, \
27 SuperMatrix *U, void *work, int lwork, \
28 SuperMatrix *B, SuperMatrix *X, \
29 FLOATTYPE *recip_pivot_growth, \
30 FLOATTYPE *rcond, FLOATTYPE *ferr, FLOATTYPE *berr, \
31 SuperLUStat_t *stats, int *info, KEYTYPE) { \
32 mem_usage_t mem_usage; \
34 PREFIX##gssvx(options, A, perm_c, perm_r, etree, equed, R, C, L, \
35 U, work, lwork, B, X, recip_pivot_growth, rcond, \
36 ferr, berr, &gLU, &mem_usage, stats, info); \
37 return mem_usage.for_lu; \
39 #else // version < 5.0
40 #define DECL_GSSVX(PREFIX,FLOATTYPE,KEYTYPE) \
42 extern void PREFIX##gssvx(superlu_options_t *, SuperMatrix *, int *, int *, int *, \
43 char *, FLOATTYPE *, FLOATTYPE *, SuperMatrix *, SuperMatrix *, \
44 void *, int, SuperMatrix *, SuperMatrix *, \
45 FLOATTYPE *, FLOATTYPE *, FLOATTYPE *, FLOATTYPE *, \
46 mem_usage_t *, SuperLUStat_t *, int *); \
48 inline float SuperLU_gssvx(superlu_options_t *options, SuperMatrix *A, \
49 int *perm_c, int *perm_r, int *etree, char *equed, \
50 FLOATTYPE *R, FLOATTYPE *C, SuperMatrix *L, \
51 SuperMatrix *U, void *work, int lwork, \
52 SuperMatrix *B, SuperMatrix *X, \
53 FLOATTYPE *recip_pivot_growth, \
54 FLOATTYPE *rcond, FLOATTYPE *ferr, FLOATTYPE *berr, \
55 SuperLUStat_t *stats, int *info, KEYTYPE) { \
56 mem_usage_t mem_usage; \
57 PREFIX##gssvx(options, A, perm_c, perm_r, etree, equed, R, C, L, \
58 U, work, lwork, B, X, recip_pivot_growth, rcond, \
59 ferr, berr, &mem_usage, stats, info); \
60 return mem_usage.for_lu; \
70 #define EIGEN_SUPERLU_HAS_ILU
73 #ifdef EIGEN_SUPERLU_HAS_ILU
76 #define DECL_GSISX(PREFIX,FLOATTYPE,KEYTYPE) \
78 extern void PREFIX##gsisx(superlu_options_t *, SuperMatrix *, int *, int *, int *, \
79 char *, FLOATTYPE *, FLOATTYPE *, SuperMatrix *, SuperMatrix *, \
80 void *, int, SuperMatrix *, SuperMatrix *, FLOATTYPE *, FLOATTYPE *, \
81 mem_usage_t *, SuperLUStat_t *, int *); \
83 inline float SuperLU_gsisx(superlu_options_t *options, SuperMatrix *A, \
84 int *perm_c, int *perm_r, int *etree, char *equed, \
85 FLOATTYPE *R, FLOATTYPE *C, SuperMatrix *L, \
86 SuperMatrix *U, void *work, int lwork, \
87 SuperMatrix *B, SuperMatrix *X, \
88 FLOATTYPE *recip_pivot_growth, \
90 SuperLUStat_t *stats, int *info, KEYTYPE) { \
91 mem_usage_t mem_usage; \
92 PREFIX##gsisx(options, A, perm_c, perm_r, etree, equed, R, C, L, \
93 U, work, lwork, B, X, recip_pivot_growth, rcond, \
94 &mem_usage, stats, info); \
95 return mem_usage.for_lu; \
98 DECL_GSISX(
s,
float,
float)
99 DECL_GSISX(
c,
float,std::complex<float>)
100 DECL_GSISX(
d,
double,
double)
101 DECL_GSISX(
z,
double,std::complex<double>)
105 template<
typename MatrixType>
131 SuperMatrix::operator=(
static_cast<const SuperMatrix&
>(
other));
148 if (
t==SLU_NC ||
t==SLU_NR ||
t==SLU_DN)
157 template<
typename Scalar>
170 eigen_assert(
false &&
"Scalar type not supported by SuperLU");
174 template<
typename MatrixType>
180 res.setStorageType(SLU_DN);
184 res.nrow = internal::convert_index<int>(
mat.rows());
185 res.ncol = internal::convert_index<int>(
mat.cols());
187 res.storage.lda = internal::convert_index<int>(MatrixType::IsVectorAtCompileTime ?
mat.size() :
mat.outerStride());
188 res.storage.values = (
void*)(
mat.data());
192 template<
typename MatrixType>
199 res.setStorageType(SLU_NR);
200 res.nrow = internal::convert_index<int>(
mat.cols());
201 res.ncol = internal::convert_index<int>(
mat.rows());
205 res.setStorageType(SLU_NC);
206 res.nrow = internal::convert_index<int>(
mat.rows());
207 res.ncol = internal::convert_index<int>(
mat.cols());
212 res.storage.nnz = internal::convert_index<int>(
mat.nonZeros());
213 res.storage.values =
mat.valuePtr();
214 res.storage.innerInd =
mat.innerIndexPtr();
215 res.storage.outerInd =
mat.outerIndexPtr();
220 if (
int(MatrixType::Flags) &
int(
Upper))
225 eigen_assert(((
int(MatrixType::Flags) &
int(
SelfAdjoint))==0) &&
"SelfAdjoint matrix shape not supported by SuperLU");
231 template<
typename Scalar,
int Rows,
int Cols,
int Options,
int MRows,
int MCols>
238 res.setStorageType(SLU_DN);
245 res.storage.lda =
mat.outerStride();
246 res.storage.values =
mat.data();
250 template<
typename Derived>
258 res.setStorageType(SLU_NR);
264 res.setStorageType(SLU_NC);
271 res.storage.nnz =
mat.nonZeros();
272 res.storage.values =
mat.valuePtr();
273 res.storage.innerInd =
mat.innerIndexPtr();
274 res.storage.outerInd =
mat.outerIndexPtr();
279 if (MatrixType::Flags &
Upper)
281 if (MatrixType::Flags &
Lower)
290 template<
typename MatrixType>
297 template<
typename Scalar,
int Flags,
typename Index>
316 template<
typename _MatrixType,
typename Derived>
385 template<
typename Stream>
433 Destroy_SuperNode_Matrix(&
m_sluL);
435 Destroy_CompCol_Matrix(&
m_sluU);
487 template<
typename _MatrixType>
540 template<
typename Rhs,
typename Dest>
612 template<
typename MatrixType>
615 eigen_assert(m_analysisIsOk &&
"You must first call analyzePattern()");
622 this->initFactorization(
a);
624 m_sluOptions.ColPerm = COLAMD;
629 StatInit(&m_sluStat);
630 SuperLU_gssvx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0],
631 &m_sluEqued, &m_sluRscale[0], &m_sluCscale[0],
635 &recip_pivot_growth, &rcond,
638 StatFree(&m_sluStat);
640 m_extractedDataAreDirty =
true;
644 m_factorizationIsOk =
true;
647 template<
typename MatrixType>
648 template<
typename Rhs,
typename Dest>
651 eigen_assert(m_factorizationIsOk &&
"The decomposition is not in a valid state for solving, you must first call either compute() or analyzePattern()/factorize()");
653 const Index rhsCols =
b.cols();
656 m_sluOptions.Trans = NOTRANS;
657 m_sluOptions.Fact = FACTORED;
658 m_sluOptions.IterRefine = NOREFINE;
661 m_sluFerr.resize(rhsCols);
662 m_sluBerr.resize(rhsCols);
670 typename Rhs::PlainObject b_cpy;
677 StatInit(&m_sluStat);
680 SuperLU_gssvx(&m_sluOptions, &m_sluA,
681 m_q.data(), m_p.data(),
682 &m_sluEtree[0], &m_sluEqued,
683 &m_sluRscale[0], &m_sluCscale[0],
687 &recip_pivot_growth, &rcond,
688 &m_sluFerr[0], &m_sluBerr[0],
690 StatFree(&m_sluStat);
692 if(
x.derived().data() != x_ref.data())
705 template<
typename MatrixType,
typename Derived>
708 eigen_assert(m_factorizationIsOk &&
"The decomposition is not in a valid state for extracting factors, you must first call either compute() or analyzePattern()/factorize()");
709 if (m_extractedDataAreDirty)
712 int fsupc, istart, nsupr;
713 int lastl = 0, lastu = 0;
714 SCformat *Lstore =
static_cast<SCformat*
>(m_sluL.Store);
715 NCformat *Ustore =
static_cast<NCformat*
>(m_sluU.Store);
720 m_l.resizeNonZeros(Lstore->nnz);
722 m_u.resizeNonZeros(Ustore->nnz);
724 int* Lcol = m_l.outerIndexPtr();
725 int* Lrow = m_l.innerIndexPtr();
726 Scalar* Lval = m_l.valuePtr();
728 int* Ucol = m_u.outerIndexPtr();
729 int* Urow = m_u.innerIndexPtr();
730 Scalar* Uval = m_u.valuePtr();
736 for (
int k = 0; k <= Lstore->nsuper; ++k)
738 fsupc = L_FST_SUPC(k);
739 istart = L_SUB_START(fsupc);
740 nsupr = L_SUB_START(fsupc+1) - istart;
744 for (
int j = fsupc;
j < L_FST_SUPC(k+1); ++
j)
746 SNptr = &((
Scalar*)Lstore->nzval)[L_NZ_START(
j)];
749 for (
int i = U_NZ_START(
j);
i < U_NZ_START(
j+1); ++
i)
751 Uval[lastu] = ((
Scalar*)Ustore->nzval)[
i];
753 if (Uval[lastu] != 0.0)
754 Urow[lastu++] = U_SUB(
i);
756 for (
int i = 0;
i < upper; ++
i)
759 Uval[lastu] = SNptr[
i];
761 if (Uval[lastu] != 0.0)
762 Urow[lastu++] = L_SUB(istart+
i);
768 Lrow[lastl++] = L_SUB(istart + upper - 1);
769 for (
int i = upper;
i < nsupr; ++
i)
771 Lval[lastl] = SNptr[
i];
773 if (Lval[lastl] != 0.0)
774 Lrow[lastl++] = L_SUB(istart+
i);
784 m_l.resizeNonZeros(lastl);
785 m_u.resizeNonZeros(lastu);
787 m_extractedDataAreDirty =
false;
791 template<
typename MatrixType>
794 eigen_assert(m_factorizationIsOk &&
"The decomposition is not in a valid state for computing the determinant, you must first call either compute() or analyzePattern()/factorize()");
796 if (m_extractedDataAreDirty)
800 for (
int j=0;
j<m_u.cols(); ++
j)
802 if (m_u.outerIndexPtr()[
j+1]-m_u.outerIndexPtr()[
j] > 0)
804 int lastId = m_u.outerIndexPtr()[
j+1]-1;
806 if (m_u.innerIndexPtr()[lastId]==
j)
807 det *= m_u.valuePtr()[lastId];
813 return det/m_sluRscale.prod()/m_sluCscale.prod();
818 #ifdef EIGEN_PARSED_BY_DOXYGEN
819 #define EIGEN_SUPERLU_HAS_ILU
822 #ifdef EIGEN_SUPERLU_HAS_ILU
840 template<
typename _MatrixType>
841 class SuperILU :
public SuperLUBase<_MatrixType,SuperILU<_MatrixType> >
850 using Base::_solve_impl;
872 Base::analyzePattern(
matrix);
883 #ifndef EIGEN_PARSED_BY_DOXYGEN
885 template<
typename Rhs,
typename Dest>
886 void _solve_impl(
const MatrixBase<Rhs> &
b, MatrixBase<Dest> &dest)
const;
887 #endif // EIGEN_PARSED_BY_DOXYGEN
891 using Base::m_matrix;
892 using Base::m_sluOptions;
898 using Base::m_sluEtree;
899 using Base::m_sluEqued;
900 using Base::m_sluRscale;
901 using Base::m_sluCscale;
904 using Base::m_sluStat;
905 using Base::m_sluFerr;
906 using Base::m_sluBerr;
910 using Base::m_analysisIsOk;
911 using Base::m_factorizationIsOk;
912 using Base::m_extractedDataAreDirty;
913 using Base::m_isInitialized;
920 ilu_set_default_options(&m_sluOptions);
921 m_sluOptions.PrintStat = NO;
922 m_sluOptions.ConditionNumber = NO;
923 m_sluOptions.Trans = NOTRANS;
924 m_sluOptions.ColPerm = MMD_AT_PLUS_A;
927 m_sluOptions.ILU_MILU = SILU;
931 m_sluOptions.ILU_DropRule = DROP_BASIC;
932 m_sluOptions.ILU_DropTol = NumTraits<Scalar>::dummy_precision()*10;
936 SuperILU(SuperILU& ) { }
939 template<
typename MatrixType>
940 void SuperILU<MatrixType>::factorize(
const MatrixType&
a)
942 eigen_assert(m_analysisIsOk &&
"You must first call analyzePattern()");
949 this->initFactorization(
a);
954 StatInit(&m_sluStat);
955 SuperLU_gsisx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0],
956 &m_sluEqued, &m_sluRscale[0], &m_sluCscale[0],
960 &recip_pivot_growth, &rcond,
962 StatFree(&m_sluStat);
966 m_factorizationIsOk =
true;
969 #ifndef EIGEN_PARSED_BY_DOXYGEN
970 template<
typename MatrixType>
971 template<
typename Rhs,
typename Dest>
972 void SuperILU<MatrixType>::_solve_impl(
const MatrixBase<Rhs> &
b, MatrixBase<Dest>&
x)
const
974 eigen_assert(m_factorizationIsOk &&
"The decomposition is not in a valid state for solving, you must first call either compute() or analyzePattern()/factorize()");
976 const int rhsCols =
b.cols();
979 m_sluOptions.Trans = NOTRANS;
980 m_sluOptions.Fact = FACTORED;
981 m_sluOptions.IterRefine = NOREFINE;
983 m_sluFerr.resize(rhsCols);
984 m_sluBerr.resize(rhsCols);
986 Ref<const Matrix<typename Rhs::Scalar,Dynamic,Dynamic,ColMajor> > b_ref(
b);
987 Ref<const Matrix<typename Dest::Scalar,Dynamic,Dynamic,ColMajor> > x_ref(
x);
992 typename Rhs::PlainObject b_cpy;
1002 StatInit(&m_sluStat);
1003 SuperLU_gsisx(&m_sluOptions, &m_sluA,
1004 m_q.data(), m_p.data(),
1005 &m_sluEtree[0], &m_sluEqued,
1006 &m_sluRscale[0], &m_sluCscale[0],
1010 &recip_pivot_growth, &rcond,
1012 StatFree(&m_sluStat);
1014 if(
x.derived().data() != x_ref.data())
1025 #endif // EIGEN_SUPERLUSUPPORT_H