10 #ifndef EIGEN_SUPERLUSUPPORT_H 11 #define EIGEN_SUPERLUSUPPORT_H 15 #define DECL_GSSVX(PREFIX,FLOATTYPE,KEYTYPE) \ 17 extern void PREFIX##gssvx(superlu_options_t *, SuperMatrix *, int *, int *, int *, \ 18 char *, FLOATTYPE *, FLOATTYPE *, SuperMatrix *, SuperMatrix *, \ 19 void *, int, SuperMatrix *, SuperMatrix *, \ 20 FLOATTYPE *, FLOATTYPE *, FLOATTYPE *, FLOATTYPE *, \ 21 mem_usage_t *, SuperLUStat_t *, int *); \ 23 inline float SuperLU_gssvx(superlu_options_t *options, SuperMatrix *A, \ 24 int *perm_c, int *perm_r, int *etree, char *equed, \ 25 FLOATTYPE *R, FLOATTYPE *C, SuperMatrix *L, \ 26 SuperMatrix *U, void *work, int lwork, \ 27 SuperMatrix *B, SuperMatrix *X, \ 28 FLOATTYPE *recip_pivot_growth, \ 29 FLOATTYPE *rcond, FLOATTYPE *ferr, FLOATTYPE *berr, \ 30 SuperLUStat_t *stats, int *info, KEYTYPE) { \ 31 mem_usage_t mem_usage; \ 32 PREFIX##gssvx(options, A, perm_c, perm_r, etree, equed, R, C, L, \ 33 U, work, lwork, B, X, recip_pivot_growth, rcond, \ 34 ferr, berr, &mem_usage, stats, info); \ 35 return mem_usage.for_lu; \ 38 DECL_GSSVX(s,
float,
float)
39 DECL_GSSVX(c,
float,
std::complex<
float>)
40 DECL_GSSVX(d,
double,
double)
41 DECL_GSSVX(z,
double,
std::complex<
double>)
44 #define EIGEN_SUPERLU_HAS_ILU 47 #ifdef EIGEN_SUPERLU_HAS_ILU 50 #define DECL_GSISX(PREFIX,FLOATTYPE,KEYTYPE) \ 52 extern void PREFIX##gsisx(superlu_options_t *, SuperMatrix *, int *, int *, int *, \ 53 char *, FLOATTYPE *, FLOATTYPE *, SuperMatrix *, SuperMatrix *, \ 54 void *, int, SuperMatrix *, SuperMatrix *, FLOATTYPE *, FLOATTYPE *, \ 55 mem_usage_t *, SuperLUStat_t *, int *); \ 57 inline float SuperLU_gsisx(superlu_options_t *options, SuperMatrix *A, \ 58 int *perm_c, int *perm_r, int *etree, char *equed, \ 59 FLOATTYPE *R, FLOATTYPE *C, SuperMatrix *L, \ 60 SuperMatrix *U, void *work, int lwork, \ 61 SuperMatrix *B, SuperMatrix *X, \ 62 FLOATTYPE *recip_pivot_growth, \ 64 SuperLUStat_t *stats, int *info, KEYTYPE) { \ 65 mem_usage_t mem_usage; \ 66 PREFIX##gsisx(options, A, perm_c, perm_r, etree, equed, R, C, L, \ 67 U, work, lwork, B, X, recip_pivot_growth, rcond, \ 68 &mem_usage, stats, info); \ 69 return mem_usage.for_lu; \ 72 DECL_GSISX(s,
float,
float)
73 DECL_GSISX(c,
float,
std::complex<
float>)
74 DECL_GSISX(d,
double,
double)
75 DECL_GSISX(z,
double,
std::complex<
double>)
79 template<
typename MatrixType>
80 struct SluMatrixMapHelper;
89 struct SluMatrix : SuperMatrix
96 SluMatrix(
const SluMatrix& other)
100 storage = other.storage;
103 SluMatrix& operator=(
const SluMatrix& other)
105 SuperMatrix::operator=(static_cast<const SuperMatrix&>(other));
107 storage = other.storage;
113 union {
int nnz;
int lda;};
119 void setStorageType(Stype_t t)
122 if (t==SLU_NC || t==SLU_NR || t==SLU_DN)
126 eigen_assert(
false &&
"storage type not supported");
131 template<
typename Scalar>
134 if (internal::is_same<Scalar,float>::value)
136 else if (internal::is_same<Scalar,double>::value)
138 else if (internal::is_same<Scalar,std::complex<float> >::value)
140 else if (internal::is_same<Scalar,std::complex<double> >::value)
144 eigen_assert(
false &&
"Scalar type not supported by SuperLU");
148 template<
typename MatrixType>
149 static SluMatrix Map(MatrixBase<MatrixType>& _mat)
151 MatrixType& mat(_mat.derived());
152 eigen_assert( ((MatrixType::Flags&
RowMajorBit)!=RowMajorBit) &&
"row-major dense matrices are not supported by SuperLU");
154 res.setStorageType(SLU_DN);
155 res.setScalarType<
typename MatrixType::Scalar>();
158 res.nrow = internal::convert_index<int>(mat.rows());
159 res.ncol = internal::convert_index<int>(mat.cols());
161 res.storage.lda = internal::convert_index<int>(MatrixType::IsVectorAtCompileTime ? mat.size() : mat.outerStride());
162 res.storage.values = (
void*)(mat.data());
166 template<
typename MatrixType>
167 static SluMatrix Map(SparseMatrixBase<MatrixType>& a_mat)
169 MatrixType &mat(a_mat.derived());
173 res.setStorageType(SLU_NR);
174 res.nrow = internal::convert_index<int>(mat.cols());
175 res.ncol = internal::convert_index<int>(mat.rows());
179 res.setStorageType(SLU_NC);
180 res.nrow = internal::convert_index<int>(mat.rows());
181 res.ncol = internal::convert_index<int>(mat.cols());
186 res.storage.nnz = internal::convert_index<int>(mat.nonZeros());
187 res.storage.values = mat.valuePtr();
188 res.storage.innerInd = mat.innerIndexPtr();
189 res.storage.outerInd = mat.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>
272 MappedSparseMatrix<Scalar,Flags,Index> map_superlu(SluMatrix& sluMat)
274 eigen_assert((Flags&
RowMajor)==RowMajor && sluMat.Stype == SLU_NR
275 || (Flags&
ColMajor)==ColMajor && sluMat.Stype == SLU_NC);
279 return MappedSparseMatrix<Scalar,Flags,Index>(
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>
296 using Base::m_isInitialized;
298 typedef _MatrixType MatrixType;
299 typedef typename MatrixType::Scalar Scalar;
300 typedef typename MatrixType::RealScalar RealScalar;
301 typedef typename MatrixType::StorageIndex StorageIndex;
308 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
309 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
321 inline Index rows()
const {
return m_matrix.rows(); }
322 inline Index cols()
const {
return m_matrix.cols(); }
325 inline superlu_options_t&
options() {
return m_sluOptions; }
334 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
341 derived().analyzePattern(matrix);
342 derived().factorize(matrix);
353 m_isInitialized =
true;
355 m_analysisIsOk =
true;
356 m_factorizationIsOk =
false;
359 template<
typename Stream>
360 void dumpMemory(Stream& )
365 void initFactorization(
const MatrixType& a)
367 set_default_options(&this->m_sluOptions);
369 const Index size = a.rows();
372 m_sluA = internal::asSluMatrix(m_matrix);
377 m_sluRscale.resize(size);
378 m_sluCscale.resize(size);
379 m_sluEtree.resize(size);
382 m_sluB.setStorageType(SLU_DN);
383 m_sluB.setScalarType<Scalar>();
384 m_sluB.Mtype = SLU_GE;
385 m_sluB.storage.values = 0;
388 m_sluB.storage.lda = internal::convert_index<int>(size);
391 m_extractedDataAreDirty =
true;
397 m_isInitialized =
false;
402 void extractData()
const;
407 Destroy_SuperNode_Matrix(&m_sluL);
409 Destroy_CompCol_Matrix(&m_sluU);
414 memset(&m_sluL,0,
sizeof m_sluL);
415 memset(&m_sluU,0,
sizeof m_sluU);
419 mutable LUMatrixType m_l;
420 mutable LUMatrixType m_u;
421 mutable IntColVectorType m_p;
422 mutable IntRowVectorType m_q;
424 mutable LUMatrixType m_matrix;
425 mutable SluMatrix m_sluA;
426 mutable SuperMatrix m_sluL, m_sluU;
427 mutable SluMatrix m_sluB, m_sluX;
428 mutable SuperLUStat_t m_sluStat;
429 mutable superlu_options_t m_sluOptions;
430 mutable std::vector<int> m_sluEtree;
433 mutable char m_sluEqued;
436 int m_factorizationIsOk;
438 mutable bool m_extractedDataAreDirty;
461 template<
typename _MatrixType>
466 typedef _MatrixType MatrixType;
467 typedef typename Base::Scalar Scalar;
468 typedef typename Base::RealScalar RealScalar;
469 typedef typename Base::StorageIndex StorageIndex;
478 using Base::_solve_impl;
482 explicit SuperLU(
const MatrixType& matrix) : Base()
485 Base::compute(matrix);
501 m_isInitialized =
false;
502 Base::analyzePattern(matrix);
511 void factorize(
const MatrixType& matrix);
514 template<
typename Rhs,
typename Dest>
517 inline const LMatrixType& matrixL()
const 519 if (m_extractedDataAreDirty) this->extractData();
523 inline const UMatrixType& matrixU()
const 525 if (m_extractedDataAreDirty) this->extractData();
529 inline const IntColVectorType& permutationP()
const 531 if (m_extractedDataAreDirty) this->extractData();
535 inline const IntRowVectorType& permutationQ()
const 537 if (m_extractedDataAreDirty) this->extractData();
541 Scalar determinant()
const;
545 using Base::m_matrix;
546 using Base::m_sluOptions;
552 using Base::m_sluEtree;
553 using Base::m_sluEqued;
554 using Base::m_sluRscale;
555 using Base::m_sluCscale;
558 using Base::m_sluStat;
559 using Base::m_sluFerr;
560 using Base::m_sluBerr;
564 using Base::m_analysisIsOk;
565 using Base::m_factorizationIsOk;
566 using Base::m_extractedDataAreDirty;
567 using Base::m_isInitialized;
574 set_default_options(&this->m_sluOptions);
575 m_sluOptions.PrintStat = NO;
576 m_sluOptions.ConditionNumber = NO;
577 m_sluOptions.Trans = NOTRANS;
578 m_sluOptions.ColPerm = COLAMD;
586 template<
typename MatrixType>
589 eigen_assert(m_analysisIsOk &&
"You must first call analyzePattern()");
596 this->initFactorization(a);
598 m_sluOptions.ColPerm = COLAMD;
600 RealScalar recip_pivot_growth, rcond;
601 RealScalar ferr, berr;
603 StatInit(&m_sluStat);
604 SuperLU_gssvx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0],
605 &m_sluEqued, &m_sluRscale[0], &m_sluCscale[0],
609 &recip_pivot_growth, &rcond,
611 &m_sluStat, &info, Scalar());
612 StatFree(&m_sluStat);
614 m_extractedDataAreDirty =
true;
618 m_factorizationIsOk =
true;
621 template<
typename MatrixType>
622 template<
typename Rhs,
typename Dest>
625 eigen_assert(m_factorizationIsOk &&
"The decomposition is not in a valid state for solving, you must first call either compute() or analyzePattern()/factorize()");
627 const Index size = m_matrix.rows();
628 const Index rhsCols = b.cols();
629 eigen_assert(size==b.rows());
631 m_sluOptions.Trans = NOTRANS;
632 m_sluOptions.Fact = FACTORED;
633 m_sluOptions.IterRefine = NOREFINE;
636 m_sluFerr.resize(rhsCols);
637 m_sluBerr.resize(rhsCols);
642 m_sluB = SluMatrix::Map(b_ref.const_cast_derived());
643 m_sluX = SluMatrix::Map(x_ref.const_cast_derived());
645 typename Rhs::PlainObject b_cpy;
649 m_sluB = SluMatrix::Map(b_cpy.const_cast_derived());
652 StatInit(&m_sluStat);
654 RealScalar recip_pivot_growth, rcond;
655 SuperLU_gssvx(&m_sluOptions, &m_sluA,
656 m_q.data(), m_p.data(),
657 &m_sluEtree[0], &m_sluEqued,
658 &m_sluRscale[0], &m_sluCscale[0],
662 &recip_pivot_growth, &rcond,
663 &m_sluFerr[0], &m_sluBerr[0],
664 &m_sluStat, &info, Scalar());
665 StatFree(&m_sluStat);
667 if(x.derived().data() != x_ref.data())
680 template<
typename MatrixType,
typename Derived>
683 eigen_assert(m_factorizationIsOk &&
"The decomposition is not in a valid state for extracting factors, you must first call either compute() or analyzePattern()/factorize()");
684 if (m_extractedDataAreDirty)
687 int fsupc, istart, nsupr;
688 int lastl = 0, lastu = 0;
689 SCformat *Lstore =
static_cast<SCformat*
>(m_sluL.Store);
690 NCformat *Ustore =
static_cast<NCformat*
>(m_sluU.Store);
693 const Index size = m_matrix.rows();
694 m_l.resize(size,size);
695 m_l.resizeNonZeros(Lstore->nnz);
696 m_u.resize(size,size);
697 m_u.resizeNonZeros(Ustore->nnz);
699 int* Lcol = m_l.outerIndexPtr();
700 int* Lrow = m_l.innerIndexPtr();
701 Scalar* Lval = m_l.valuePtr();
703 int* Ucol = m_u.outerIndexPtr();
704 int* Urow = m_u.innerIndexPtr();
705 Scalar* Uval = m_u.valuePtr();
711 for (
int k = 0; k <= Lstore->nsuper; ++k)
713 fsupc = L_FST_SUPC(k);
714 istart = L_SUB_START(fsupc);
715 nsupr = L_SUB_START(fsupc+1) - istart;
719 for (
int j = fsupc; j < L_FST_SUPC(k+1); ++j)
721 SNptr = &((Scalar*)Lstore->nzval)[L_NZ_START(j)];
724 for (
int i = U_NZ_START(j); i < U_NZ_START(j+1); ++i)
726 Uval[lastu] = ((Scalar*)Ustore->nzval)[i];
728 if (Uval[lastu] != 0.0)
729 Urow[lastu++] = U_SUB(i);
731 for (
int i = 0; i < upper; ++i)
734 Uval[lastu] = SNptr[i];
736 if (Uval[lastu] != 0.0)
737 Urow[lastu++] = L_SUB(istart+i);
743 Lrow[lastl++] = L_SUB(istart + upper - 1);
744 for (
int i = upper; i < nsupr; ++i)
746 Lval[lastl] = SNptr[i];
748 if (Lval[lastl] != 0.0)
749 Lrow[lastl++] = L_SUB(istart+i);
759 m_l.resizeNonZeros(lastl);
760 m_u.resizeNonZeros(lastu);
762 m_extractedDataAreDirty =
false;
766 template<
typename MatrixType>
769 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()");
771 if (m_extractedDataAreDirty)
774 Scalar det = Scalar(1);
775 for (
int j=0; j<m_u.cols(); ++j)
777 if (m_u.outerIndexPtr()[j+1]-m_u.outerIndexPtr()[j] > 0)
779 int lastId = m_u.outerIndexPtr()[j+1]-1;
780 eigen_assert(m_u.innerIndexPtr()[lastId]<=j);
781 if (m_u.innerIndexPtr()[lastId]==j)
782 det *= m_u.valuePtr()[lastId];
788 return det/m_sluRscale.prod()/m_sluCscale.prod();
793 #ifdef EIGEN_PARSED_BY_DOXYGEN 794 #define EIGEN_SUPERLU_HAS_ILU 797 #ifdef EIGEN_SUPERLU_HAS_ILU 815 template<
typename _MatrixType>
820 typedef _MatrixType MatrixType;
821 typedef typename Base::Scalar Scalar;
822 typedef typename Base::RealScalar RealScalar;
825 using Base::_solve_impl;
829 SuperILU(
const MatrixType& matrix) : Base()
832 Base::compute(matrix);
847 Base::analyzePattern(matrix);
856 void factorize(
const MatrixType& matrix);
858 #ifndef EIGEN_PARSED_BY_DOXYGEN 860 template<
typename Rhs,
typename Dest>
862 #endif // EIGEN_PARSED_BY_DOXYGEN 866 using Base::m_matrix;
867 using Base::m_sluOptions;
873 using Base::m_sluEtree;
874 using Base::m_sluEqued;
875 using Base::m_sluRscale;
876 using Base::m_sluCscale;
879 using Base::m_sluStat;
880 using Base::m_sluFerr;
881 using Base::m_sluBerr;
885 using Base::m_analysisIsOk;
886 using Base::m_factorizationIsOk;
887 using Base::m_extractedDataAreDirty;
888 using Base::m_isInitialized;
895 ilu_set_default_options(&m_sluOptions);
896 m_sluOptions.PrintStat = NO;
897 m_sluOptions.ConditionNumber = NO;
898 m_sluOptions.Trans = NOTRANS;
899 m_sluOptions.ColPerm = MMD_AT_PLUS_A;
902 m_sluOptions.ILU_MILU = SILU;
906 m_sluOptions.ILU_DropRule = DROP_BASIC;
914 template<
typename MatrixType>
917 eigen_assert(m_analysisIsOk &&
"You must first call analyzePattern()");
924 this->initFactorization(a);
927 RealScalar recip_pivot_growth, rcond;
929 StatInit(&m_sluStat);
930 SuperLU_gsisx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0],
931 &m_sluEqued, &m_sluRscale[0], &m_sluCscale[0],
935 &recip_pivot_growth, &rcond,
936 &m_sluStat, &info, Scalar());
937 StatFree(&m_sluStat);
941 m_factorizationIsOk =
true;
944 template<
typename MatrixType>
945 template<
typename Rhs,
typename Dest>
948 eigen_assert(m_factorizationIsOk &&
"The decomposition is not in a valid state for solving, you must first call either compute() or analyzePattern()/factorize()");
950 const int size = m_matrix.rows();
951 const int rhsCols = b.cols();
952 eigen_assert(size==b.rows());
954 m_sluOptions.Trans = NOTRANS;
955 m_sluOptions.Fact = FACTORED;
956 m_sluOptions.IterRefine = NOREFINE;
958 m_sluFerr.resize(rhsCols);
959 m_sluBerr.resize(rhsCols);
964 m_sluB = SluMatrix::Map(b_ref.const_cast_derived());
965 m_sluX = SluMatrix::Map(x_ref.const_cast_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);
989 if(&x.coeffRef(0) != x_ref.data())
998 #endif // EIGEN_SUPERLUSUPPORT_H void compute(const MatrixType &matrix)
Definition: SuperLUSupport.h:339
void analyzePattern(const MatrixType &matrix)
Definition: SuperLUSupport.h:845
A sparse direct LU factorization and solver based on the SuperLU library.
Definition: SuperLUSupport.h:462
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:816
A base class for sparse solvers.
Definition: SparseSolverBase.h:53
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:332
const unsigned int RowMajorBit
Definition: Constants.h:61
void analyzePattern(const MatrixType &)
Definition: SuperLUSupport.h:351
Definition: Constants.h:320
Definition: Constants.h:204
Definition: Constants.h:322
void factorize(const MatrixType &matrix)
Definition: SuperLUSupport.h:915
Definition: Constants.h:434
The base class for the direct and incomplete LU factorization of SuperLU.
Definition: SuperLUSupport.h:291
Definition: Constants.h:439
Definition: Constants.h:432
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
void factorize(const MatrixType &matrix)
Definition: SuperLUSupport.h:587
Definition: Constants.h:220
superlu_options_t & options()
Definition: SuperLUSupport.h:325
Definition: Constants.h:206
ComputationInfo
Definition: Constants.h:430
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
void analyzePattern(const MatrixType &matrix)
Definition: SuperLUSupport.h:498