$treeview $search $mathjax
Eigen
3.2.5
$projectbrief
|
$projectbrief
|
$searchbox |
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_UMFPACKSUPPORT_H 00011 #define EIGEN_UMFPACKSUPPORT_H 00012 00013 namespace Eigen { 00014 00015 /* TODO extract L, extract U, compute det, etc... */ 00016 00017 // generic double/complex<double> wrapper functions: 00018 00019 inline void umfpack_free_numeric(void **Numeric, double) 00020 { umfpack_di_free_numeric(Numeric); *Numeric = 0; } 00021 00022 inline void umfpack_free_numeric(void **Numeric, std::complex<double>) 00023 { umfpack_zi_free_numeric(Numeric); *Numeric = 0; } 00024 00025 inline void umfpack_free_symbolic(void **Symbolic, double) 00026 { umfpack_di_free_symbolic(Symbolic); *Symbolic = 0; } 00027 00028 inline void umfpack_free_symbolic(void **Symbolic, std::complex<double>) 00029 { umfpack_zi_free_symbolic(Symbolic); *Symbolic = 0; } 00030 00031 inline int umfpack_symbolic(int n_row,int n_col, 00032 const int Ap[], const int Ai[], const double Ax[], void **Symbolic, 00033 const double Control [UMFPACK_CONTROL], double Info [UMFPACK_INFO]) 00034 { 00035 return umfpack_di_symbolic(n_row,n_col,Ap,Ai,Ax,Symbolic,Control,Info); 00036 } 00037 00038 inline int umfpack_symbolic(int n_row,int n_col, 00039 const int Ap[], const int Ai[], const std::complex<double> Ax[], void **Symbolic, 00040 const double Control [UMFPACK_CONTROL], double Info [UMFPACK_INFO]) 00041 { 00042 return umfpack_zi_symbolic(n_row,n_col,Ap,Ai,&numext::real_ref(Ax[0]),0,Symbolic,Control,Info); 00043 } 00044 00045 inline int umfpack_numeric( const int Ap[], const int Ai[], const double Ax[], 00046 void *Symbolic, void **Numeric, 00047 const double Control[UMFPACK_CONTROL],double Info [UMFPACK_INFO]) 00048 { 00049 return umfpack_di_numeric(Ap,Ai,Ax,Symbolic,Numeric,Control,Info); 00050 } 00051 00052 inline int umfpack_numeric( const int Ap[], const int Ai[], const std::complex<double> Ax[], 00053 void *Symbolic, void **Numeric, 00054 const double Control[UMFPACK_CONTROL],double Info [UMFPACK_INFO]) 00055 { 00056 return umfpack_zi_numeric(Ap,Ai,&numext::real_ref(Ax[0]),0,Symbolic,Numeric,Control,Info); 00057 } 00058 00059 inline int umfpack_solve( int sys, const int Ap[], const int Ai[], const double Ax[], 00060 double X[], const double B[], void *Numeric, 00061 const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO]) 00062 { 00063 return umfpack_di_solve(sys,Ap,Ai,Ax,X,B,Numeric,Control,Info); 00064 } 00065 00066 inline int umfpack_solve( int sys, const int Ap[], const int Ai[], const std::complex<double> Ax[], 00067 std::complex<double> X[], const std::complex<double> B[], void *Numeric, 00068 const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO]) 00069 { 00070 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); 00071 } 00072 00073 inline int umfpack_get_lunz(int *lnz, int *unz, int *n_row, int *n_col, int *nz_udiag, void *Numeric, double) 00074 { 00075 return umfpack_di_get_lunz(lnz,unz,n_row,n_col,nz_udiag,Numeric); 00076 } 00077 00078 inline int umfpack_get_lunz(int *lnz, int *unz, int *n_row, int *n_col, int *nz_udiag, void *Numeric, std::complex<double>) 00079 { 00080 return umfpack_zi_get_lunz(lnz,unz,n_row,n_col,nz_udiag,Numeric); 00081 } 00082 00083 inline int umfpack_get_numeric(int Lp[], int Lj[], double Lx[], int Up[], int Ui[], double Ux[], 00084 int P[], int Q[], double Dx[], int *do_recip, double Rs[], void *Numeric) 00085 { 00086 return umfpack_di_get_numeric(Lp,Lj,Lx,Up,Ui,Ux,P,Q,Dx,do_recip,Rs,Numeric); 00087 } 00088 00089 inline int umfpack_get_numeric(int Lp[], int Lj[], std::complex<double> Lx[], int Up[], int Ui[], std::complex<double> Ux[], 00090 int P[], int Q[], std::complex<double> Dx[], int *do_recip, double Rs[], void *Numeric) 00091 { 00092 double& lx0_real = numext::real_ref(Lx[0]); 00093 double& ux0_real = numext::real_ref(Ux[0]); 00094 double& dx0_real = numext::real_ref(Dx[0]); 00095 return umfpack_zi_get_numeric(Lp,Lj,Lx?&lx0_real:0,0,Up,Ui,Ux?&ux0_real:0,0,P,Q, 00096 Dx?&dx0_real:0,0,do_recip,Rs,Numeric); 00097 } 00098 00099 inline int umfpack_get_determinant(double *Mx, double *Ex, void *NumericHandle, double User_Info [UMFPACK_INFO]) 00100 { 00101 return umfpack_di_get_determinant(Mx,Ex,NumericHandle,User_Info); 00102 } 00103 00104 inline int umfpack_get_determinant(std::complex<double> *Mx, double *Ex, void *NumericHandle, double User_Info [UMFPACK_INFO]) 00105 { 00106 double& mx_real = numext::real_ref(*Mx); 00107 return umfpack_zi_get_determinant(&mx_real,0,Ex,NumericHandle,User_Info); 00108 } 00109 00110 namespace internal { 00111 template<typename T> struct umfpack_helper_is_sparse_plain : false_type {}; 00112 template<typename Scalar, int Options, typename StorageIndex> 00113 struct umfpack_helper_is_sparse_plain<SparseMatrix<Scalar,Options,StorageIndex> > 00114 : true_type {}; 00115 template<typename Scalar, int Options, typename StorageIndex> 00116 struct umfpack_helper_is_sparse_plain<MappedSparseMatrix<Scalar,Options,StorageIndex> > 00117 : true_type {}; 00118 } 00119 00133 template<typename _MatrixType> 00134 class UmfPackLU : internal::noncopyable 00135 { 00136 public: 00137 typedef _MatrixType MatrixType; 00138 typedef typename MatrixType::Scalar Scalar; 00139 typedef typename MatrixType::RealScalar RealScalar; 00140 typedef typename MatrixType::Index Index; 00141 typedef Matrix<Scalar,Dynamic,1> Vector; 00142 typedef Matrix<int, 1, MatrixType::ColsAtCompileTime> IntRowVectorType; 00143 typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType; 00144 typedef SparseMatrix<Scalar> LUMatrixType; 00145 typedef SparseMatrix<Scalar,ColMajor,int> UmfpackMatrixType; 00146 00147 public: 00148 00149 UmfPackLU() { init(); } 00150 00151 UmfPackLU(const MatrixType& matrix) 00152 { 00153 init(); 00154 compute(matrix); 00155 } 00156 00157 ~UmfPackLU() 00158 { 00159 if(m_symbolic) umfpack_free_symbolic(&m_symbolic,Scalar()); 00160 if(m_numeric) umfpack_free_numeric(&m_numeric,Scalar()); 00161 } 00162 00163 inline Index rows() const { return m_copyMatrix.rows(); } 00164 inline Index cols() const { return m_copyMatrix.cols(); } 00165 00171 ComputationInfo info() const 00172 { 00173 eigen_assert(m_isInitialized && "Decomposition is not initialized."); 00174 return m_info; 00175 } 00176 00177 inline const LUMatrixType& matrixL() const 00178 { 00179 if (m_extractedDataAreDirty) extractData(); 00180 return m_l; 00181 } 00182 00183 inline const LUMatrixType& matrixU() const 00184 { 00185 if (m_extractedDataAreDirty) extractData(); 00186 return m_u; 00187 } 00188 00189 inline const IntColVectorType& permutationP() const 00190 { 00191 if (m_extractedDataAreDirty) extractData(); 00192 return m_p; 00193 } 00194 00195 inline const IntRowVectorType& permutationQ() const 00196 { 00197 if (m_extractedDataAreDirty) extractData(); 00198 return m_q; 00199 } 00200 00205 template<typename InputMatrixType> 00206 void compute(const InputMatrixType& matrix) 00207 { 00208 if(m_symbolic) umfpack_free_symbolic(&m_symbolic,Scalar()); 00209 if(m_numeric) umfpack_free_numeric(&m_numeric,Scalar()); 00210 grapInput(matrix.derived()); 00211 analyzePattern_impl(); 00212 factorize_impl(); 00213 } 00214 00219 template<typename Rhs> 00220 inline const internal::solve_retval<UmfPackLU, Rhs> solve(const MatrixBase<Rhs>& b) const 00221 { 00222 eigen_assert(m_isInitialized && "UmfPackLU is not initialized."); 00223 eigen_assert(rows()==b.rows() 00224 && "UmfPackLU::solve(): invalid number of rows of the right hand side matrix b"); 00225 return internal::solve_retval<UmfPackLU, Rhs>(*this, b.derived()); 00226 } 00227 00232 template<typename Rhs> 00233 inline const internal::sparse_solve_retval<UmfPackLU, Rhs> solve(const SparseMatrixBase<Rhs>& b) const 00234 { 00235 eigen_assert(m_isInitialized && "UmfPackLU is not initialized."); 00236 eigen_assert(rows()==b.rows() 00237 && "UmfPackLU::solve(): invalid number of rows of the right hand side matrix b"); 00238 return internal::sparse_solve_retval<UmfPackLU, Rhs>(*this, b.derived()); 00239 } 00240 00247 template<typename InputMatrixType> 00248 void analyzePattern(const InputMatrixType& matrix) 00249 { 00250 if(m_symbolic) umfpack_free_symbolic(&m_symbolic,Scalar()); 00251 if(m_numeric) umfpack_free_numeric(&m_numeric,Scalar()); 00252 00253 grapInput(matrix.derived()); 00254 00255 analyzePattern_impl(); 00256 } 00257 00264 template<typename InputMatrixType> 00265 void factorize(const InputMatrixType& matrix) 00266 { 00267 eigen_assert(m_analysisIsOk && "UmfPackLU: you must first call analyzePattern()"); 00268 if(m_numeric) 00269 umfpack_free_numeric(&m_numeric,Scalar()); 00270 00271 grapInput(matrix.derived()); 00272 00273 factorize_impl(); 00274 } 00275 00276 #ifndef EIGEN_PARSED_BY_DOXYGEN 00277 00278 template<typename BDerived,typename XDerived> 00279 bool _solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived> &x) const; 00280 #endif 00281 00282 Scalar determinant() const; 00283 00284 void extractData() const; 00285 00286 protected: 00287 00288 void init() 00289 { 00290 m_info = InvalidInput; 00291 m_isInitialized = false; 00292 m_numeric = 0; 00293 m_symbolic = 0; 00294 m_outerIndexPtr = 0; 00295 m_innerIndexPtr = 0; 00296 m_valuePtr = 0; 00297 m_extractedDataAreDirty = true; 00298 } 00299 00300 template<typename InputMatrixType> 00301 void grapInput_impl(const InputMatrixType& mat, internal::true_type) 00302 { 00303 m_copyMatrix.resize(mat.rows(), mat.cols()); 00304 if( ((MatrixType::Flags&RowMajorBit)==RowMajorBit) || sizeof(typename MatrixType::Index)!=sizeof(int) || !mat.isCompressed() ) 00305 { 00306 // non supported input -> copy 00307 m_copyMatrix = mat; 00308 m_outerIndexPtr = m_copyMatrix.outerIndexPtr(); 00309 m_innerIndexPtr = m_copyMatrix.innerIndexPtr(); 00310 m_valuePtr = m_copyMatrix.valuePtr(); 00311 } 00312 else 00313 { 00314 m_outerIndexPtr = mat.outerIndexPtr(); 00315 m_innerIndexPtr = mat.innerIndexPtr(); 00316 m_valuePtr = mat.valuePtr(); 00317 } 00318 } 00319 00320 template<typename InputMatrixType> 00321 void grapInput_impl(const InputMatrixType& mat, internal::false_type) 00322 { 00323 m_copyMatrix = mat; 00324 m_outerIndexPtr = m_copyMatrix.outerIndexPtr(); 00325 m_innerIndexPtr = m_copyMatrix.innerIndexPtr(); 00326 m_valuePtr = m_copyMatrix.valuePtr(); 00327 } 00328 00329 template<typename InputMatrixType> 00330 void grapInput(const InputMatrixType& mat) 00331 { 00332 grapInput_impl(mat, internal::umfpack_helper_is_sparse_plain<InputMatrixType>()); 00333 } 00334 00335 void analyzePattern_impl() 00336 { 00337 int errorCode = 0; 00338 errorCode = umfpack_symbolic(m_copyMatrix.rows(), m_copyMatrix.cols(), m_outerIndexPtr, m_innerIndexPtr, m_valuePtr, 00339 &m_symbolic, 0, 0); 00340 00341 m_isInitialized = true; 00342 m_info = errorCode ? InvalidInput : Success; 00343 m_analysisIsOk = true; 00344 m_factorizationIsOk = false; 00345 m_extractedDataAreDirty = true; 00346 } 00347 00348 void factorize_impl() 00349 { 00350 int errorCode; 00351 errorCode = umfpack_numeric(m_outerIndexPtr, m_innerIndexPtr, m_valuePtr, 00352 m_symbolic, &m_numeric, 0, 0); 00353 00354 m_info = errorCode ? NumericalIssue : Success; 00355 m_factorizationIsOk = true; 00356 m_extractedDataAreDirty = true; 00357 } 00358 00359 // cached data to reduce reallocation, etc. 00360 mutable LUMatrixType m_l; 00361 mutable LUMatrixType m_u; 00362 mutable IntColVectorType m_p; 00363 mutable IntRowVectorType m_q; 00364 00365 UmfpackMatrixType m_copyMatrix; 00366 const Scalar* m_valuePtr; 00367 const int* m_outerIndexPtr; 00368 const int* m_innerIndexPtr; 00369 void* m_numeric; 00370 void* m_symbolic; 00371 00372 mutable ComputationInfo m_info; 00373 bool m_isInitialized; 00374 int m_factorizationIsOk; 00375 int m_analysisIsOk; 00376 mutable bool m_extractedDataAreDirty; 00377 00378 private: 00379 UmfPackLU(UmfPackLU& ) { } 00380 }; 00381 00382 00383 template<typename MatrixType> 00384 void UmfPackLU<MatrixType>::extractData() const 00385 { 00386 if (m_extractedDataAreDirty) 00387 { 00388 // get size of the data 00389 int lnz, unz, rows, cols, nz_udiag; 00390 umfpack_get_lunz(&lnz, &unz, &rows, &cols, &nz_udiag, m_numeric, Scalar()); 00391 00392 // allocate data 00393 m_l.resize(rows,(std::min)(rows,cols)); 00394 m_l.resizeNonZeros(lnz); 00395 00396 m_u.resize((std::min)(rows,cols),cols); 00397 m_u.resizeNonZeros(unz); 00398 00399 m_p.resize(rows); 00400 m_q.resize(cols); 00401 00402 // extract 00403 umfpack_get_numeric(m_l.outerIndexPtr(), m_l.innerIndexPtr(), m_l.valuePtr(), 00404 m_u.outerIndexPtr(), m_u.innerIndexPtr(), m_u.valuePtr(), 00405 m_p.data(), m_q.data(), 0, 0, 0, m_numeric); 00406 00407 m_extractedDataAreDirty = false; 00408 } 00409 } 00410 00411 template<typename MatrixType> 00412 typename UmfPackLU<MatrixType>::Scalar UmfPackLU<MatrixType>::determinant() const 00413 { 00414 Scalar det; 00415 umfpack_get_determinant(&det, 0, m_numeric, 0); 00416 return det; 00417 } 00418 00419 template<typename MatrixType> 00420 template<typename BDerived,typename XDerived> 00421 bool UmfPackLU<MatrixType>::_solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived> &x) const 00422 { 00423 const int rhsCols = b.cols(); 00424 eigen_assert((BDerived::Flags&RowMajorBit)==0 && "UmfPackLU backend does not support non col-major rhs yet"); 00425 eigen_assert((XDerived::Flags&RowMajorBit)==0 && "UmfPackLU backend does not support non col-major result yet"); 00426 eigen_assert(b.derived().data() != x.derived().data() && " Umfpack does not support inplace solve"); 00427 00428 int errorCode; 00429 for (int j=0; j<rhsCols; ++j) 00430 { 00431 errorCode = umfpack_solve(UMFPACK_A, 00432 m_outerIndexPtr, m_innerIndexPtr, m_valuePtr, 00433 &x.col(j).coeffRef(0), &b.const_cast_derived().col(j).coeffRef(0), m_numeric, 0, 0); 00434 if (errorCode!=0) 00435 return false; 00436 } 00437 00438 return true; 00439 } 00440 00441 00442 namespace internal { 00443 00444 template<typename _MatrixType, typename Rhs> 00445 struct solve_retval<UmfPackLU<_MatrixType>, Rhs> 00446 : solve_retval_base<UmfPackLU<_MatrixType>, Rhs> 00447 { 00448 typedef UmfPackLU<_MatrixType> Dec; 00449 EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) 00450 00451 template<typename Dest> void evalTo(Dest& dst) const 00452 { 00453 dec()._solve(rhs(),dst); 00454 } 00455 }; 00456 00457 template<typename _MatrixType, typename Rhs> 00458 struct sparse_solve_retval<UmfPackLU<_MatrixType>, Rhs> 00459 : sparse_solve_retval_base<UmfPackLU<_MatrixType>, Rhs> 00460 { 00461 typedef UmfPackLU<_MatrixType> Dec; 00462 EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs) 00463 00464 template<typename Dest> void evalTo(Dest& dst) const 00465 { 00466 this->defaultEvalTo(dst); 00467 } 00468 }; 00469 00470 } // end namespace internal 00471 00472 } // end namespace Eigen 00473 00474 #endif // EIGEN_UMFPACKSUPPORT_H