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