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


win_eigen
Author(s): Daniel Stonier
autogenerated on Mon Oct 6 2014 12:26:39