Back to home page

EIC code displayed by LXR
 
 

    

Warning, file /include/eigen3/Eigen/src/SuperLUSupport/SuperLUSupport.h was not indexed or was modified since last indexation (in which case cross-reference links may be missing, inaccurate or erroneous).

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2008-2015 Gael Guennebaud <gael.guennebaud@inria.fr>
0005 //
0006 // This Source Code Form is subject to the terms of the Mozilla
0007 // Public License v. 2.0. If a copy of the MPL was not distributed
0008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0009 
0010 #ifndef EIGEN_SUPERLUSUPPORT_H
0011 #define EIGEN_SUPERLUSUPPORT_H
0012 
0013 namespace Eigen {
0014 
0015 #if defined(SUPERLU_MAJOR_VERSION) && (SUPERLU_MAJOR_VERSION >= 5)
0016 #define DECL_GSSVX(PREFIX,FLOATTYPE,KEYTYPE)        \
0017     extern "C" {                                                                                          \
0018       extern void PREFIX##gssvx(superlu_options_t *, SuperMatrix *, int *, int *, int *,                  \
0019                                 char *, FLOATTYPE *, FLOATTYPE *, SuperMatrix *, SuperMatrix *,           \
0020                                 void *, int, SuperMatrix *, SuperMatrix *,                                \
0021                                 FLOATTYPE *, FLOATTYPE *, FLOATTYPE *, FLOATTYPE *,                       \
0022                                 GlobalLU_t *, mem_usage_t *, SuperLUStat_t *, int *);                     \
0023     }                                                                                                     \
0024     inline float SuperLU_gssvx(superlu_options_t *options, SuperMatrix *A,                                \
0025          int *perm_c, int *perm_r, int *etree, char *equed,                                               \
0026          FLOATTYPE *R, FLOATTYPE *C, SuperMatrix *L,                                                      \
0027          SuperMatrix *U, void *work, int lwork,                                                           \
0028          SuperMatrix *B, SuperMatrix *X,                                                                  \
0029          FLOATTYPE *recip_pivot_growth,                                                                   \
0030          FLOATTYPE *rcond, FLOATTYPE *ferr, FLOATTYPE *berr,                                              \
0031          SuperLUStat_t *stats, int *info, KEYTYPE) {                                                      \
0032     mem_usage_t mem_usage;                                                                                \
0033     GlobalLU_t gLU;                                                                                       \
0034     PREFIX##gssvx(options, A, perm_c, perm_r, etree, equed, R, C, L,                                      \
0035          U, work, lwork, B, X, recip_pivot_growth, rcond,                                                 \
0036          ferr, berr, &gLU, &mem_usage, stats, info);                                                      \
0037     return mem_usage.for_lu; /* bytes used by the factor storage */                                       \
0038   }
0039 #else // version < 5.0
0040 #define DECL_GSSVX(PREFIX,FLOATTYPE,KEYTYPE)        \
0041     extern "C" {                                                                                          \
0042       extern void PREFIX##gssvx(superlu_options_t *, SuperMatrix *, int *, int *, int *,                  \
0043                                 char *, FLOATTYPE *, FLOATTYPE *, SuperMatrix *, SuperMatrix *,           \
0044                                 void *, int, SuperMatrix *, SuperMatrix *,                                \
0045                                 FLOATTYPE *, FLOATTYPE *, FLOATTYPE *, FLOATTYPE *,                       \
0046                                 mem_usage_t *, SuperLUStat_t *, int *);                                   \
0047     }                                                                                                     \
0048     inline float SuperLU_gssvx(superlu_options_t *options, SuperMatrix *A,                                \
0049          int *perm_c, int *perm_r, int *etree, char *equed,                                               \
0050          FLOATTYPE *R, FLOATTYPE *C, SuperMatrix *L,                                                      \
0051          SuperMatrix *U, void *work, int lwork,                                                           \
0052          SuperMatrix *B, SuperMatrix *X,                                                                  \
0053          FLOATTYPE *recip_pivot_growth,                                                                   \
0054          FLOATTYPE *rcond, FLOATTYPE *ferr, FLOATTYPE *berr,                                              \
0055          SuperLUStat_t *stats, int *info, KEYTYPE) {                                                      \
0056     mem_usage_t mem_usage;                                                                                \
0057     PREFIX##gssvx(options, A, perm_c, perm_r, etree, equed, R, C, L,                                      \
0058          U, work, lwork, B, X, recip_pivot_growth, rcond,                                                 \
0059          ferr, berr, &mem_usage, stats, info);                                                            \
0060     return mem_usage.for_lu; /* bytes used by the factor storage */                                       \
0061   }
0062 #endif
0063 
0064 DECL_GSSVX(s,float,float)
0065 DECL_GSSVX(c,float,std::complex<float>)
0066 DECL_GSSVX(d,double,double)
0067 DECL_GSSVX(z,double,std::complex<double>)
0068 
0069 #ifdef MILU_ALPHA
0070 #define EIGEN_SUPERLU_HAS_ILU
0071 #endif
0072 
0073 #ifdef EIGEN_SUPERLU_HAS_ILU
0074 
0075 // similarly for the incomplete factorization using gsisx
0076 #define DECL_GSISX(PREFIX,FLOATTYPE,KEYTYPE)                                                    \
0077     extern "C" {                                                                                \
0078       extern void PREFIX##gsisx(superlu_options_t *, SuperMatrix *, int *, int *, int *,        \
0079                          char *, FLOATTYPE *, FLOATTYPE *, SuperMatrix *, SuperMatrix *,        \
0080                          void *, int, SuperMatrix *, SuperMatrix *, FLOATTYPE *, FLOATTYPE *,   \
0081                          mem_usage_t *, SuperLUStat_t *, int *);                        \
0082     }                                                                                           \
0083     inline float SuperLU_gsisx(superlu_options_t *options, SuperMatrix *A,                      \
0084          int *perm_c, int *perm_r, int *etree, char *equed,                                     \
0085          FLOATTYPE *R, FLOATTYPE *C, SuperMatrix *L,                                            \
0086          SuperMatrix *U, void *work, int lwork,                                                 \
0087          SuperMatrix *B, SuperMatrix *X,                                                        \
0088          FLOATTYPE *recip_pivot_growth,                                                         \
0089          FLOATTYPE *rcond,                                                                      \
0090          SuperLUStat_t *stats, int *info, KEYTYPE) {                                            \
0091     mem_usage_t mem_usage;                                                              \
0092     PREFIX##gsisx(options, A, perm_c, perm_r, etree, equed, R, C, L,                            \
0093          U, work, lwork, B, X, recip_pivot_growth, rcond,                                       \
0094          &mem_usage, stats, info);                                                              \
0095     return mem_usage.for_lu; /* bytes used by the factor storage */                             \
0096   }
0097 
0098 DECL_GSISX(s,float,float)
0099 DECL_GSISX(c,float,std::complex<float>)
0100 DECL_GSISX(d,double,double)
0101 DECL_GSISX(z,double,std::complex<double>)
0102 
0103 #endif
0104 
0105 template<typename MatrixType>
0106 struct SluMatrixMapHelper;
0107 
0108 /** \internal
0109   *
0110   * A wrapper class for SuperLU matrices. It supports only compressed sparse matrices
0111   * and dense matrices. Supernodal and other fancy format are not supported by this wrapper.
0112   *
0113   * This wrapper class mainly aims to avoids the need of dynamic allocation of the storage structure.
0114   */
0115 struct SluMatrix : SuperMatrix
0116 {
0117   SluMatrix()
0118   {
0119     Store = &storage;
0120   }
0121 
0122   SluMatrix(const SluMatrix& other)
0123     : SuperMatrix(other)
0124   {
0125     Store = &storage;
0126     storage = other.storage;
0127   }
0128 
0129   SluMatrix& operator=(const SluMatrix& other)
0130   {
0131     SuperMatrix::operator=(static_cast<const SuperMatrix&>(other));
0132     Store = &storage;
0133     storage = other.storage;
0134     return *this;
0135   }
0136 
0137   struct
0138   {
0139     union {int nnz;int lda;};
0140     void *values;
0141     int *innerInd;
0142     int *outerInd;
0143   } storage;
0144 
0145   void setStorageType(Stype_t t)
0146   {
0147     Stype = t;
0148     if (t==SLU_NC || t==SLU_NR || t==SLU_DN)
0149       Store = &storage;
0150     else
0151     {
0152       eigen_assert(false && "storage type not supported");
0153       Store = 0;
0154     }
0155   }
0156 
0157   template<typename Scalar>
0158   void setScalarType()
0159   {
0160     if (internal::is_same<Scalar,float>::value)
0161       Dtype = SLU_S;
0162     else if (internal::is_same<Scalar,double>::value)
0163       Dtype = SLU_D;
0164     else if (internal::is_same<Scalar,std::complex<float> >::value)
0165       Dtype = SLU_C;
0166     else if (internal::is_same<Scalar,std::complex<double> >::value)
0167       Dtype = SLU_Z;
0168     else
0169     {
0170       eigen_assert(false && "Scalar type not supported by SuperLU");
0171     }
0172   }
0173 
0174   template<typename MatrixType>
0175   static SluMatrix Map(MatrixBase<MatrixType>& _mat)
0176   {
0177     MatrixType& mat(_mat.derived());
0178     eigen_assert( ((MatrixType::Flags&RowMajorBit)!=RowMajorBit) && "row-major dense matrices are not supported by SuperLU");
0179     SluMatrix res;
0180     res.setStorageType(SLU_DN);
0181     res.setScalarType<typename MatrixType::Scalar>();
0182     res.Mtype     = SLU_GE;
0183 
0184     res.nrow      = internal::convert_index<int>(mat.rows());
0185     res.ncol      = internal::convert_index<int>(mat.cols());
0186 
0187     res.storage.lda       = internal::convert_index<int>(MatrixType::IsVectorAtCompileTime ? mat.size() : mat.outerStride());
0188     res.storage.values    = (void*)(mat.data());
0189     return res;
0190   }
0191 
0192   template<typename MatrixType>
0193   static SluMatrix Map(SparseMatrixBase<MatrixType>& a_mat)
0194   {
0195     MatrixType &mat(a_mat.derived());
0196     SluMatrix res;
0197     if ((MatrixType::Flags&RowMajorBit)==RowMajorBit)
0198     {
0199       res.setStorageType(SLU_NR);
0200       res.nrow      = internal::convert_index<int>(mat.cols());
0201       res.ncol      = internal::convert_index<int>(mat.rows());
0202     }
0203     else
0204     {
0205       res.setStorageType(SLU_NC);
0206       res.nrow      = internal::convert_index<int>(mat.rows());
0207       res.ncol      = internal::convert_index<int>(mat.cols());
0208     }
0209 
0210     res.Mtype       = SLU_GE;
0211 
0212     res.storage.nnz       = internal::convert_index<int>(mat.nonZeros());
0213     res.storage.values    = mat.valuePtr();
0214     res.storage.innerInd  = mat.innerIndexPtr();
0215     res.storage.outerInd  = mat.outerIndexPtr();
0216 
0217     res.setScalarType<typename MatrixType::Scalar>();
0218 
0219     // FIXME the following is not very accurate
0220     if (int(MatrixType::Flags) & int(Upper))
0221       res.Mtype = SLU_TRU;
0222     if (int(MatrixType::Flags) & int(Lower))
0223       res.Mtype = SLU_TRL;
0224 
0225     eigen_assert(((int(MatrixType::Flags) & int(SelfAdjoint))==0) && "SelfAdjoint matrix shape not supported by SuperLU");
0226 
0227     return res;
0228   }
0229 };
0230 
0231 template<typename Scalar, int Rows, int Cols, int Options, int MRows, int MCols>
0232 struct SluMatrixMapHelper<Matrix<Scalar,Rows,Cols,Options,MRows,MCols> >
0233 {
0234   typedef Matrix<Scalar,Rows,Cols,Options,MRows,MCols> MatrixType;
0235   static void run(MatrixType& mat, SluMatrix& res)
0236   {
0237     eigen_assert( ((Options&RowMajor)!=RowMajor) && "row-major dense matrices is not supported by SuperLU");
0238     res.setStorageType(SLU_DN);
0239     res.setScalarType<Scalar>();
0240     res.Mtype     = SLU_GE;
0241 
0242     res.nrow      = mat.rows();
0243     res.ncol      = mat.cols();
0244 
0245     res.storage.lda       = mat.outerStride();
0246     res.storage.values    = mat.data();
0247   }
0248 };
0249 
0250 template<typename Derived>
0251 struct SluMatrixMapHelper<SparseMatrixBase<Derived> >
0252 {
0253   typedef Derived MatrixType;
0254   static void run(MatrixType& mat, SluMatrix& res)
0255   {
0256     if ((MatrixType::Flags&RowMajorBit)==RowMajorBit)
0257     {
0258       res.setStorageType(SLU_NR);
0259       res.nrow      = mat.cols();
0260       res.ncol      = mat.rows();
0261     }
0262     else
0263     {
0264       res.setStorageType(SLU_NC);
0265       res.nrow      = mat.rows();
0266       res.ncol      = mat.cols();
0267     }
0268 
0269     res.Mtype       = SLU_GE;
0270 
0271     res.storage.nnz       = mat.nonZeros();
0272     res.storage.values    = mat.valuePtr();
0273     res.storage.innerInd  = mat.innerIndexPtr();
0274     res.storage.outerInd  = mat.outerIndexPtr();
0275 
0276     res.setScalarType<typename MatrixType::Scalar>();
0277 
0278     // FIXME the following is not very accurate
0279     if (MatrixType::Flags & Upper)
0280       res.Mtype = SLU_TRU;
0281     if (MatrixType::Flags & Lower)
0282       res.Mtype = SLU_TRL;
0283 
0284     eigen_assert(((MatrixType::Flags & SelfAdjoint)==0) && "SelfAdjoint matrix shape not supported by SuperLU");
0285   }
0286 };
0287 
0288 namespace internal {
0289 
0290 template<typename MatrixType>
0291 SluMatrix asSluMatrix(MatrixType& mat)
0292 {
0293   return SluMatrix::Map(mat);
0294 }
0295 
0296 /** View a Super LU matrix as an Eigen expression */
0297 template<typename Scalar, int Flags, typename Index>
0298 MappedSparseMatrix<Scalar,Flags,Index> map_superlu(SluMatrix& sluMat)
0299 {
0300   eigen_assert(((Flags&RowMajor)==RowMajor && sluMat.Stype == SLU_NR)
0301          || ((Flags&ColMajor)==ColMajor && sluMat.Stype == SLU_NC));
0302 
0303   Index outerSize = (Flags&RowMajor)==RowMajor ? sluMat.ncol : sluMat.nrow;
0304 
0305   return MappedSparseMatrix<Scalar,Flags,Index>(
0306     sluMat.nrow, sluMat.ncol, sluMat.storage.outerInd[outerSize],
0307     sluMat.storage.outerInd, sluMat.storage.innerInd, reinterpret_cast<Scalar*>(sluMat.storage.values) );
0308 }
0309 
0310 } // end namespace internal
0311 
0312 /** \ingroup SuperLUSupport_Module
0313   * \class SuperLUBase
0314   * \brief The base class for the direct and incomplete LU factorization of SuperLU
0315   */
0316 template<typename _MatrixType, typename Derived>
0317 class SuperLUBase : public SparseSolverBase<Derived>
0318 {
0319   protected:
0320     typedef SparseSolverBase<Derived> Base;
0321     using Base::derived;
0322     using Base::m_isInitialized;
0323   public:
0324     typedef _MatrixType MatrixType;
0325     typedef typename MatrixType::Scalar Scalar;
0326     typedef typename MatrixType::RealScalar RealScalar;
0327     typedef typename MatrixType::StorageIndex StorageIndex;
0328     typedef Matrix<Scalar,Dynamic,1> Vector;
0329     typedef Matrix<int, 1, MatrixType::ColsAtCompileTime> IntRowVectorType;
0330     typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType;    
0331     typedef Map<PermutationMatrix<Dynamic,Dynamic,int> > PermutationMap;
0332     typedef SparseMatrix<Scalar> LUMatrixType;
0333     enum {
0334       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
0335       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
0336     };
0337 
0338   public:
0339 
0340     SuperLUBase() {}
0341 
0342     ~SuperLUBase()
0343     {
0344       clearFactors();
0345     }
0346     
0347     inline Index rows() const { return m_matrix.rows(); }
0348     inline Index cols() const { return m_matrix.cols(); }
0349     
0350     /** \returns a reference to the Super LU option object to configure the  Super LU algorithms. */
0351     inline superlu_options_t& options() { return m_sluOptions; }
0352     
0353     /** \brief Reports whether previous computation was successful.
0354       *
0355       * \returns \c Success if computation was successful,
0356       *          \c NumericalIssue if the matrix.appears to be negative.
0357       */
0358     ComputationInfo info() const
0359     {
0360       eigen_assert(m_isInitialized && "Decomposition is not initialized.");
0361       return m_info;
0362     }
0363 
0364     /** Computes the sparse Cholesky decomposition of \a matrix */
0365     void compute(const MatrixType& matrix)
0366     {
0367       derived().analyzePattern(matrix);
0368       derived().factorize(matrix);
0369     }
0370 
0371     /** Performs a symbolic decomposition on the sparcity of \a matrix.
0372       *
0373       * This function is particularly useful when solving for several problems having the same structure.
0374       * 
0375       * \sa factorize()
0376       */
0377     void analyzePattern(const MatrixType& /*matrix*/)
0378     {
0379       m_isInitialized = true;
0380       m_info = Success;
0381       m_analysisIsOk = true;
0382       m_factorizationIsOk = false;
0383     }
0384     
0385     template<typename Stream>
0386     void dumpMemory(Stream& /*s*/)
0387     {}
0388     
0389   protected:
0390     
0391     void initFactorization(const MatrixType& a)
0392     {
0393       set_default_options(&this->m_sluOptions);
0394       
0395       const Index size = a.rows();
0396       m_matrix = a;
0397 
0398       m_sluA = internal::asSluMatrix(m_matrix);
0399       clearFactors();
0400 
0401       m_p.resize(size);
0402       m_q.resize(size);
0403       m_sluRscale.resize(size);
0404       m_sluCscale.resize(size);
0405       m_sluEtree.resize(size);
0406 
0407       // set empty B and X
0408       m_sluB.setStorageType(SLU_DN);
0409       m_sluB.setScalarType<Scalar>();
0410       m_sluB.Mtype          = SLU_GE;
0411       m_sluB.storage.values = 0;
0412       m_sluB.nrow           = 0;
0413       m_sluB.ncol           = 0;
0414       m_sluB.storage.lda    = internal::convert_index<int>(size);
0415       m_sluX                = m_sluB;
0416       
0417       m_extractedDataAreDirty = true;
0418     }
0419     
0420     void init()
0421     {
0422       m_info = InvalidInput;
0423       m_isInitialized = false;
0424       m_sluL.Store = 0;
0425       m_sluU.Store = 0;
0426     }
0427     
0428     void extractData() const;
0429 
0430     void clearFactors()
0431     {
0432       if(m_sluL.Store)
0433         Destroy_SuperNode_Matrix(&m_sluL);
0434       if(m_sluU.Store)
0435         Destroy_CompCol_Matrix(&m_sluU);
0436 
0437       m_sluL.Store = 0;
0438       m_sluU.Store = 0;
0439 
0440       memset(&m_sluL,0,sizeof m_sluL);
0441       memset(&m_sluU,0,sizeof m_sluU);
0442     }
0443 
0444     // cached data to reduce reallocation, etc.
0445     mutable LUMatrixType m_l;
0446     mutable LUMatrixType m_u;
0447     mutable IntColVectorType m_p;
0448     mutable IntRowVectorType m_q;
0449 
0450     mutable LUMatrixType m_matrix;  // copy of the factorized matrix
0451     mutable SluMatrix m_sluA;
0452     mutable SuperMatrix m_sluL, m_sluU;
0453     mutable SluMatrix m_sluB, m_sluX;
0454     mutable SuperLUStat_t m_sluStat;
0455     mutable superlu_options_t m_sluOptions;
0456     mutable std::vector<int> m_sluEtree;
0457     mutable Matrix<RealScalar,Dynamic,1> m_sluRscale, m_sluCscale;
0458     mutable Matrix<RealScalar,Dynamic,1> m_sluFerr, m_sluBerr;
0459     mutable char m_sluEqued;
0460 
0461     mutable ComputationInfo m_info;
0462     int m_factorizationIsOk;
0463     int m_analysisIsOk;
0464     mutable bool m_extractedDataAreDirty;
0465     
0466   private:
0467     SuperLUBase(SuperLUBase& ) { }
0468 };
0469 
0470 
0471 /** \ingroup SuperLUSupport_Module
0472   * \class SuperLU
0473   * \brief A sparse direct LU factorization and solver based on the SuperLU library
0474   *
0475   * This class allows to solve for A.X = B sparse linear problems via a direct LU factorization
0476   * using the SuperLU library. The sparse matrix A must be squared and invertible. The vectors or matrices
0477   * X and B can be either dense or sparse.
0478   *
0479   * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<>
0480   *
0481   * \warning This class is only for the 4.x versions of SuperLU. The 3.x and 5.x versions are not supported.
0482   *
0483   * \implsparsesolverconcept
0484   *
0485   * \sa \ref TutorialSparseSolverConcept, class SparseLU
0486   */
0487 template<typename _MatrixType>
0488 class SuperLU : public SuperLUBase<_MatrixType,SuperLU<_MatrixType> >
0489 {
0490   public:
0491     typedef SuperLUBase<_MatrixType,SuperLU> Base;
0492     typedef _MatrixType MatrixType;
0493     typedef typename Base::Scalar Scalar;
0494     typedef typename Base::RealScalar RealScalar;
0495     typedef typename Base::StorageIndex StorageIndex;
0496     typedef typename Base::IntRowVectorType IntRowVectorType;
0497     typedef typename Base::IntColVectorType IntColVectorType;   
0498     typedef typename Base::PermutationMap PermutationMap;
0499     typedef typename Base::LUMatrixType LUMatrixType;
0500     typedef TriangularView<LUMatrixType, Lower|UnitDiag>  LMatrixType;
0501     typedef TriangularView<LUMatrixType,  Upper>          UMatrixType;
0502 
0503   public:
0504     using Base::_solve_impl;
0505 
0506     SuperLU() : Base() { init(); }
0507 
0508     explicit SuperLU(const MatrixType& matrix) : Base()
0509     {
0510       init();
0511       Base::compute(matrix);
0512     }
0513 
0514     ~SuperLU()
0515     {
0516     }
0517     
0518     /** Performs a symbolic decomposition on the sparcity of \a matrix.
0519       *
0520       * This function is particularly useful when solving for several problems having the same structure.
0521       * 
0522       * \sa factorize()
0523       */
0524     void analyzePattern(const MatrixType& matrix)
0525     {
0526       m_info = InvalidInput;
0527       m_isInitialized = false;
0528       Base::analyzePattern(matrix);
0529     }
0530     
0531     /** Performs a numeric decomposition of \a matrix
0532       *
0533       * The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been performed.
0534       *
0535       * \sa analyzePattern()
0536       */
0537     void factorize(const MatrixType& matrix);
0538     
0539     /** \internal */
0540     template<typename Rhs,typename Dest>
0541     void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const;
0542     
0543     inline const LMatrixType& matrixL() const
0544     {
0545       if (m_extractedDataAreDirty) this->extractData();
0546       return m_l;
0547     }
0548 
0549     inline const UMatrixType& matrixU() const
0550     {
0551       if (m_extractedDataAreDirty) this->extractData();
0552       return m_u;
0553     }
0554 
0555     inline const IntColVectorType& permutationP() const
0556     {
0557       if (m_extractedDataAreDirty) this->extractData();
0558       return m_p;
0559     }
0560 
0561     inline const IntRowVectorType& permutationQ() const
0562     {
0563       if (m_extractedDataAreDirty) this->extractData();
0564       return m_q;
0565     }
0566     
0567     Scalar determinant() const;
0568     
0569   protected:
0570     
0571     using Base::m_matrix;
0572     using Base::m_sluOptions;
0573     using Base::m_sluA;
0574     using Base::m_sluB;
0575     using Base::m_sluX;
0576     using Base::m_p;
0577     using Base::m_q;
0578     using Base::m_sluEtree;
0579     using Base::m_sluEqued;
0580     using Base::m_sluRscale;
0581     using Base::m_sluCscale;
0582     using Base::m_sluL;
0583     using Base::m_sluU;
0584     using Base::m_sluStat;
0585     using Base::m_sluFerr;
0586     using Base::m_sluBerr;
0587     using Base::m_l;
0588     using Base::m_u;
0589     
0590     using Base::m_analysisIsOk;
0591     using Base::m_factorizationIsOk;
0592     using Base::m_extractedDataAreDirty;
0593     using Base::m_isInitialized;
0594     using Base::m_info;
0595     
0596     void init()
0597     {
0598       Base::init();
0599       
0600       set_default_options(&this->m_sluOptions);
0601       m_sluOptions.PrintStat        = NO;
0602       m_sluOptions.ConditionNumber  = NO;
0603       m_sluOptions.Trans            = NOTRANS;
0604       m_sluOptions.ColPerm          = COLAMD;
0605     }
0606     
0607     
0608   private:
0609     SuperLU(SuperLU& ) { }
0610 };
0611 
0612 template<typename MatrixType>
0613 void SuperLU<MatrixType>::factorize(const MatrixType& a)
0614 {
0615   eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
0616   if(!m_analysisIsOk)
0617   {
0618     m_info = InvalidInput;
0619     return;
0620   }
0621   
0622   this->initFactorization(a);
0623   
0624   m_sluOptions.ColPerm = COLAMD;
0625   int info = 0;
0626   RealScalar recip_pivot_growth, rcond;
0627   RealScalar ferr, berr;
0628 
0629   StatInit(&m_sluStat);
0630   SuperLU_gssvx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0],
0631                 &m_sluEqued, &m_sluRscale[0], &m_sluCscale[0],
0632                 &m_sluL, &m_sluU,
0633                 NULL, 0,
0634                 &m_sluB, &m_sluX,
0635                 &recip_pivot_growth, &rcond,
0636                 &ferr, &berr,
0637                 &m_sluStat, &info, Scalar());
0638   StatFree(&m_sluStat);
0639 
0640   m_extractedDataAreDirty = true;
0641 
0642   // FIXME how to better check for errors ???
0643   m_info = info == 0 ? Success : NumericalIssue;
0644   m_factorizationIsOk = true;
0645 }
0646 
0647 template<typename MatrixType>
0648 template<typename Rhs,typename Dest>
0649 void SuperLU<MatrixType>::_solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest>& x) const
0650 {
0651   eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or analyzePattern()/factorize()");
0652 
0653   const Index rhsCols = b.cols();
0654   eigen_assert(m_matrix.rows()==b.rows());
0655 
0656   m_sluOptions.Trans = NOTRANS;
0657   m_sluOptions.Fact = FACTORED;
0658   m_sluOptions.IterRefine = NOREFINE;
0659   
0660 
0661   m_sluFerr.resize(rhsCols);
0662   m_sluBerr.resize(rhsCols);
0663   
0664   Ref<const Matrix<typename Rhs::Scalar,Dynamic,Dynamic,ColMajor> > b_ref(b);
0665   Ref<const Matrix<typename Dest::Scalar,Dynamic,Dynamic,ColMajor> > x_ref(x);
0666   
0667   m_sluB = SluMatrix::Map(b_ref.const_cast_derived());
0668   m_sluX = SluMatrix::Map(x_ref.const_cast_derived());
0669   
0670   typename Rhs::PlainObject b_cpy;
0671   if(m_sluEqued!='N')
0672   {
0673     b_cpy = b;
0674     m_sluB = SluMatrix::Map(b_cpy.const_cast_derived());  
0675   }
0676 
0677   StatInit(&m_sluStat);
0678   int info = 0;
0679   RealScalar recip_pivot_growth, rcond;
0680   SuperLU_gssvx(&m_sluOptions, &m_sluA,
0681                 m_q.data(), m_p.data(),
0682                 &m_sluEtree[0], &m_sluEqued,
0683                 &m_sluRscale[0], &m_sluCscale[0],
0684                 &m_sluL, &m_sluU,
0685                 NULL, 0,
0686                 &m_sluB, &m_sluX,
0687                 &recip_pivot_growth, &rcond,
0688                 &m_sluFerr[0], &m_sluBerr[0],
0689                 &m_sluStat, &info, Scalar());
0690   StatFree(&m_sluStat);
0691   
0692   if(x.derived().data() != x_ref.data())
0693     x = x_ref;
0694   
0695   m_info = info==0 ? Success : NumericalIssue;
0696 }
0697 
0698 // the code of this extractData() function has been adapted from the SuperLU's Matlab support code,
0699 //
0700 //  Copyright (c) 1994 by Xerox Corporation.  All rights reserved.
0701 //
0702 //  THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
0703 //  EXPRESSED OR IMPLIED.  ANY USE IS AT YOUR OWN RISK.
0704 //
0705 template<typename MatrixType, typename Derived>
0706 void SuperLUBase<MatrixType,Derived>::extractData() const
0707 {
0708   eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for extracting factors, you must first call either compute() or analyzePattern()/factorize()");
0709   if (m_extractedDataAreDirty)
0710   {
0711     int         upper;
0712     int         fsupc, istart, nsupr;
0713     int         lastl = 0, lastu = 0;
0714     SCformat    *Lstore = static_cast<SCformat*>(m_sluL.Store);
0715     NCformat    *Ustore = static_cast<NCformat*>(m_sluU.Store);
0716     Scalar      *SNptr;
0717 
0718     const Index size = m_matrix.rows();
0719     m_l.resize(size,size);
0720     m_l.resizeNonZeros(Lstore->nnz);
0721     m_u.resize(size,size);
0722     m_u.resizeNonZeros(Ustore->nnz);
0723 
0724     int* Lcol = m_l.outerIndexPtr();
0725     int* Lrow = m_l.innerIndexPtr();
0726     Scalar* Lval = m_l.valuePtr();
0727 
0728     int* Ucol = m_u.outerIndexPtr();
0729     int* Urow = m_u.innerIndexPtr();
0730     Scalar* Uval = m_u.valuePtr();
0731 
0732     Ucol[0] = 0;
0733     Ucol[0] = 0;
0734 
0735     /* for each supernode */
0736     for (int k = 0; k <= Lstore->nsuper; ++k)
0737     {
0738       fsupc   = L_FST_SUPC(k);
0739       istart  = L_SUB_START(fsupc);
0740       nsupr   = L_SUB_START(fsupc+1) - istart;
0741       upper   = 1;
0742 
0743       /* for each column in the supernode */
0744       for (int j = fsupc; j < L_FST_SUPC(k+1); ++j)
0745       {
0746         SNptr = &((Scalar*)Lstore->nzval)[L_NZ_START(j)];
0747 
0748         /* Extract U */
0749         for (int i = U_NZ_START(j); i < U_NZ_START(j+1); ++i)
0750         {
0751           Uval[lastu] = ((Scalar*)Ustore->nzval)[i];
0752           /* Matlab doesn't like explicit zero. */
0753           if (Uval[lastu] != 0.0)
0754             Urow[lastu++] = U_SUB(i);
0755         }
0756         for (int i = 0; i < upper; ++i)
0757         {
0758           /* upper triangle in the supernode */
0759           Uval[lastu] = SNptr[i];
0760           /* Matlab doesn't like explicit zero. */
0761           if (Uval[lastu] != 0.0)
0762             Urow[lastu++] = L_SUB(istart+i);
0763         }
0764         Ucol[j+1] = lastu;
0765 
0766         /* Extract L */
0767         Lval[lastl] = 1.0; /* unit diagonal */
0768         Lrow[lastl++] = L_SUB(istart + upper - 1);
0769         for (int i = upper; i < nsupr; ++i)
0770         {
0771           Lval[lastl] = SNptr[i];
0772           /* Matlab doesn't like explicit zero. */
0773           if (Lval[lastl] != 0.0)
0774             Lrow[lastl++] = L_SUB(istart+i);
0775         }
0776         Lcol[j+1] = lastl;
0777 
0778         ++upper;
0779       } /* for j ... */
0780 
0781     } /* for k ... */
0782 
0783     // squeeze the matrices :
0784     m_l.resizeNonZeros(lastl);
0785     m_u.resizeNonZeros(lastu);
0786 
0787     m_extractedDataAreDirty = false;
0788   }
0789 }
0790 
0791 template<typename MatrixType>
0792 typename SuperLU<MatrixType>::Scalar SuperLU<MatrixType>::determinant() const
0793 {
0794   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()");
0795   
0796   if (m_extractedDataAreDirty)
0797     this->extractData();
0798 
0799   Scalar det = Scalar(1);
0800   for (int j=0; j<m_u.cols(); ++j)
0801   {
0802     if (m_u.outerIndexPtr()[j+1]-m_u.outerIndexPtr()[j] > 0)
0803     {
0804       int lastId = m_u.outerIndexPtr()[j+1]-1;
0805       eigen_assert(m_u.innerIndexPtr()[lastId]<=j);
0806       if (m_u.innerIndexPtr()[lastId]==j)
0807         det *= m_u.valuePtr()[lastId];
0808     }
0809   }
0810   if(PermutationMap(m_p.data(),m_p.size()).determinant()*PermutationMap(m_q.data(),m_q.size()).determinant()<0)
0811     det = -det;
0812   if(m_sluEqued!='N')
0813     return det/m_sluRscale.prod()/m_sluCscale.prod();
0814   else
0815     return det;
0816 }
0817 
0818 #ifdef EIGEN_PARSED_BY_DOXYGEN
0819 #define EIGEN_SUPERLU_HAS_ILU
0820 #endif
0821 
0822 #ifdef EIGEN_SUPERLU_HAS_ILU
0823 
0824 /** \ingroup SuperLUSupport_Module
0825   * \class SuperILU
0826   * \brief A sparse direct \b incomplete LU factorization and solver based on the SuperLU library
0827   *
0828   * This class allows to solve for an approximate solution of A.X = B sparse linear problems via an incomplete LU factorization
0829   * using the SuperLU library. This class is aimed to be used as a preconditioner of the iterative linear solvers.
0830   *
0831   * \warning This class is only for the 4.x versions of SuperLU. The 3.x and 5.x versions are not supported.
0832   *
0833   * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<>
0834   *
0835   * \implsparsesolverconcept
0836   *
0837   * \sa \ref TutorialSparseSolverConcept, class IncompleteLUT, class ConjugateGradient, class BiCGSTAB
0838   */
0839 
0840 template<typename _MatrixType>
0841 class SuperILU : public SuperLUBase<_MatrixType,SuperILU<_MatrixType> >
0842 {
0843   public:
0844     typedef SuperLUBase<_MatrixType,SuperILU> Base;
0845     typedef _MatrixType MatrixType;
0846     typedef typename Base::Scalar Scalar;
0847     typedef typename Base::RealScalar RealScalar;
0848 
0849   public:
0850     using Base::_solve_impl;
0851 
0852     SuperILU() : Base() { init(); }
0853 
0854     SuperILU(const MatrixType& matrix) : Base()
0855     {
0856       init();
0857       Base::compute(matrix);
0858     }
0859 
0860     ~SuperILU()
0861     {
0862     }
0863     
0864     /** Performs a symbolic decomposition on the sparcity of \a matrix.
0865       *
0866       * This function is particularly useful when solving for several problems having the same structure.
0867       * 
0868       * \sa factorize()
0869       */
0870     void analyzePattern(const MatrixType& matrix)
0871     {
0872       Base::analyzePattern(matrix);
0873     }
0874     
0875     /** Performs a numeric decomposition of \a matrix
0876       *
0877       * The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been performed.
0878       *
0879       * \sa analyzePattern()
0880       */
0881     void factorize(const MatrixType& matrix);
0882     
0883     #ifndef EIGEN_PARSED_BY_DOXYGEN
0884     /** \internal */
0885     template<typename Rhs,typename Dest>
0886     void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const;
0887     #endif // EIGEN_PARSED_BY_DOXYGEN
0888     
0889   protected:
0890     
0891     using Base::m_matrix;
0892     using Base::m_sluOptions;
0893     using Base::m_sluA;
0894     using Base::m_sluB;
0895     using Base::m_sluX;
0896     using Base::m_p;
0897     using Base::m_q;
0898     using Base::m_sluEtree;
0899     using Base::m_sluEqued;
0900     using Base::m_sluRscale;
0901     using Base::m_sluCscale;
0902     using Base::m_sluL;
0903     using Base::m_sluU;
0904     using Base::m_sluStat;
0905     using Base::m_sluFerr;
0906     using Base::m_sluBerr;
0907     using Base::m_l;
0908     using Base::m_u;
0909     
0910     using Base::m_analysisIsOk;
0911     using Base::m_factorizationIsOk;
0912     using Base::m_extractedDataAreDirty;
0913     using Base::m_isInitialized;
0914     using Base::m_info;
0915 
0916     void init()
0917     {
0918       Base::init();
0919       
0920       ilu_set_default_options(&m_sluOptions);
0921       m_sluOptions.PrintStat        = NO;
0922       m_sluOptions.ConditionNumber  = NO;
0923       m_sluOptions.Trans            = NOTRANS;
0924       m_sluOptions.ColPerm          = MMD_AT_PLUS_A;
0925       
0926       // no attempt to preserve column sum
0927       m_sluOptions.ILU_MILU = SILU;
0928       // only basic ILU(k) support -- no direct control over memory consumption
0929       // better to use ILU_DropRule = DROP_BASIC | DROP_AREA
0930       // and set ILU_FillFactor to max memory growth
0931       m_sluOptions.ILU_DropRule = DROP_BASIC;
0932       m_sluOptions.ILU_DropTol = NumTraits<Scalar>::dummy_precision()*10;
0933     }
0934     
0935   private:
0936     SuperILU(SuperILU& ) { }
0937 };
0938 
0939 template<typename MatrixType>
0940 void SuperILU<MatrixType>::factorize(const MatrixType& a)
0941 {
0942   eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
0943   if(!m_analysisIsOk)
0944   {
0945     m_info = InvalidInput;
0946     return;
0947   }
0948   
0949   this->initFactorization(a);
0950 
0951   int info = 0;
0952   RealScalar recip_pivot_growth, rcond;
0953 
0954   StatInit(&m_sluStat);
0955   SuperLU_gsisx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0],
0956                 &m_sluEqued, &m_sluRscale[0], &m_sluCscale[0],
0957                 &m_sluL, &m_sluU,
0958                 NULL, 0,
0959                 &m_sluB, &m_sluX,
0960                 &recip_pivot_growth, &rcond,
0961                 &m_sluStat, &info, Scalar());
0962   StatFree(&m_sluStat);
0963 
0964   // FIXME how to better check for errors ???
0965   m_info = info == 0 ? Success : NumericalIssue;
0966   m_factorizationIsOk = true;
0967 }
0968 
0969 #ifndef EIGEN_PARSED_BY_DOXYGEN
0970 template<typename MatrixType>
0971 template<typename Rhs,typename Dest>
0972 void SuperILU<MatrixType>::_solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest>& x) const
0973 {
0974   eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or analyzePattern()/factorize()");
0975 
0976   const int rhsCols = b.cols();
0977   eigen_assert(m_matrix.rows()==b.rows());
0978 
0979   m_sluOptions.Trans = NOTRANS;
0980   m_sluOptions.Fact = FACTORED;
0981   m_sluOptions.IterRefine = NOREFINE;
0982 
0983   m_sluFerr.resize(rhsCols);
0984   m_sluBerr.resize(rhsCols);
0985   
0986   Ref<const Matrix<typename Rhs::Scalar,Dynamic,Dynamic,ColMajor> > b_ref(b);
0987   Ref<const Matrix<typename Dest::Scalar,Dynamic,Dynamic,ColMajor> > x_ref(x);
0988   
0989   m_sluB = SluMatrix::Map(b_ref.const_cast_derived());
0990   m_sluX = SluMatrix::Map(x_ref.const_cast_derived());
0991 
0992   typename Rhs::PlainObject b_cpy;
0993   if(m_sluEqued!='N')
0994   {
0995     b_cpy = b;
0996     m_sluB = SluMatrix::Map(b_cpy.const_cast_derived());  
0997   }
0998   
0999   int info = 0;
1000   RealScalar recip_pivot_growth, rcond;
1001 
1002   StatInit(&m_sluStat);
1003   SuperLU_gsisx(&m_sluOptions, &m_sluA,
1004                 m_q.data(), m_p.data(),
1005                 &m_sluEtree[0], &m_sluEqued,
1006                 &m_sluRscale[0], &m_sluCscale[0],
1007                 &m_sluL, &m_sluU,
1008                 NULL, 0,
1009                 &m_sluB, &m_sluX,
1010                 &recip_pivot_growth, &rcond,
1011                 &m_sluStat, &info, Scalar());
1012   StatFree(&m_sluStat);
1013   
1014   if(x.derived().data() != x_ref.data())
1015     x = x_ref;
1016 
1017   m_info = info==0 ? Success : NumericalIssue;
1018 }
1019 #endif
1020 
1021 #endif
1022 
1023 } // end namespace Eigen
1024 
1025 #endif // EIGEN_SUPERLUSUPPORT_H
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!