Back to home page

EIC code displayed by LXR
 
 

    

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

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2008-2011 Gael Guennebaud <gael.guennebaud@inria.fr>
0005 //
0006 // This Source Code Form is subject to the terms of the Mozilla
0007 // Public License v. 2.0. If a copy of the MPL was not distributed
0008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0009 
0010 #ifndef EIGEN_UMFPACKSUPPORT_H
0011 #define EIGEN_UMFPACKSUPPORT_H
0012 
0013 // for compatibility with super old version of umfpack,
0014 // not sure this is really needed, but this is harmless.
0015 #ifndef SuiteSparse_long
0016 #ifdef UF_long
0017 #define SuiteSparse_long UF_long
0018 #else
0019 #error neither SuiteSparse_long nor UF_long are defined
0020 #endif
0021 #endif
0022 
0023 namespace Eigen {
0024 
0025 /* TODO extract L, extract U, compute det, etc... */
0026 
0027 // generic double/complex<double> wrapper functions:
0028 
0029 
0030  // Defaults
0031 inline void umfpack_defaults(double control[UMFPACK_CONTROL], double, int)
0032 { umfpack_di_defaults(control); }
0033 
0034 inline void umfpack_defaults(double control[UMFPACK_CONTROL], std::complex<double>, int)
0035 { umfpack_zi_defaults(control); }
0036 
0037 inline void umfpack_defaults(double control[UMFPACK_CONTROL], double, SuiteSparse_long)
0038 { umfpack_dl_defaults(control); }
0039 
0040 inline void umfpack_defaults(double control[UMFPACK_CONTROL], std::complex<double>, SuiteSparse_long)
0041 { umfpack_zl_defaults(control); }
0042 
0043 // Report info
0044 inline void umfpack_report_info(double control[UMFPACK_CONTROL], double info[UMFPACK_INFO], double, int)
0045 { umfpack_di_report_info(control, info);}
0046 
0047 inline void umfpack_report_info(double control[UMFPACK_CONTROL], double info[UMFPACK_INFO], std::complex<double>, int)
0048 { umfpack_zi_report_info(control, info);}
0049 
0050 inline void umfpack_report_info(double control[UMFPACK_CONTROL], double info[UMFPACK_INFO], double, SuiteSparse_long)
0051 { umfpack_dl_report_info(control, info);}
0052 
0053 inline void umfpack_report_info(double control[UMFPACK_CONTROL], double info[UMFPACK_INFO], std::complex<double>, SuiteSparse_long)
0054 { umfpack_zl_report_info(control, info);}
0055 
0056 // Report status
0057 inline void umfpack_report_status(double control[UMFPACK_CONTROL], int status, double, int)
0058 { umfpack_di_report_status(control, status);}
0059 
0060 inline void umfpack_report_status(double control[UMFPACK_CONTROL], int status, std::complex<double>, int)
0061 { umfpack_zi_report_status(control, status);}
0062 
0063 inline void umfpack_report_status(double control[UMFPACK_CONTROL], int status, double, SuiteSparse_long)
0064 { umfpack_dl_report_status(control, status);}
0065 
0066 inline void umfpack_report_status(double control[UMFPACK_CONTROL], int status, std::complex<double>, SuiteSparse_long)
0067 { umfpack_zl_report_status(control, status);}
0068 
0069 // report control
0070 inline void umfpack_report_control(double control[UMFPACK_CONTROL], double, int)
0071 { umfpack_di_report_control(control);}
0072 
0073 inline void umfpack_report_control(double control[UMFPACK_CONTROL], std::complex<double>, int)
0074 { umfpack_zi_report_control(control);}
0075 
0076 inline void umfpack_report_control(double control[UMFPACK_CONTROL], double, SuiteSparse_long)
0077 { umfpack_dl_report_control(control);}
0078 
0079 inline void umfpack_report_control(double control[UMFPACK_CONTROL], std::complex<double>, SuiteSparse_long)
0080 { umfpack_zl_report_control(control);}
0081 
0082 // Free numeric
0083 inline void umfpack_free_numeric(void **Numeric, double, int)
0084 { umfpack_di_free_numeric(Numeric); *Numeric = 0; }
0085 
0086 inline void umfpack_free_numeric(void **Numeric, std::complex<double>, int)
0087 { umfpack_zi_free_numeric(Numeric); *Numeric = 0; }
0088 
0089 inline void umfpack_free_numeric(void **Numeric, double, SuiteSparse_long)
0090 { umfpack_dl_free_numeric(Numeric); *Numeric = 0; }
0091 
0092 inline void umfpack_free_numeric(void **Numeric, std::complex<double>, SuiteSparse_long)
0093 { umfpack_zl_free_numeric(Numeric); *Numeric = 0; }
0094 
0095 // Free symbolic
0096 inline void umfpack_free_symbolic(void **Symbolic, double, int)
0097 { umfpack_di_free_symbolic(Symbolic); *Symbolic = 0; }
0098 
0099 inline void umfpack_free_symbolic(void **Symbolic, std::complex<double>, int)
0100 { umfpack_zi_free_symbolic(Symbolic); *Symbolic = 0; }
0101 
0102 inline void umfpack_free_symbolic(void **Symbolic, double, SuiteSparse_long)
0103 { umfpack_dl_free_symbolic(Symbolic); *Symbolic = 0; }
0104 
0105 inline void umfpack_free_symbolic(void **Symbolic, std::complex<double>, SuiteSparse_long)
0106 { umfpack_zl_free_symbolic(Symbolic); *Symbolic = 0; }
0107 
0108 // Symbolic
0109 inline int umfpack_symbolic(int n_row,int n_col,
0110                             const int Ap[], const int Ai[], const double Ax[], void **Symbolic,
0111                             const double Control [UMFPACK_CONTROL], double Info [UMFPACK_INFO])
0112 {
0113   return umfpack_di_symbolic(n_row,n_col,Ap,Ai,Ax,Symbolic,Control,Info);
0114 }
0115 
0116 inline int umfpack_symbolic(int n_row,int n_col,
0117                             const int Ap[], const int Ai[], const std::complex<double> Ax[], void **Symbolic,
0118                             const double Control [UMFPACK_CONTROL], double Info [UMFPACK_INFO])
0119 {
0120   return umfpack_zi_symbolic(n_row,n_col,Ap,Ai,&numext::real_ref(Ax[0]),0,Symbolic,Control,Info);
0121 }
0122 inline SuiteSparse_long umfpack_symbolic( SuiteSparse_long n_row,SuiteSparse_long n_col,
0123                                           const SuiteSparse_long Ap[], const SuiteSparse_long Ai[], const double Ax[], void **Symbolic,
0124                                           const double Control [UMFPACK_CONTROL], double Info [UMFPACK_INFO])
0125 {
0126   return umfpack_dl_symbolic(n_row,n_col,Ap,Ai,Ax,Symbolic,Control,Info);
0127 }
0128 
0129 inline SuiteSparse_long umfpack_symbolic( SuiteSparse_long n_row,SuiteSparse_long n_col,
0130                                           const SuiteSparse_long Ap[], const SuiteSparse_long Ai[], const std::complex<double> Ax[], void **Symbolic,
0131                                           const double Control [UMFPACK_CONTROL], double Info [UMFPACK_INFO])
0132 {
0133   return umfpack_zl_symbolic(n_row,n_col,Ap,Ai,&numext::real_ref(Ax[0]),0,Symbolic,Control,Info);
0134 }
0135 
0136 // Numeric
0137 inline int umfpack_numeric( const int Ap[], const int Ai[], const double Ax[],
0138                             void *Symbolic, void **Numeric,
0139                             const double Control[UMFPACK_CONTROL],double Info [UMFPACK_INFO])
0140 {
0141   return umfpack_di_numeric(Ap,Ai,Ax,Symbolic,Numeric,Control,Info);
0142 }
0143 
0144 inline int umfpack_numeric( const int Ap[], const int Ai[], const std::complex<double> Ax[],
0145                             void *Symbolic, void **Numeric,
0146                             const double Control[UMFPACK_CONTROL],double Info [UMFPACK_INFO])
0147 {
0148   return umfpack_zi_numeric(Ap,Ai,&numext::real_ref(Ax[0]),0,Symbolic,Numeric,Control,Info);
0149 }
0150 inline SuiteSparse_long umfpack_numeric(const SuiteSparse_long Ap[], const SuiteSparse_long Ai[], const double Ax[],
0151                                         void *Symbolic, void **Numeric,
0152                                         const double Control[UMFPACK_CONTROL],double Info [UMFPACK_INFO])
0153 {
0154   return umfpack_dl_numeric(Ap,Ai,Ax,Symbolic,Numeric,Control,Info);
0155 }
0156 
0157 inline SuiteSparse_long umfpack_numeric(const SuiteSparse_long Ap[], const SuiteSparse_long Ai[], const std::complex<double> Ax[],
0158                                         void *Symbolic, void **Numeric,
0159                                         const double Control[UMFPACK_CONTROL],double Info [UMFPACK_INFO])
0160 {
0161   return umfpack_zl_numeric(Ap,Ai,&numext::real_ref(Ax[0]),0,Symbolic,Numeric,Control,Info);
0162 }
0163 
0164 // solve
0165 inline int umfpack_solve( int sys, const int Ap[], const int Ai[], const double Ax[],
0166                           double X[], const double B[], void *Numeric,
0167                           const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO])
0168 {
0169   return umfpack_di_solve(sys,Ap,Ai,Ax,X,B,Numeric,Control,Info);
0170 }
0171 
0172 inline int umfpack_solve( int sys, const int Ap[], const int Ai[], const std::complex<double> Ax[],
0173                           std::complex<double> X[], const std::complex<double> B[], void *Numeric,
0174                           const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO])
0175 {
0176   return umfpack_zi_solve(sys,Ap,Ai,&numext::real_ref(Ax[0]),0,&numext::real_ref(X[0]),0,&numext::real_ref(B[0]),0,Numeric,Control,Info);
0177 }
0178 
0179 inline SuiteSparse_long umfpack_solve(int sys, const SuiteSparse_long Ap[], const SuiteSparse_long Ai[], const double Ax[],
0180                                       double X[], const double B[], void *Numeric,
0181                                       const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO])
0182 {
0183   return umfpack_dl_solve(sys,Ap,Ai,Ax,X,B,Numeric,Control,Info);
0184 }
0185 
0186 inline SuiteSparse_long umfpack_solve(int sys, const SuiteSparse_long Ap[], const SuiteSparse_long Ai[], const std::complex<double> Ax[],
0187                                       std::complex<double> X[], const std::complex<double> B[], void *Numeric,
0188                                       const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO])
0189 {
0190   return umfpack_zl_solve(sys,Ap,Ai,&numext::real_ref(Ax[0]),0,&numext::real_ref(X[0]),0,&numext::real_ref(B[0]),0,Numeric,Control,Info);
0191 }
0192 
0193 // Get Lunz
0194 inline int umfpack_get_lunz(int *lnz, int *unz, int *n_row, int *n_col, int *nz_udiag, void *Numeric, double)
0195 {
0196   return umfpack_di_get_lunz(lnz,unz,n_row,n_col,nz_udiag,Numeric);
0197 }
0198 
0199 inline int umfpack_get_lunz(int *lnz, int *unz, int *n_row, int *n_col, int *nz_udiag, void *Numeric, std::complex<double>)
0200 {
0201   return umfpack_zi_get_lunz(lnz,unz,n_row,n_col,nz_udiag,Numeric);
0202 }
0203 
0204 inline SuiteSparse_long umfpack_get_lunz( SuiteSparse_long *lnz, SuiteSparse_long *unz, SuiteSparse_long *n_row, SuiteSparse_long *n_col,
0205                                           SuiteSparse_long *nz_udiag, void *Numeric, double)
0206 {
0207   return umfpack_dl_get_lunz(lnz,unz,n_row,n_col,nz_udiag,Numeric);
0208 }
0209 
0210 inline SuiteSparse_long umfpack_get_lunz( SuiteSparse_long *lnz, SuiteSparse_long *unz, SuiteSparse_long *n_row, SuiteSparse_long *n_col,
0211                                           SuiteSparse_long *nz_udiag, void *Numeric, std::complex<double>)
0212 {
0213   return umfpack_zl_get_lunz(lnz,unz,n_row,n_col,nz_udiag,Numeric);
0214 }
0215 
0216 // Get Numeric
0217 inline int umfpack_get_numeric(int Lp[], int Lj[], double Lx[], int Up[], int Ui[], double Ux[],
0218                                int P[], int Q[], double Dx[], int *do_recip, double Rs[], void *Numeric)
0219 {
0220   return umfpack_di_get_numeric(Lp,Lj,Lx,Up,Ui,Ux,P,Q,Dx,do_recip,Rs,Numeric);
0221 }
0222 
0223 inline int umfpack_get_numeric(int Lp[], int Lj[], std::complex<double> Lx[], int Up[], int Ui[], std::complex<double> Ux[],
0224                                int P[], int Q[], std::complex<double> Dx[], int *do_recip, double Rs[], void *Numeric)
0225 {
0226   double& lx0_real = numext::real_ref(Lx[0]);
0227   double& ux0_real = numext::real_ref(Ux[0]);
0228   double& dx0_real = numext::real_ref(Dx[0]);
0229   return umfpack_zi_get_numeric(Lp,Lj,Lx?&lx0_real:0,0,Up,Ui,Ux?&ux0_real:0,0,P,Q,
0230                                 Dx?&dx0_real:0,0,do_recip,Rs,Numeric);
0231 }
0232 inline SuiteSparse_long umfpack_get_numeric(SuiteSparse_long Lp[], SuiteSparse_long Lj[], double Lx[], SuiteSparse_long Up[], SuiteSparse_long Ui[], double Ux[],
0233                                             SuiteSparse_long P[], SuiteSparse_long Q[], double Dx[], SuiteSparse_long *do_recip, double Rs[], void *Numeric)
0234 {
0235   return umfpack_dl_get_numeric(Lp,Lj,Lx,Up,Ui,Ux,P,Q,Dx,do_recip,Rs,Numeric);
0236 }
0237 
0238 inline SuiteSparse_long umfpack_get_numeric(SuiteSparse_long Lp[], SuiteSparse_long Lj[], std::complex<double> Lx[], SuiteSparse_long Up[], SuiteSparse_long Ui[], std::complex<double> Ux[],
0239                                             SuiteSparse_long P[], SuiteSparse_long Q[], std::complex<double> Dx[], SuiteSparse_long *do_recip, double Rs[], void *Numeric)
0240 {
0241   double& lx0_real = numext::real_ref(Lx[0]);
0242   double& ux0_real = numext::real_ref(Ux[0]);
0243   double& dx0_real = numext::real_ref(Dx[0]);
0244   return umfpack_zl_get_numeric(Lp,Lj,Lx?&lx0_real:0,0,Up,Ui,Ux?&ux0_real:0,0,P,Q,
0245                                 Dx?&dx0_real:0,0,do_recip,Rs,Numeric);
0246 }
0247 
0248 // Get Determinant
0249 inline int umfpack_get_determinant(double *Mx, double *Ex, void *NumericHandle, double User_Info [UMFPACK_INFO], int)
0250 {
0251   return umfpack_di_get_determinant(Mx,Ex,NumericHandle,User_Info);
0252 }
0253 
0254 inline int umfpack_get_determinant(std::complex<double> *Mx, double *Ex, void *NumericHandle, double User_Info [UMFPACK_INFO], int)
0255 {
0256   double& mx_real = numext::real_ref(*Mx);
0257   return umfpack_zi_get_determinant(&mx_real,0,Ex,NumericHandle,User_Info);
0258 }
0259 
0260 inline SuiteSparse_long umfpack_get_determinant(double *Mx, double *Ex, void *NumericHandle, double User_Info [UMFPACK_INFO], SuiteSparse_long)
0261 {
0262   return umfpack_dl_get_determinant(Mx,Ex,NumericHandle,User_Info);
0263 }
0264 
0265 inline SuiteSparse_long umfpack_get_determinant(std::complex<double> *Mx, double *Ex, void *NumericHandle, double User_Info [UMFPACK_INFO], SuiteSparse_long)
0266 {
0267   double& mx_real = numext::real_ref(*Mx);
0268   return umfpack_zl_get_determinant(&mx_real,0,Ex,NumericHandle,User_Info);
0269 }
0270 
0271 
0272 /** \ingroup UmfPackSupport_Module
0273   * \brief A sparse LU factorization and solver based on UmfPack
0274   *
0275   * This class allows to solve for A.X = B sparse linear problems via a LU factorization
0276   * using the UmfPack library. The sparse matrix A must be squared and full rank.
0277   * The vectors or matrices X and B can be either dense or sparse.
0278   *
0279   * \warning The input matrix A should be in a \b compressed and \b column-major form.
0280   * Otherwise an expensive copy will be made. You can call the inexpensive makeCompressed() to get a compressed matrix.
0281   * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<>
0282   *
0283   * \implsparsesolverconcept
0284   *
0285   * \sa \ref TutorialSparseSolverConcept, class SparseLU
0286   */
0287 template<typename _MatrixType>
0288 class UmfPackLU : public SparseSolverBase<UmfPackLU<_MatrixType> >
0289 {
0290   protected:
0291     typedef SparseSolverBase<UmfPackLU<_MatrixType> > Base;
0292     using Base::m_isInitialized;
0293   public:
0294     using Base::_solve_impl;
0295     typedef _MatrixType MatrixType;
0296     typedef typename MatrixType::Scalar Scalar;
0297     typedef typename MatrixType::RealScalar RealScalar;
0298     typedef typename MatrixType::StorageIndex StorageIndex;
0299     typedef Matrix<Scalar,Dynamic,1> Vector;
0300     typedef Matrix<int, 1, MatrixType::ColsAtCompileTime> IntRowVectorType;
0301     typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType;
0302     typedef SparseMatrix<Scalar> LUMatrixType;
0303     typedef SparseMatrix<Scalar,ColMajor,StorageIndex> UmfpackMatrixType;
0304     typedef Ref<const UmfpackMatrixType, StandardCompressedFormat> UmfpackMatrixRef;
0305     enum {
0306       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
0307       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
0308     };
0309 
0310   public:
0311 
0312     typedef Array<double, UMFPACK_CONTROL, 1> UmfpackControl;
0313     typedef Array<double, UMFPACK_INFO, 1> UmfpackInfo;
0314 
0315     UmfPackLU()
0316       : m_dummy(0,0), mp_matrix(m_dummy)
0317     {
0318       init();
0319     }
0320 
0321     template<typename InputMatrixType>
0322     explicit UmfPackLU(const InputMatrixType& matrix)
0323       : mp_matrix(matrix)
0324     {
0325       init();
0326       compute(matrix);
0327     }
0328 
0329     ~UmfPackLU()
0330     {
0331       if(m_symbolic) umfpack_free_symbolic(&m_symbolic,Scalar(), StorageIndex());
0332       if(m_numeric)  umfpack_free_numeric(&m_numeric,Scalar(), StorageIndex());
0333     }
0334 
0335     inline Index rows() const { return mp_matrix.rows(); }
0336     inline Index cols() const { return mp_matrix.cols(); }
0337 
0338     /** \brief Reports whether previous computation was successful.
0339       *
0340       * \returns \c Success if computation was successful,
0341       *          \c NumericalIssue if the matrix.appears to be negative.
0342       */
0343     ComputationInfo info() const
0344     {
0345       eigen_assert(m_isInitialized && "Decomposition is not initialized.");
0346       return m_info;
0347     }
0348 
0349     inline const LUMatrixType& matrixL() const
0350     {
0351       if (m_extractedDataAreDirty) extractData();
0352       return m_l;
0353     }
0354 
0355     inline const LUMatrixType& matrixU() const
0356     {
0357       if (m_extractedDataAreDirty) extractData();
0358       return m_u;
0359     }
0360 
0361     inline const IntColVectorType& permutationP() const
0362     {
0363       if (m_extractedDataAreDirty) extractData();
0364       return m_p;
0365     }
0366 
0367     inline const IntRowVectorType& permutationQ() const
0368     {
0369       if (m_extractedDataAreDirty) extractData();
0370       return m_q;
0371     }
0372 
0373     /** Computes the sparse Cholesky decomposition of \a matrix
0374      *  Note that the matrix should be column-major, and in compressed format for best performance.
0375      *  \sa SparseMatrix::makeCompressed().
0376      */
0377     template<typename InputMatrixType>
0378     void compute(const InputMatrixType& matrix)
0379     {
0380       if(m_symbolic) umfpack_free_symbolic(&m_symbolic,Scalar(),StorageIndex());
0381       if(m_numeric)  umfpack_free_numeric(&m_numeric,Scalar(),StorageIndex());
0382       grab(matrix.derived());
0383       analyzePattern_impl();
0384       factorize_impl();
0385     }
0386 
0387     /** Performs a symbolic decomposition on the sparcity of \a matrix.
0388       *
0389       * This function is particularly useful when solving for several problems having the same structure.
0390       *
0391       * \sa factorize(), compute()
0392       */
0393     template<typename InputMatrixType>
0394     void analyzePattern(const InputMatrixType& matrix)
0395     {
0396       if(m_symbolic) umfpack_free_symbolic(&m_symbolic,Scalar(),StorageIndex());
0397       if(m_numeric)  umfpack_free_numeric(&m_numeric,Scalar(),StorageIndex());
0398 
0399       grab(matrix.derived());
0400 
0401       analyzePattern_impl();
0402     }
0403 
0404     /** Provides the return status code returned by UmfPack during the numeric
0405       * factorization.
0406       *
0407       * \sa factorize(), compute()
0408       */
0409     inline int umfpackFactorizeReturncode() const
0410     {
0411       eigen_assert(m_numeric && "UmfPackLU: you must first call factorize()");
0412       return m_fact_errorCode;
0413     }
0414 
0415     /** Provides access to the control settings array used by UmfPack.
0416       *
0417       * If this array contains NaN's, the default values are used.
0418       *
0419       * See UMFPACK documentation for details.
0420       */
0421     inline const UmfpackControl& umfpackControl() const
0422     {
0423       return m_control;
0424     }
0425 
0426     /** Provides access to the control settings array used by UmfPack.
0427       *
0428       * If this array contains NaN's, the default values are used.
0429       *
0430       * See UMFPACK documentation for details.
0431       */
0432     inline UmfpackControl& umfpackControl()
0433     {
0434       return m_control;
0435     }
0436 
0437     /** Performs a numeric decomposition of \a matrix
0438       *
0439       * The given matrix must has the same sparcity than the matrix on which the pattern anylysis has been performed.
0440       *
0441       * \sa analyzePattern(), compute()
0442       */
0443     template<typename InputMatrixType>
0444     void factorize(const InputMatrixType& matrix)
0445     {
0446       eigen_assert(m_analysisIsOk && "UmfPackLU: you must first call analyzePattern()");
0447       if(m_numeric)
0448         umfpack_free_numeric(&m_numeric,Scalar(),StorageIndex());
0449 
0450       grab(matrix.derived());
0451 
0452       factorize_impl();
0453     }
0454 
0455     /** Prints the current UmfPack control settings.
0456       *
0457       * \sa umfpackControl()
0458       */
0459     void printUmfpackControl()
0460     {
0461       umfpack_report_control(m_control.data(), Scalar(),StorageIndex());
0462     }
0463 
0464     /** Prints statistics collected by UmfPack.
0465       *
0466       * \sa analyzePattern(), compute()
0467       */
0468     void printUmfpackInfo()
0469     {
0470       eigen_assert(m_analysisIsOk && "UmfPackLU: you must first call analyzePattern()");
0471       umfpack_report_info(m_control.data(), m_umfpackInfo.data(), Scalar(),StorageIndex());
0472     }
0473 
0474     /** Prints the status of the previous factorization operation performed by UmfPack (symbolic or numerical factorization).
0475       *
0476       * \sa analyzePattern(), compute()
0477       */
0478     void printUmfpackStatus() {
0479       eigen_assert(m_analysisIsOk && "UmfPackLU: you must first call analyzePattern()");
0480       umfpack_report_status(m_control.data(), m_fact_errorCode, Scalar(),StorageIndex());
0481     }
0482 
0483     /** \internal */
0484     template<typename BDerived,typename XDerived>
0485     bool _solve_impl(const MatrixBase<BDerived> &b, MatrixBase<XDerived> &x) const;
0486 
0487     Scalar determinant() const;
0488 
0489     void extractData() const;
0490 
0491   protected:
0492 
0493     void init()
0494     {
0495       m_info                  = InvalidInput;
0496       m_isInitialized         = false;
0497       m_numeric               = 0;
0498       m_symbolic              = 0;
0499       m_extractedDataAreDirty = true;
0500 
0501       umfpack_defaults(m_control.data(), Scalar(),StorageIndex());
0502     }
0503 
0504     void analyzePattern_impl()
0505     {
0506       m_fact_errorCode = umfpack_symbolic(internal::convert_index<StorageIndex>(mp_matrix.rows()),
0507                                           internal::convert_index<StorageIndex>(mp_matrix.cols()),
0508                                           mp_matrix.outerIndexPtr(), mp_matrix.innerIndexPtr(), mp_matrix.valuePtr(),
0509                                           &m_symbolic, m_control.data(), m_umfpackInfo.data());
0510 
0511       m_isInitialized = true;
0512       m_info = m_fact_errorCode ? InvalidInput : Success;
0513       m_analysisIsOk = true;
0514       m_factorizationIsOk = false;
0515       m_extractedDataAreDirty = true;
0516     }
0517 
0518     void factorize_impl()
0519     {
0520 
0521       m_fact_errorCode = umfpack_numeric(mp_matrix.outerIndexPtr(), mp_matrix.innerIndexPtr(), mp_matrix.valuePtr(),
0522                                          m_symbolic, &m_numeric, m_control.data(), m_umfpackInfo.data());
0523 
0524       m_info = m_fact_errorCode == UMFPACK_OK ? Success : NumericalIssue;
0525       m_factorizationIsOk = true;
0526       m_extractedDataAreDirty = true;
0527     }
0528 
0529     template<typename MatrixDerived>
0530     void grab(const EigenBase<MatrixDerived> &A)
0531     {
0532       mp_matrix.~UmfpackMatrixRef();
0533       ::new (&mp_matrix) UmfpackMatrixRef(A.derived());
0534     }
0535 
0536     void grab(const UmfpackMatrixRef &A)
0537     {
0538       if(&(A.derived()) != &mp_matrix)
0539       {
0540         mp_matrix.~UmfpackMatrixRef();
0541         ::new (&mp_matrix) UmfpackMatrixRef(A);
0542       }
0543     }
0544 
0545     // cached data to reduce reallocation, etc.
0546     mutable LUMatrixType m_l;
0547     StorageIndex m_fact_errorCode;
0548     UmfpackControl m_control;
0549     mutable UmfpackInfo m_umfpackInfo;
0550 
0551     mutable LUMatrixType m_u;
0552     mutable IntColVectorType m_p;
0553     mutable IntRowVectorType m_q;
0554 
0555     UmfpackMatrixType m_dummy;
0556     UmfpackMatrixRef mp_matrix;
0557 
0558     void* m_numeric;
0559     void* m_symbolic;
0560 
0561     mutable ComputationInfo m_info;
0562     int m_factorizationIsOk;
0563     int m_analysisIsOk;
0564     mutable bool m_extractedDataAreDirty;
0565 
0566   private:
0567     UmfPackLU(const UmfPackLU& ) { }
0568 };
0569 
0570 
0571 template<typename MatrixType>
0572 void UmfPackLU<MatrixType>::extractData() const
0573 {
0574   if (m_extractedDataAreDirty)
0575   {
0576     // get size of the data
0577     StorageIndex lnz, unz, rows, cols, nz_udiag;
0578     umfpack_get_lunz(&lnz, &unz, &rows, &cols, &nz_udiag, m_numeric, Scalar());
0579 
0580     // allocate data
0581     m_l.resize(rows,(std::min)(rows,cols));
0582     m_l.resizeNonZeros(lnz);
0583 
0584     m_u.resize((std::min)(rows,cols),cols);
0585     m_u.resizeNonZeros(unz);
0586 
0587     m_p.resize(rows);
0588     m_q.resize(cols);
0589 
0590     // extract
0591     umfpack_get_numeric(m_l.outerIndexPtr(), m_l.innerIndexPtr(), m_l.valuePtr(),
0592                         m_u.outerIndexPtr(), m_u.innerIndexPtr(), m_u.valuePtr(),
0593                         m_p.data(), m_q.data(), 0, 0, 0, m_numeric);
0594 
0595     m_extractedDataAreDirty = false;
0596   }
0597 }
0598 
0599 template<typename MatrixType>
0600 typename UmfPackLU<MatrixType>::Scalar UmfPackLU<MatrixType>::determinant() const
0601 {
0602   Scalar det;
0603   umfpack_get_determinant(&det, 0, m_numeric, 0, StorageIndex());
0604   return det;
0605 }
0606 
0607 template<typename MatrixType>
0608 template<typename BDerived,typename XDerived>
0609 bool UmfPackLU<MatrixType>::_solve_impl(const MatrixBase<BDerived> &b, MatrixBase<XDerived> &x) const
0610 {
0611   Index rhsCols = b.cols();
0612   eigen_assert((BDerived::Flags&RowMajorBit)==0 && "UmfPackLU backend does not support non col-major rhs yet");
0613   eigen_assert((XDerived::Flags&RowMajorBit)==0 && "UmfPackLU backend does not support non col-major result yet");
0614   eigen_assert(b.derived().data() != x.derived().data() && " Umfpack does not support inplace solve");
0615 
0616   Scalar* x_ptr = 0;
0617   Matrix<Scalar,Dynamic,1> x_tmp;
0618   if(x.innerStride()!=1)
0619   {
0620     x_tmp.resize(x.rows());
0621     x_ptr = x_tmp.data();
0622   }
0623   for (int j=0; j<rhsCols; ++j)
0624   {
0625     if(x.innerStride()==1)
0626       x_ptr = &x.col(j).coeffRef(0);
0627     StorageIndex errorCode = umfpack_solve(UMFPACK_A,
0628                                 mp_matrix.outerIndexPtr(), mp_matrix.innerIndexPtr(), mp_matrix.valuePtr(),
0629                                 x_ptr, &b.const_cast_derived().col(j).coeffRef(0),
0630                                 m_numeric, m_control.data(), m_umfpackInfo.data());
0631     if(x.innerStride()!=1)
0632       x.col(j) = x_tmp;
0633     if (errorCode!=0)
0634       return false;
0635   }
0636 
0637   return true;
0638 }
0639 
0640 } // end namespace Eigen
0641 
0642 #endif // EIGEN_UMFPACKSUPPORT_H
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!