Eigen  3.3.3
SuperLUSupport.h
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
 All Classes Functions Variables Typedefs Enumerations Enumerator Friends