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