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


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Thu Aug 27 2015 12:01:09