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