25 #ifndef EIGEN_SUPERLUSUPPORT_H
26 #define EIGEN_SUPERLUSUPPORT_H
30 #define DECL_GSSVX(PREFIX,FLOATTYPE,KEYTYPE) \
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 *); \
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; \
60 #define EIGEN_SUPERLU_HAS_ILU
63 #ifdef EIGEN_SUPERLU_HAS_ILU
66 #define DECL_GSISX(PREFIX,FLOATTYPE,KEYTYPE) \
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 *); \
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, \
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; \
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>)
95 template<
typename MatrixType>
96 struct SluMatrixMapHelper;
138 if (t==SLU_NC || t==SLU_NR || t==SLU_DN)
147 template<
typename Scalar>
150 if (internal::is_same<Scalar,float>::value)
152 else if (internal::is_same<Scalar,double>::value)
154 else if (internal::is_same<Scalar,std::complex<float> >::value)
156 else if (internal::is_same<Scalar,std::complex<double> >::value)
160 eigen_assert(
false &&
"Scalar type not supported by SuperLU");
164 template<
typename MatrixType>
167 MatrixType& mat(_mat.derived());
168 eigen_assert( ((MatrixType::Flags&
RowMajorBit)!=RowMajorBit) &&
"row-major dense matrices are not supported by SuperLU");
174 res.nrow = mat.rows();
175 res.ncol = mat.cols();
177 res.
storage.
lda = MatrixType::IsVectorAtCompileTime ? mat.size() : mat.outerStride();
182 template<
typename MatrixType>
189 res.nrow = mat.
cols();
190 res.ncol = mat.
rows();
195 res.nrow = mat.
rows();
196 res.ncol = mat.
cols();
209 if (MatrixType::Flags &
Upper)
211 if (MatrixType::Flags &
Lower)
220 template<
typename Scalar,
int Rows,
int Cols,
int Options,
int MRows,
int MCols>
221 struct SluMatrixMapHelper<
Matrix<Scalar,Rows,Cols,Options,MRows,MCols> >
226 eigen_assert( ((Options&
RowMajor)!=RowMajor) &&
"row-major dense matrices is not supported by SuperLU");
231 res.nrow = mat.
rows();
232 res.ncol = mat.
cols();
239 template<
typename Derived>
248 res.nrow = mat.cols();
249 res.ncol = mat.rows();
254 res.nrow = mat.rows();
255 res.ncol = mat.cols();
268 if (MatrixType::Flags &
Upper)
270 if (MatrixType::Flags &
Lower)
279 template<
typename MatrixType>
286 template<
typename Scalar,
int Flags,
typename Index>
290 || (Flags&
ColMajor)==ColMajor && sluMat.Stype == SLU_NC);
292 Index outerSize = (Flags&
RowMajor)==RowMajor ? sluMat.ncol : sluMat.nrow;
305 template<
typename _MatrixType,
typename Derived>
310 typedef typename MatrixType::Scalar
Scalar;
312 typedef typename MatrixType::Index
Index;
327 Derived&
derived() {
return *
static_cast<Derived*
>(
this); }
328 const Derived&
derived()
const {
return *
static_cast<const Derived*
>(
this); }
350 derived().analyzePattern(matrix);
358 template<
typename Rhs>
363 &&
"SuperLU::solve(): invalid number of rows of the right hand side matrix b");
364 return internal::solve_retval<SuperLUBase, Rhs>(*
this, b.derived());
394 template<
typename Stream>
404 const int size = a.rows();
442 Destroy_SuperNode_Matrix(&
m_sluL);
444 Destroy_CompCol_Matrix(&
m_sluU);
493 template<
typename _MatrixType>
543 #ifndef EIGEN_PARSED_BY_DOXYGEN
545 template<
typename Rhs,
typename Dest>
547 #endif // EIGEN_PARSED_BY_DOXYGEN
618 template<
typename MatrixType>
621 eigen_assert(m_analysisIsOk &&
"You must first call analyzePattern()");
628 this->initFactorization(a);
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],
640 &recip_pivot_growth, &rcond,
642 &m_sluStat, &info,
Scalar());
643 StatFree(&m_sluStat);
645 m_extractedDataAreDirty =
true;
649 m_factorizationIsOk =
true;
652 template<
typename MatrixType>
653 template<
typename Rhs,
typename Dest>
656 eigen_assert(m_factorizationIsOk &&
"The decomposition is not in a valid state for solving, you must first call either compute() or analyzePattern()/factorize()");
658 const int size = m_matrix.rows();
659 const int rhsCols = b.cols();
662 m_sluOptions.Trans = NOTRANS;
663 m_sluOptions.Fact = FACTORED;
664 m_sluOptions.IterRefine = NOREFINE;
667 m_sluFerr.resize(rhsCols);
668 m_sluBerr.resize(rhsCols);
672 typename Rhs::PlainObject b_cpy;
679 StatInit(&m_sluStat);
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],
689 &recip_pivot_growth, &rcond,
690 &m_sluFerr[0], &m_sluBerr[0],
691 &m_sluStat, &info, Scalar());
692 StatFree(&m_sluStat);
703 template<
typename MatrixType,
typename Derived>
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)
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);
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);
722 int* Lcol = m_l.outerIndexPtr();
723 int* Lrow = m_l.innerIndexPtr();
724 Scalar* Lval = m_l.valuePtr();
726 int* Ucol = m_u.outerIndexPtr();
727 int* Urow = m_u.innerIndexPtr();
728 Scalar* Uval = m_u.valuePtr();
734 for (
int k = 0; k <= Lstore->nsuper; ++k)
736 fsupc = L_FST_SUPC(k);
737 istart = L_SUB_START(fsupc);
738 nsupr = L_SUB_START(fsupc+1) - istart;
742 for (
int j = fsupc; j < L_FST_SUPC(k+1); ++j)
744 SNptr = &((
Scalar*)Lstore->nzval)[L_NZ_START(j)];
747 for (
int i = U_NZ_START(j); i < U_NZ_START(j+1); ++i)
749 Uval[lastu] = ((
Scalar*)Ustore->nzval)[i];
751 if (Uval[lastu] != 0.0)
752 Urow[lastu++] = U_SUB(i);
754 for (
int i = 0; i < upper; ++i)
757 Uval[lastu] = SNptr[i];
759 if (Uval[lastu] != 0.0)
760 Urow[lastu++] = L_SUB(istart+i);
766 Lrow[lastl++] = L_SUB(istart + upper - 1);
767 for (
int i = upper; i < nsupr; ++i)
769 Lval[lastl] = SNptr[i];
771 if (Lval[lastl] != 0.0)
772 Lrow[lastl++] = L_SUB(istart+i);
782 m_l.resizeNonZeros(lastl);
783 m_u.resizeNonZeros(lastu);
785 m_extractedDataAreDirty =
false;
789 template<
typename MatrixType>
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()");
794 if (m_extractedDataAreDirty)
798 for (
int j=0; j<m_u.cols(); ++j)
800 if (m_u.outerIndexPtr()[j+1]-m_u.outerIndexPtr()[j] > 0)
802 int lastId = m_u.outerIndexPtr()[j+1]-1;
804 if (m_u.innerIndexPtr()[lastId]==j)
805 det *= m_u.valuePtr()[lastId];
809 return det/m_sluRscale.prod()/m_sluCscale.prod();
814 #ifdef EIGEN_PARSED_BY_DOXYGEN
815 #define EIGEN_SUPERLU_HAS_ILU
818 #ifdef EIGEN_SUPERLU_HAS_ILU
834 template<
typename _MatrixType>
877 #ifndef EIGEN_PARSED_BY_DOXYGEN
879 template<
typename Rhs,
typename Dest>
881 #endif // EIGEN_PARSED_BY_DOXYGEN
933 template<
typename MatrixType>
936 eigen_assert(m_analysisIsOk &&
"You must first call analyzePattern()");
943 this->initFactorization(a);
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],
954 &recip_pivot_growth, &rcond,
955 &m_sluStat, &info,
Scalar());
956 StatFree(&m_sluStat);
960 m_factorizationIsOk =
true;
963 template<
typename MatrixType>
964 template<
typename Rhs,
typename Dest>
967 eigen_assert(m_factorizationIsOk &&
"The decomposition is not in a valid state for solving, you must first call either compute() or analyzePattern()/factorize()");
969 const int size = m_matrix.rows();
970 const int rhsCols = b.cols();
973 m_sluOptions.Trans = NOTRANS;
974 m_sluOptions.Fact = FACTORED;
975 m_sluOptions.IterRefine = NOREFINE;
977 m_sluFerr.resize(rhsCols);
978 m_sluBerr.resize(rhsCols);
982 typename Rhs::PlainObject b_cpy;
990 RealScalar recip_pivot_growth, rcond;
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],
1000 &recip_pivot_growth, &rcond,
1001 &m_sluStat, &info, Scalar());
1002 StatFree(&m_sluStat);
1008 namespace internal {
1010 template<
typename _MatrixType,
typename Derived,
typename Rhs>
1011 struct solve_retval<SuperLUBase<_MatrixType,Derived>, Rhs>
1012 : solve_retval_base<SuperLUBase<_MatrixType,Derived>, Rhs>
1014 typedef SuperLUBase<_MatrixType,Derived> Dec;
1017 template<typename Dest>
void evalTo(Dest& dst)
const
1019 dec().derived()._solve(rhs(),dst);
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>
1027 typedef SuperLUBase<_MatrixType,Derived> Dec;
1030 template<typename Dest>
void evalTo(Dest& dst)
const
1032 dec().derived()._solve(rhs(),dst);
1040 #endif // EIGEN_SUPERLUSUPPORT_H