00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025 #ifndef EIGEN_SUPERLUSUPPORT_H
00026 #define EIGEN_SUPERLUSUPPORT_H
00027
00028 namespace Eigen {
00029
00030 #define DECL_GSSVX(PREFIX,FLOATTYPE,KEYTYPE) \
00031 extern "C" { \
00032 typedef struct { FLOATTYPE for_lu; FLOATTYPE total_needed; int expansions; } PREFIX##mem_usage_t; \
00033 extern void PREFIX##gssvx(superlu_options_t *, SuperMatrix *, int *, int *, int *, \
00034 char *, FLOATTYPE *, FLOATTYPE *, SuperMatrix *, SuperMatrix *, \
00035 void *, int, SuperMatrix *, SuperMatrix *, \
00036 FLOATTYPE *, FLOATTYPE *, FLOATTYPE *, FLOATTYPE *, \
00037 PREFIX##mem_usage_t *, SuperLUStat_t *, int *); \
00038 } \
00039 inline float SuperLU_gssvx(superlu_options_t *options, SuperMatrix *A, \
00040 int *perm_c, int *perm_r, int *etree, char *equed, \
00041 FLOATTYPE *R, FLOATTYPE *C, SuperMatrix *L, \
00042 SuperMatrix *U, void *work, int lwork, \
00043 SuperMatrix *B, SuperMatrix *X, \
00044 FLOATTYPE *recip_pivot_growth, \
00045 FLOATTYPE *rcond, FLOATTYPE *ferr, FLOATTYPE *berr, \
00046 SuperLUStat_t *stats, int *info, KEYTYPE) { \
00047 PREFIX##mem_usage_t mem_usage; \
00048 PREFIX##gssvx(options, A, perm_c, perm_r, etree, equed, R, C, L, \
00049 U, work, lwork, B, X, recip_pivot_growth, rcond, \
00050 ferr, berr, &mem_usage, stats, info); \
00051 return mem_usage.for_lu; \
00052 }
00053
00054 DECL_GSSVX(s,float,float)
00055 DECL_GSSVX(c,float,std::complex<float>)
00056 DECL_GSSVX(d,double,double)
00057 DECL_GSSVX(z,double,std::complex<double>)
00058
00059 #ifdef MILU_ALPHA
00060 #define EIGEN_SUPERLU_HAS_ILU
00061 #endif
00062
00063 #ifdef EIGEN_SUPERLU_HAS_ILU
00064
00065
00066 #define DECL_GSISX(PREFIX,FLOATTYPE,KEYTYPE) \
00067 extern "C" { \
00068 extern void PREFIX##gsisx(superlu_options_t *, SuperMatrix *, int *, int *, int *, \
00069 char *, FLOATTYPE *, FLOATTYPE *, SuperMatrix *, SuperMatrix *, \
00070 void *, int, SuperMatrix *, SuperMatrix *, FLOATTYPE *, FLOATTYPE *, \
00071 PREFIX##mem_usage_t *, SuperLUStat_t *, int *); \
00072 } \
00073 inline float SuperLU_gsisx(superlu_options_t *options, SuperMatrix *A, \
00074 int *perm_c, int *perm_r, int *etree, char *equed, \
00075 FLOATTYPE *R, FLOATTYPE *C, SuperMatrix *L, \
00076 SuperMatrix *U, void *work, int lwork, \
00077 SuperMatrix *B, SuperMatrix *X, \
00078 FLOATTYPE *recip_pivot_growth, \
00079 FLOATTYPE *rcond, \
00080 SuperLUStat_t *stats, int *info, KEYTYPE) { \
00081 PREFIX##mem_usage_t mem_usage; \
00082 PREFIX##gsisx(options, A, perm_c, perm_r, etree, equed, R, C, L, \
00083 U, work, lwork, B, X, recip_pivot_growth, rcond, \
00084 &mem_usage, stats, info); \
00085 return mem_usage.for_lu; \
00086 }
00087
00088 DECL_GSISX(s,float,float)
00089 DECL_GSISX(c,float,std::complex<float>)
00090 DECL_GSISX(d,double,double)
00091 DECL_GSISX(z,double,std::complex<double>)
00092
00093 #endif
00094
00095 template<typename MatrixType>
00096 struct SluMatrixMapHelper;
00097
00105 struct SluMatrix : SuperMatrix
00106 {
00107 SluMatrix()
00108 {
00109 Store = &storage;
00110 }
00111
00112 SluMatrix(const SluMatrix& other)
00113 : SuperMatrix(other)
00114 {
00115 Store = &storage;
00116 storage = other.storage;
00117 }
00118
00119 SluMatrix& operator=(const SluMatrix& other)
00120 {
00121 SuperMatrix::operator=(static_cast<const SuperMatrix&>(other));
00122 Store = &storage;
00123 storage = other.storage;
00124 return *this;
00125 }
00126
00127 struct
00128 {
00129 union {int nnz;int lda;};
00130 void *values;
00131 int *innerInd;
00132 int *outerInd;
00133 } storage;
00134
00135 void setStorageType(Stype_t t)
00136 {
00137 Stype = t;
00138 if (t==SLU_NC || t==SLU_NR || t==SLU_DN)
00139 Store = &storage;
00140 else
00141 {
00142 eigen_assert(false && "storage type not supported");
00143 Store = 0;
00144 }
00145 }
00146
00147 template<typename Scalar>
00148 void setScalarType()
00149 {
00150 if (internal::is_same<Scalar,float>::value)
00151 Dtype = SLU_S;
00152 else if (internal::is_same<Scalar,double>::value)
00153 Dtype = SLU_D;
00154 else if (internal::is_same<Scalar,std::complex<float> >::value)
00155 Dtype = SLU_C;
00156 else if (internal::is_same<Scalar,std::complex<double> >::value)
00157 Dtype = SLU_Z;
00158 else
00159 {
00160 eigen_assert(false && "Scalar type not supported by SuperLU");
00161 }
00162 }
00163
00164 template<typename MatrixType>
00165 static SluMatrix Map(MatrixBase<MatrixType>& _mat)
00166 {
00167 MatrixType& mat(_mat.derived());
00168 eigen_assert( ((MatrixType::Flags&RowMajorBit)!=RowMajorBit) && "row-major dense matrices are not supported by SuperLU");
00169 SluMatrix res;
00170 res.setStorageType(SLU_DN);
00171 res.setScalarType<typename MatrixType::Scalar>();
00172 res.Mtype = SLU_GE;
00173
00174 res.nrow = mat.rows();
00175 res.ncol = mat.cols();
00176
00177 res.storage.lda = MatrixType::IsVectorAtCompileTime ? mat.size() : mat.outerStride();
00178 res.storage.values = mat.data();
00179 return res;
00180 }
00181
00182 template<typename MatrixType>
00183 static SluMatrix Map(SparseMatrixBase<MatrixType>& mat)
00184 {
00185 SluMatrix res;
00186 if ((MatrixType::Flags&RowMajorBit)==RowMajorBit)
00187 {
00188 res.setStorageType(SLU_NR);
00189 res.nrow = mat.cols();
00190 res.ncol = mat.rows();
00191 }
00192 else
00193 {
00194 res.setStorageType(SLU_NC);
00195 res.nrow = mat.rows();
00196 res.ncol = mat.cols();
00197 }
00198
00199 res.Mtype = SLU_GE;
00200
00201 res.storage.nnz = mat.nonZeros();
00202 res.storage.values = mat.derived().valuePtr();
00203 res.storage.innerInd = mat.derived().innerIndexPtr();
00204 res.storage.outerInd = mat.derived().outerIndexPtr();
00205
00206 res.setScalarType<typename MatrixType::Scalar>();
00207
00208
00209 if (MatrixType::Flags & Upper)
00210 res.Mtype = SLU_TRU;
00211 if (MatrixType::Flags & Lower)
00212 res.Mtype = SLU_TRL;
00213
00214 eigen_assert(((MatrixType::Flags & SelfAdjoint)==0) && "SelfAdjoint matrix shape not supported by SuperLU");
00215
00216 return res;
00217 }
00218 };
00219
00220 template<typename Scalar, int Rows, int Cols, int Options, int MRows, int MCols>
00221 struct SluMatrixMapHelper<Matrix<Scalar,Rows,Cols,Options,MRows,MCols> >
00222 {
00223 typedef Matrix<Scalar,Rows,Cols,Options,MRows,MCols> MatrixType;
00224 static void run(MatrixType& mat, SluMatrix& res)
00225 {
00226 eigen_assert( ((Options&RowMajor)!=RowMajor) && "row-major dense matrices is not supported by SuperLU");
00227 res.setStorageType(SLU_DN);
00228 res.setScalarType<Scalar>();
00229 res.Mtype = SLU_GE;
00230
00231 res.nrow = mat.rows();
00232 res.ncol = mat.cols();
00233
00234 res.storage.lda = mat.outerStride();
00235 res.storage.values = mat.data();
00236 }
00237 };
00238
00239 template<typename Derived>
00240 struct SluMatrixMapHelper<SparseMatrixBase<Derived> >
00241 {
00242 typedef Derived MatrixType;
00243 static void run(MatrixType& mat, SluMatrix& res)
00244 {
00245 if ((MatrixType::Flags&RowMajorBit)==RowMajorBit)
00246 {
00247 res.setStorageType(SLU_NR);
00248 res.nrow = mat.cols();
00249 res.ncol = mat.rows();
00250 }
00251 else
00252 {
00253 res.setStorageType(SLU_NC);
00254 res.nrow = mat.rows();
00255 res.ncol = mat.cols();
00256 }
00257
00258 res.Mtype = SLU_GE;
00259
00260 res.storage.nnz = mat.nonZeros();
00261 res.storage.values = mat.valuePtr();
00262 res.storage.innerInd = mat.innerIndexPtr();
00263 res.storage.outerInd = mat.outerIndexPtr();
00264
00265 res.setScalarType<typename MatrixType::Scalar>();
00266
00267
00268 if (MatrixType::Flags & Upper)
00269 res.Mtype = SLU_TRU;
00270 if (MatrixType::Flags & Lower)
00271 res.Mtype = SLU_TRL;
00272
00273 eigen_assert(((MatrixType::Flags & SelfAdjoint)==0) && "SelfAdjoint matrix shape not supported by SuperLU");
00274 }
00275 };
00276
00277 namespace internal {
00278
00279 template<typename MatrixType>
00280 SluMatrix asSluMatrix(MatrixType& mat)
00281 {
00282 return SluMatrix::Map(mat);
00283 }
00284
00286 template<typename Scalar, int Flags, typename Index>
00287 MappedSparseMatrix<Scalar,Flags,Index> map_superlu(SluMatrix& sluMat)
00288 {
00289 eigen_assert((Flags&RowMajor)==RowMajor && sluMat.Stype == SLU_NR
00290 || (Flags&ColMajor)==ColMajor && sluMat.Stype == SLU_NC);
00291
00292 Index outerSize = (Flags&RowMajor)==RowMajor ? sluMat.ncol : sluMat.nrow;
00293
00294 return MappedSparseMatrix<Scalar,Flags,Index>(
00295 sluMat.nrow, sluMat.ncol, sluMat.storage.outerInd[outerSize],
00296 sluMat.storage.outerInd, sluMat.storage.innerInd, reinterpret_cast<Scalar*>(sluMat.storage.values) );
00297 }
00298
00299 }
00300
00305 template<typename _MatrixType, typename Derived>
00306 class SuperLUBase : internal::noncopyable
00307 {
00308 public:
00309 typedef _MatrixType MatrixType;
00310 typedef typename MatrixType::Scalar Scalar;
00311 typedef typename MatrixType::RealScalar RealScalar;
00312 typedef typename MatrixType::Index Index;
00313 typedef Matrix<Scalar,Dynamic,1> Vector;
00314 typedef Matrix<int, 1, MatrixType::ColsAtCompileTime> IntRowVectorType;
00315 typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType;
00316 typedef SparseMatrix<Scalar> LUMatrixType;
00317
00318 public:
00319
00320 SuperLUBase() {}
00321
00322 ~SuperLUBase()
00323 {
00324 clearFactors();
00325 }
00326
00327 Derived& derived() { return *static_cast<Derived*>(this); }
00328 const Derived& derived() const { return *static_cast<const Derived*>(this); }
00329
00330 inline Index rows() const { return m_matrix.rows(); }
00331 inline Index cols() const { return m_matrix.cols(); }
00332
00334 inline superlu_options_t& options() { return m_sluOptions; }
00335
00341 ComputationInfo info() const
00342 {
00343 eigen_assert(m_isInitialized && "Decomposition is not initialized.");
00344 return m_info;
00345 }
00346
00348 void compute(const MatrixType& matrix)
00349 {
00350 derived().analyzePattern(matrix);
00351 derived().factorize(matrix);
00352 }
00353
00358 template<typename Rhs>
00359 inline const internal::solve_retval<SuperLUBase, Rhs> solve(const MatrixBase<Rhs>& b) const
00360 {
00361 eigen_assert(m_isInitialized && "SuperLU is not initialized.");
00362 eigen_assert(rows()==b.rows()
00363 && "SuperLU::solve(): invalid number of rows of the right hand side matrix b");
00364 return internal::solve_retval<SuperLUBase, Rhs>(*this, b.derived());
00365 }
00366
00371
00372
00373
00374
00375
00376
00377
00378
00379
00386 void analyzePattern(const MatrixType& )
00387 {
00388 m_isInitialized = true;
00389 m_info = Success;
00390 m_analysisIsOk = true;
00391 m_factorizationIsOk = false;
00392 }
00393
00394 template<typename Stream>
00395 void dumpMemory(Stream& s)
00396 {}
00397
00398 protected:
00399
00400 void initFactorization(const MatrixType& a)
00401 {
00402 set_default_options(&this->m_sluOptions);
00403
00404 const int size = a.rows();
00405 m_matrix = a;
00406
00407 m_sluA = internal::asSluMatrix(m_matrix);
00408 clearFactors();
00409
00410 m_p.resize(size);
00411 m_q.resize(size);
00412 m_sluRscale.resize(size);
00413 m_sluCscale.resize(size);
00414 m_sluEtree.resize(size);
00415
00416
00417 m_sluB.setStorageType(SLU_DN);
00418 m_sluB.setScalarType<Scalar>();
00419 m_sluB.Mtype = SLU_GE;
00420 m_sluB.storage.values = 0;
00421 m_sluB.nrow = 0;
00422 m_sluB.ncol = 0;
00423 m_sluB.storage.lda = size;
00424 m_sluX = m_sluB;
00425
00426 m_extractedDataAreDirty = true;
00427 }
00428
00429 void init()
00430 {
00431 m_info = InvalidInput;
00432 m_isInitialized = false;
00433 m_sluL.Store = 0;
00434 m_sluU.Store = 0;
00435 }
00436
00437 void extractData() const;
00438
00439 void clearFactors()
00440 {
00441 if(m_sluL.Store)
00442 Destroy_SuperNode_Matrix(&m_sluL);
00443 if(m_sluU.Store)
00444 Destroy_CompCol_Matrix(&m_sluU);
00445
00446 m_sluL.Store = 0;
00447 m_sluU.Store = 0;
00448
00449 memset(&m_sluL,0,sizeof m_sluL);
00450 memset(&m_sluU,0,sizeof m_sluU);
00451 }
00452
00453
00454 mutable LUMatrixType m_l;
00455 mutable LUMatrixType m_u;
00456 mutable IntColVectorType m_p;
00457 mutable IntRowVectorType m_q;
00458
00459 mutable LUMatrixType m_matrix;
00460 mutable SluMatrix m_sluA;
00461 mutable SuperMatrix m_sluL, m_sluU;
00462 mutable SluMatrix m_sluB, m_sluX;
00463 mutable SuperLUStat_t m_sluStat;
00464 mutable superlu_options_t m_sluOptions;
00465 mutable std::vector<int> m_sluEtree;
00466 mutable Matrix<RealScalar,Dynamic,1> m_sluRscale, m_sluCscale;
00467 mutable Matrix<RealScalar,Dynamic,1> m_sluFerr, m_sluBerr;
00468 mutable char m_sluEqued;
00469
00470 mutable ComputationInfo m_info;
00471 bool m_isInitialized;
00472 int m_factorizationIsOk;
00473 int m_analysisIsOk;
00474 mutable bool m_extractedDataAreDirty;
00475
00476 private:
00477 SuperLUBase(SuperLUBase& ) { }
00478 };
00479
00480
00493 template<typename _MatrixType>
00494 class SuperLU : public SuperLUBase<_MatrixType,SuperLU<_MatrixType> >
00495 {
00496 public:
00497 typedef SuperLUBase<_MatrixType,SuperLU> Base;
00498 typedef _MatrixType MatrixType;
00499 typedef typename Base::Scalar Scalar;
00500 typedef typename Base::RealScalar RealScalar;
00501 typedef typename Base::Index Index;
00502 typedef typename Base::IntRowVectorType IntRowVectorType;
00503 typedef typename Base::IntColVectorType IntColVectorType;
00504 typedef typename Base::LUMatrixType LUMatrixType;
00505 typedef TriangularView<LUMatrixType, Lower|UnitDiag> LMatrixType;
00506 typedef TriangularView<LUMatrixType, Upper> UMatrixType;
00507
00508 public:
00509
00510 SuperLU() : Base() { init(); }
00511
00512 SuperLU(const MatrixType& matrix) : Base()
00513 {
00514 Base::init();
00515 compute(matrix);
00516 }
00517
00518 ~SuperLU()
00519 {
00520 }
00521
00528 void analyzePattern(const MatrixType& matrix)
00529 {
00530 m_info = InvalidInput;
00531 m_isInitialized = false;
00532 Base::analyzePattern(matrix);
00533 }
00534
00541 void factorize(const MatrixType& matrix);
00542
00543 #ifndef EIGEN_PARSED_BY_DOXYGEN
00544
00545 template<typename Rhs,typename Dest>
00546 void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const;
00547 #endif // EIGEN_PARSED_BY_DOXYGEN
00548
00549 inline const LMatrixType& matrixL() const
00550 {
00551 if (m_extractedDataAreDirty) this->extractData();
00552 return m_l;
00553 }
00554
00555 inline const UMatrixType& matrixU() const
00556 {
00557 if (m_extractedDataAreDirty) this->extractData();
00558 return m_u;
00559 }
00560
00561 inline const IntColVectorType& permutationP() const
00562 {
00563 if (m_extractedDataAreDirty) this->extractData();
00564 return m_p;
00565 }
00566
00567 inline const IntRowVectorType& permutationQ() const
00568 {
00569 if (m_extractedDataAreDirty) this->extractData();
00570 return m_q;
00571 }
00572
00573 Scalar determinant() const;
00574
00575 protected:
00576
00577 using Base::m_matrix;
00578 using Base::m_sluOptions;
00579 using Base::m_sluA;
00580 using Base::m_sluB;
00581 using Base::m_sluX;
00582 using Base::m_p;
00583 using Base::m_q;
00584 using Base::m_sluEtree;
00585 using Base::m_sluEqued;
00586 using Base::m_sluRscale;
00587 using Base::m_sluCscale;
00588 using Base::m_sluL;
00589 using Base::m_sluU;
00590 using Base::m_sluStat;
00591 using Base::m_sluFerr;
00592 using Base::m_sluBerr;
00593 using Base::m_l;
00594 using Base::m_u;
00595
00596 using Base::m_analysisIsOk;
00597 using Base::m_factorizationIsOk;
00598 using Base::m_extractedDataAreDirty;
00599 using Base::m_isInitialized;
00600 using Base::m_info;
00601
00602 void init()
00603 {
00604 Base::init();
00605
00606 set_default_options(&this->m_sluOptions);
00607 m_sluOptions.PrintStat = NO;
00608 m_sluOptions.ConditionNumber = NO;
00609 m_sluOptions.Trans = NOTRANS;
00610 m_sluOptions.ColPerm = COLAMD;
00611 }
00612
00613
00614 private:
00615 SuperLU(SuperLU& ) { }
00616 };
00617
00618 template<typename MatrixType>
00619 void SuperLU<MatrixType>::factorize(const MatrixType& a)
00620 {
00621 eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
00622 if(!m_analysisIsOk)
00623 {
00624 m_info = InvalidInput;
00625 return;
00626 }
00627
00628 this->initFactorization(a);
00629
00630 int info = 0;
00631 RealScalar recip_pivot_growth, rcond;
00632 RealScalar ferr, berr;
00633
00634 StatInit(&m_sluStat);
00635 SuperLU_gssvx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0],
00636 &m_sluEqued, &m_sluRscale[0], &m_sluCscale[0],
00637 &m_sluL, &m_sluU,
00638 NULL, 0,
00639 &m_sluB, &m_sluX,
00640 &recip_pivot_growth, &rcond,
00641 &ferr, &berr,
00642 &m_sluStat, &info, Scalar());
00643 StatFree(&m_sluStat);
00644
00645 m_extractedDataAreDirty = true;
00646
00647
00648 m_info = info == 0 ? Success : NumericalIssue;
00649 m_factorizationIsOk = true;
00650 }
00651
00652 template<typename MatrixType>
00653 template<typename Rhs,typename Dest>
00654 void SuperLU<MatrixType>::_solve(const MatrixBase<Rhs> &b, MatrixBase<Dest>& x) const
00655 {
00656 eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or analyzePattern()/factorize()");
00657
00658 const int size = m_matrix.rows();
00659 const int rhsCols = b.cols();
00660 eigen_assert(size==b.rows());
00661
00662 m_sluOptions.Trans = NOTRANS;
00663 m_sluOptions.Fact = FACTORED;
00664 m_sluOptions.IterRefine = NOREFINE;
00665
00666
00667 m_sluFerr.resize(rhsCols);
00668 m_sluBerr.resize(rhsCols);
00669 m_sluB = SluMatrix::Map(b.const_cast_derived());
00670 m_sluX = SluMatrix::Map(x.derived());
00671
00672 typename Rhs::PlainObject b_cpy;
00673 if(m_sluEqued!='N')
00674 {
00675 b_cpy = b;
00676 m_sluB = SluMatrix::Map(b_cpy.const_cast_derived());
00677 }
00678
00679 StatInit(&m_sluStat);
00680 int info = 0;
00681 RealScalar recip_pivot_growth, rcond;
00682 SuperLU_gssvx(&m_sluOptions, &m_sluA,
00683 m_q.data(), m_p.data(),
00684 &m_sluEtree[0], &m_sluEqued,
00685 &m_sluRscale[0], &m_sluCscale[0],
00686 &m_sluL, &m_sluU,
00687 NULL, 0,
00688 &m_sluB, &m_sluX,
00689 &recip_pivot_growth, &rcond,
00690 &m_sluFerr[0], &m_sluBerr[0],
00691 &m_sluStat, &info, Scalar());
00692 StatFree(&m_sluStat);
00693 m_info = info==0 ? Success : NumericalIssue;
00694 }
00695
00696
00697
00698
00699
00700
00701
00702
00703 template<typename MatrixType, typename Derived>
00704 void SuperLUBase<MatrixType,Derived>::extractData() const
00705 {
00706 eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for extracting factors, you must first call either compute() or analyzePattern()/factorize()");
00707 if (m_extractedDataAreDirty)
00708 {
00709 int upper;
00710 int fsupc, istart, nsupr;
00711 int lastl = 0, lastu = 0;
00712 SCformat *Lstore = static_cast<SCformat*>(m_sluL.Store);
00713 NCformat *Ustore = static_cast<NCformat*>(m_sluU.Store);
00714 Scalar *SNptr;
00715
00716 const int size = m_matrix.rows();
00717 m_l.resize(size,size);
00718 m_l.resizeNonZeros(Lstore->nnz);
00719 m_u.resize(size,size);
00720 m_u.resizeNonZeros(Ustore->nnz);
00721
00722 int* Lcol = m_l.outerIndexPtr();
00723 int* Lrow = m_l.innerIndexPtr();
00724 Scalar* Lval = m_l.valuePtr();
00725
00726 int* Ucol = m_u.outerIndexPtr();
00727 int* Urow = m_u.innerIndexPtr();
00728 Scalar* Uval = m_u.valuePtr();
00729
00730 Ucol[0] = 0;
00731 Ucol[0] = 0;
00732
00733
00734 for (int k = 0; k <= Lstore->nsuper; ++k)
00735 {
00736 fsupc = L_FST_SUPC(k);
00737 istart = L_SUB_START(fsupc);
00738 nsupr = L_SUB_START(fsupc+1) - istart;
00739 upper = 1;
00740
00741
00742 for (int j = fsupc; j < L_FST_SUPC(k+1); ++j)
00743 {
00744 SNptr = &((Scalar*)Lstore->nzval)[L_NZ_START(j)];
00745
00746
00747 for (int i = U_NZ_START(j); i < U_NZ_START(j+1); ++i)
00748 {
00749 Uval[lastu] = ((Scalar*)Ustore->nzval)[i];
00750
00751 if (Uval[lastu] != 0.0)
00752 Urow[lastu++] = U_SUB(i);
00753 }
00754 for (int i = 0; i < upper; ++i)
00755 {
00756
00757 Uval[lastu] = SNptr[i];
00758
00759 if (Uval[lastu] != 0.0)
00760 Urow[lastu++] = L_SUB(istart+i);
00761 }
00762 Ucol[j+1] = lastu;
00763
00764
00765 Lval[lastl] = 1.0;
00766 Lrow[lastl++] = L_SUB(istart + upper - 1);
00767 for (int i = upper; i < nsupr; ++i)
00768 {
00769 Lval[lastl] = SNptr[i];
00770
00771 if (Lval[lastl] != 0.0)
00772 Lrow[lastl++] = L_SUB(istart+i);
00773 }
00774 Lcol[j+1] = lastl;
00775
00776 ++upper;
00777 }
00778
00779 }
00780
00781
00782 m_l.resizeNonZeros(lastl);
00783 m_u.resizeNonZeros(lastu);
00784
00785 m_extractedDataAreDirty = false;
00786 }
00787 }
00788
00789 template<typename MatrixType>
00790 typename SuperLU<MatrixType>::Scalar SuperLU<MatrixType>::determinant() const
00791 {
00792 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()");
00793
00794 if (m_extractedDataAreDirty)
00795 this->extractData();
00796
00797 Scalar det = Scalar(1);
00798 for (int j=0; j<m_u.cols(); ++j)
00799 {
00800 if (m_u.outerIndexPtr()[j+1]-m_u.outerIndexPtr()[j] > 0)
00801 {
00802 int lastId = m_u.outerIndexPtr()[j+1]-1;
00803 eigen_assert(m_u.innerIndexPtr()[lastId]<=j);
00804 if (m_u.innerIndexPtr()[lastId]==j)
00805 det *= m_u.valuePtr()[lastId];
00806 }
00807 }
00808 if(m_sluEqued!='N')
00809 return det/m_sluRscale.prod()/m_sluCscale.prod();
00810 else
00811 return det;
00812 }
00813
00814 #ifdef EIGEN_PARSED_BY_DOXYGEN
00815 #define EIGEN_SUPERLU_HAS_ILU
00816 #endif
00817
00818 #ifdef EIGEN_SUPERLU_HAS_ILU
00819
00834 template<typename _MatrixType>
00835 class SuperILU : public SuperLUBase<_MatrixType,SuperILU<_MatrixType> >
00836 {
00837 public:
00838 typedef SuperLUBase<_MatrixType,SuperILU> Base;
00839 typedef _MatrixType MatrixType;
00840 typedef typename Base::Scalar Scalar;
00841 typedef typename Base::RealScalar RealScalar;
00842 typedef typename Base::Index Index;
00843
00844 public:
00845
00846 SuperILU() : Base() { init(); }
00847
00848 SuperILU(const MatrixType& matrix) : Base()
00849 {
00850 init();
00851 compute(matrix);
00852 }
00853
00854 ~SuperILU()
00855 {
00856 }
00857
00864 void analyzePattern(const MatrixType& matrix)
00865 {
00866 Base::analyzePattern(matrix);
00867 }
00868
00875 void factorize(const MatrixType& matrix);
00876
00877 #ifndef EIGEN_PARSED_BY_DOXYGEN
00878
00879 template<typename Rhs,typename Dest>
00880 void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const;
00881 #endif // EIGEN_PARSED_BY_DOXYGEN
00882
00883 protected:
00884
00885 using Base::m_matrix;
00886 using Base::m_sluOptions;
00887 using Base::m_sluA;
00888 using Base::m_sluB;
00889 using Base::m_sluX;
00890 using Base::m_p;
00891 using Base::m_q;
00892 using Base::m_sluEtree;
00893 using Base::m_sluEqued;
00894 using Base::m_sluRscale;
00895 using Base::m_sluCscale;
00896 using Base::m_sluL;
00897 using Base::m_sluU;
00898 using Base::m_sluStat;
00899 using Base::m_sluFerr;
00900 using Base::m_sluBerr;
00901 using Base::m_l;
00902 using Base::m_u;
00903
00904 using Base::m_analysisIsOk;
00905 using Base::m_factorizationIsOk;
00906 using Base::m_extractedDataAreDirty;
00907 using Base::m_isInitialized;
00908 using Base::m_info;
00909
00910 void init()
00911 {
00912 Base::init();
00913
00914 ilu_set_default_options(&m_sluOptions);
00915 m_sluOptions.PrintStat = NO;
00916 m_sluOptions.ConditionNumber = NO;
00917 m_sluOptions.Trans = NOTRANS;
00918 m_sluOptions.ColPerm = MMD_AT_PLUS_A;
00919
00920
00921 m_sluOptions.ILU_MILU = SILU;
00922
00923
00924
00925 m_sluOptions.ILU_DropRule = DROP_BASIC;
00926 m_sluOptions.ILU_DropTol = NumTraits<Scalar>::dummy_precision()*10;
00927 }
00928
00929 private:
00930 SuperILU(SuperILU& ) { }
00931 };
00932
00933 template<typename MatrixType>
00934 void SuperILU<MatrixType>::factorize(const MatrixType& a)
00935 {
00936 eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
00937 if(!m_analysisIsOk)
00938 {
00939 m_info = InvalidInput;
00940 return;
00941 }
00942
00943 this->initFactorization(a);
00944
00945 int info = 0;
00946 RealScalar recip_pivot_growth, rcond;
00947
00948 StatInit(&m_sluStat);
00949 SuperLU_gsisx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0],
00950 &m_sluEqued, &m_sluRscale[0], &m_sluCscale[0],
00951 &m_sluL, &m_sluU,
00952 NULL, 0,
00953 &m_sluB, &m_sluX,
00954 &recip_pivot_growth, &rcond,
00955 &m_sluStat, &info, Scalar());
00956 StatFree(&m_sluStat);
00957
00958
00959 m_info = info == 0 ? Success : NumericalIssue;
00960 m_factorizationIsOk = true;
00961 }
00962
00963 template<typename MatrixType>
00964 template<typename Rhs,typename Dest>
00965 void SuperILU<MatrixType>::_solve(const MatrixBase<Rhs> &b, MatrixBase<Dest>& x) const
00966 {
00967 eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or analyzePattern()/factorize()");
00968
00969 const int size = m_matrix.rows();
00970 const int rhsCols = b.cols();
00971 eigen_assert(size==b.rows());
00972
00973 m_sluOptions.Trans = NOTRANS;
00974 m_sluOptions.Fact = FACTORED;
00975 m_sluOptions.IterRefine = NOREFINE;
00976
00977 m_sluFerr.resize(rhsCols);
00978 m_sluBerr.resize(rhsCols);
00979 m_sluB = SluMatrix::Map(b.const_cast_derived());
00980 m_sluX = SluMatrix::Map(x.derived());
00981
00982 typename Rhs::PlainObject b_cpy;
00983 if(m_sluEqued!='N')
00984 {
00985 b_cpy = b;
00986 m_sluB = SluMatrix::Map(b_cpy.const_cast_derived());
00987 }
00988
00989 int info = 0;
00990 RealScalar recip_pivot_growth, rcond;
00991
00992 StatInit(&m_sluStat);
00993 SuperLU_gsisx(&m_sluOptions, &m_sluA,
00994 m_q.data(), m_p.data(),
00995 &m_sluEtree[0], &m_sluEqued,
00996 &m_sluRscale[0], &m_sluCscale[0],
00997 &m_sluL, &m_sluU,
00998 NULL, 0,
00999 &m_sluB, &m_sluX,
01000 &recip_pivot_growth, &rcond,
01001 &m_sluStat, &info, Scalar());
01002 StatFree(&m_sluStat);
01003
01004 m_info = info==0 ? Success : NumericalIssue;
01005 }
01006 #endif
01007
01008 namespace internal {
01009
01010 template<typename _MatrixType, typename Derived, typename Rhs>
01011 struct solve_retval<SuperLUBase<_MatrixType,Derived>, Rhs>
01012 : solve_retval_base<SuperLUBase<_MatrixType,Derived>, Rhs>
01013 {
01014 typedef SuperLUBase<_MatrixType,Derived> Dec;
01015 EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
01016
01017 template<typename Dest> void evalTo(Dest& dst) const
01018 {
01019 dec().derived()._solve(rhs(),dst);
01020 }
01021 };
01022
01023 template<typename _MatrixType, typename Derived, typename Rhs>
01024 struct sparse_solve_retval<SuperLUBase<_MatrixType,Derived>, Rhs>
01025 : sparse_solve_retval_base<SuperLUBase<_MatrixType,Derived>, Rhs>
01026 {
01027 typedef SuperLUBase<_MatrixType,Derived> Dec;
01028 EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
01029
01030 template<typename Dest> void evalTo(Dest& dst) const
01031 {
01032 dec().derived()._solve(rhs(),dst);
01033 }
01034 };
01035
01036 }
01037
01038 }
01039
01040 #endif // EIGEN_SUPERLUSUPPORT_H