![]() |
Eigen
3.3.3
|
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 00020 inline void umfpack_defaults(double control[UMFPACK_CONTROL], double) 00021 { umfpack_di_defaults(control); } 00022 00023 inline void umfpack_defaults(double control[UMFPACK_CONTROL], std::complex<double>) 00024 { umfpack_zi_defaults(control); } 00025 00026 inline void umfpack_free_numeric(void **Numeric, double) 00027 { umfpack_di_free_numeric(Numeric); *Numeric = 0; } 00028 00029 inline void umfpack_free_numeric(void **Numeric, std::complex<double>) 00030 { umfpack_zi_free_numeric(Numeric); *Numeric = 0; } 00031 00032 inline void umfpack_free_symbolic(void **Symbolic, double) 00033 { umfpack_di_free_symbolic(Symbolic); *Symbolic = 0; } 00034 00035 inline void umfpack_free_symbolic(void **Symbolic, std::complex<double>) 00036 { umfpack_zi_free_symbolic(Symbolic); *Symbolic = 0; } 00037 00038 inline int umfpack_symbolic(int n_row,int n_col, 00039 const int Ap[], const int Ai[], const double Ax[], void **Symbolic, 00040 const double Control [UMFPACK_CONTROL], double Info [UMFPACK_INFO]) 00041 { 00042 return umfpack_di_symbolic(n_row,n_col,Ap,Ai,Ax,Symbolic,Control,Info); 00043 } 00044 00045 inline int umfpack_symbolic(int n_row,int n_col, 00046 const int Ap[], const int Ai[], const std::complex<double> Ax[], void **Symbolic, 00047 const double Control [UMFPACK_CONTROL], double Info [UMFPACK_INFO]) 00048 { 00049 return umfpack_zi_symbolic(n_row,n_col,Ap,Ai,&numext::real_ref(Ax[0]),0,Symbolic,Control,Info); 00050 } 00051 00052 inline int umfpack_numeric( const int Ap[], const int Ai[], const double Ax[], 00053 void *Symbolic, void **Numeric, 00054 const double Control[UMFPACK_CONTROL],double Info [UMFPACK_INFO]) 00055 { 00056 return umfpack_di_numeric(Ap,Ai,Ax,Symbolic,Numeric,Control,Info); 00057 } 00058 00059 inline int umfpack_numeric( const int Ap[], const int Ai[], const std::complex<double> Ax[], 00060 void *Symbolic, void **Numeric, 00061 const double Control[UMFPACK_CONTROL],double Info [UMFPACK_INFO]) 00062 { 00063 return umfpack_zi_numeric(Ap,Ai,&numext::real_ref(Ax[0]),0,Symbolic,Numeric,Control,Info); 00064 } 00065 00066 inline int umfpack_solve( int sys, const int Ap[], const int Ai[], const double Ax[], 00067 double X[], const double B[], void *Numeric, 00068 const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO]) 00069 { 00070 return umfpack_di_solve(sys,Ap,Ai,Ax,X,B,Numeric,Control,Info); 00071 } 00072 00073 inline int umfpack_solve( int sys, const int Ap[], const int Ai[], const std::complex<double> Ax[], 00074 std::complex<double> X[], const std::complex<double> B[], void *Numeric, 00075 const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO]) 00076 { 00077 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); 00078 } 00079 00080 inline int umfpack_get_lunz(int *lnz, int *unz, int *n_row, int *n_col, int *nz_udiag, void *Numeric, double) 00081 { 00082 return umfpack_di_get_lunz(lnz,unz,n_row,n_col,nz_udiag,Numeric); 00083 } 00084 00085 inline int umfpack_get_lunz(int *lnz, int *unz, int *n_row, int *n_col, int *nz_udiag, void *Numeric, std::complex<double>) 00086 { 00087 return umfpack_zi_get_lunz(lnz,unz,n_row,n_col,nz_udiag,Numeric); 00088 } 00089 00090 inline int umfpack_get_numeric(int Lp[], int Lj[], double Lx[], int Up[], int Ui[], double Ux[], 00091 int P[], int Q[], double Dx[], int *do_recip, double Rs[], void *Numeric) 00092 { 00093 return umfpack_di_get_numeric(Lp,Lj,Lx,Up,Ui,Ux,P,Q,Dx,do_recip,Rs,Numeric); 00094 } 00095 00096 inline int umfpack_get_numeric(int Lp[], int Lj[], std::complex<double> Lx[], int Up[], int Ui[], std::complex<double> Ux[], 00097 int P[], int Q[], std::complex<double> Dx[], int *do_recip, double Rs[], void *Numeric) 00098 { 00099 double& lx0_real = numext::real_ref(Lx[0]); 00100 double& ux0_real = numext::real_ref(Ux[0]); 00101 double& dx0_real = numext::real_ref(Dx[0]); 00102 return umfpack_zi_get_numeric(Lp,Lj,Lx?&lx0_real:0,0,Up,Ui,Ux?&ux0_real:0,0,P,Q, 00103 Dx?&dx0_real:0,0,do_recip,Rs,Numeric); 00104 } 00105 00106 inline int umfpack_get_determinant(double *Mx, double *Ex, void *NumericHandle, double User_Info [UMFPACK_INFO]) 00107 { 00108 return umfpack_di_get_determinant(Mx,Ex,NumericHandle,User_Info); 00109 } 00110 00111 inline int umfpack_get_determinant(std::complex<double> *Mx, double *Ex, void *NumericHandle, double User_Info [UMFPACK_INFO]) 00112 { 00113 double& mx_real = numext::real_ref(*Mx); 00114 return umfpack_zi_get_determinant(&mx_real,0,Ex,NumericHandle,User_Info); 00115 } 00116 00117 00133 template<typename _MatrixType> 00134 class UmfPackLU : public SparseSolverBase<UmfPackLU<_MatrixType> > 00135 { 00136 protected: 00137 typedef SparseSolverBase<UmfPackLU<_MatrixType> > Base; 00138 using Base::m_isInitialized; 00139 public: 00140 using Base::_solve_impl; 00141 typedef _MatrixType MatrixType; 00142 typedef typename MatrixType::Scalar Scalar; 00143 typedef typename MatrixType::RealScalar RealScalar; 00144 typedef typename MatrixType::StorageIndex StorageIndex; 00145 typedef Matrix<Scalar,Dynamic,1> Vector; 00146 typedef Matrix<int, 1, MatrixType::ColsAtCompileTime> IntRowVectorType; 00147 typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType; 00148 typedef SparseMatrix<Scalar> LUMatrixType; 00149 typedef SparseMatrix<Scalar,ColMajor,int> UmfpackMatrixType; 00150 typedef Ref<const UmfpackMatrixType, StandardCompressedFormat> UmfpackMatrixRef; 00151 enum { 00152 ColsAtCompileTime = MatrixType::ColsAtCompileTime, 00153 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime 00154 }; 00155 00156 public: 00157 00158 typedef Array<double, UMFPACK_CONTROL, 1> UmfpackControl; 00159 00160 UmfPackLU() 00161 : m_dummy(0,0), mp_matrix(m_dummy) 00162 { 00163 init(); 00164 } 00165 00166 template<typename InputMatrixType> 00167 explicit UmfPackLU(const InputMatrixType& matrix) 00168 : mp_matrix(matrix) 00169 { 00170 init(); 00171 compute(matrix); 00172 } 00173 00174 ~UmfPackLU() 00175 { 00176 if(m_symbolic) umfpack_free_symbolic(&m_symbolic,Scalar()); 00177 if(m_numeric) umfpack_free_numeric(&m_numeric,Scalar()); 00178 } 00179 00180 inline Index rows() const { return mp_matrix.rows(); } 00181 inline Index cols() const { return mp_matrix.cols(); } 00182 00188 ComputationInfo info() const 00189 { 00190 eigen_assert(m_isInitialized && "Decomposition is not initialized."); 00191 return m_info; 00192 } 00193 00194 inline const LUMatrixType& matrixL() const 00195 { 00196 if (m_extractedDataAreDirty) extractData(); 00197 return m_l; 00198 } 00199 00200 inline const LUMatrixType& matrixU() const 00201 { 00202 if (m_extractedDataAreDirty) extractData(); 00203 return m_u; 00204 } 00205 00206 inline const IntColVectorType& permutationP() const 00207 { 00208 if (m_extractedDataAreDirty) extractData(); 00209 return m_p; 00210 } 00211 00212 inline const IntRowVectorType& permutationQ() const 00213 { 00214 if (m_extractedDataAreDirty) extractData(); 00215 return m_q; 00216 } 00217 00222 template<typename InputMatrixType> 00223 void compute(const InputMatrixType& matrix) 00224 { 00225 if(m_symbolic) umfpack_free_symbolic(&m_symbolic,Scalar()); 00226 if(m_numeric) umfpack_free_numeric(&m_numeric,Scalar()); 00227 grab(matrix.derived()); 00228 analyzePattern_impl(); 00229 factorize_impl(); 00230 } 00231 00238 template<typename InputMatrixType> 00239 void analyzePattern(const InputMatrixType& matrix) 00240 { 00241 if(m_symbolic) umfpack_free_symbolic(&m_symbolic,Scalar()); 00242 if(m_numeric) umfpack_free_numeric(&m_numeric,Scalar()); 00243 00244 grab(matrix.derived()); 00245 00246 analyzePattern_impl(); 00247 } 00248 00254 inline int umfpackFactorizeReturncode() const 00255 { 00256 eigen_assert(m_numeric && "UmfPackLU: you must first call factorize()"); 00257 return m_fact_errorCode; 00258 } 00259 00266 inline const UmfpackControl& umfpackControl() const 00267 { 00268 return m_control; 00269 } 00270 00277 inline UmfpackControl& umfpackControl() 00278 { 00279 return m_control; 00280 } 00281 00288 template<typename InputMatrixType> 00289 void factorize(const InputMatrixType& matrix) 00290 { 00291 eigen_assert(m_analysisIsOk && "UmfPackLU: you must first call analyzePattern()"); 00292 if(m_numeric) 00293 umfpack_free_numeric(&m_numeric,Scalar()); 00294 00295 grab(matrix.derived()); 00296 00297 factorize_impl(); 00298 } 00299 00301 template<typename BDerived,typename XDerived> 00302 bool _solve_impl(const MatrixBase<BDerived> &b, MatrixBase<XDerived> &x) const; 00303 00304 Scalar determinant() const; 00305 00306 void extractData() const; 00307 00308 protected: 00309 00310 void init() 00311 { 00312 m_info = InvalidInput; 00313 m_isInitialized = false; 00314 m_numeric = 0; 00315 m_symbolic = 0; 00316 m_extractedDataAreDirty = true; 00317 } 00318 00319 void analyzePattern_impl() 00320 { 00321 umfpack_defaults(m_control.data(), Scalar()); 00322 int errorCode = 0; 00323 errorCode = umfpack_symbolic(internal::convert_index<int>(mp_matrix.rows()), 00324 internal::convert_index<int>(mp_matrix.cols()), 00325 mp_matrix.outerIndexPtr(), mp_matrix.innerIndexPtr(), mp_matrix.valuePtr(), 00326 &m_symbolic, m_control.data(), 0); 00327 00328 m_isInitialized = true; 00329 m_info = errorCode ? InvalidInput : Success; 00330 m_analysisIsOk = true; 00331 m_factorizationIsOk = false; 00332 m_extractedDataAreDirty = true; 00333 } 00334 00335 void factorize_impl() 00336 { 00337 m_fact_errorCode = umfpack_numeric(mp_matrix.outerIndexPtr(), mp_matrix.innerIndexPtr(), mp_matrix.valuePtr(), 00338 m_symbolic, &m_numeric, m_control.data(), 0); 00339 00340 m_info = m_fact_errorCode == UMFPACK_OK ? Success : NumericalIssue; 00341 m_factorizationIsOk = true; 00342 m_extractedDataAreDirty = true; 00343 } 00344 00345 template<typename MatrixDerived> 00346 void grab(const EigenBase<MatrixDerived> &A) 00347 { 00348 mp_matrix.~UmfpackMatrixRef(); 00349 ::new (&mp_matrix) UmfpackMatrixRef(A.derived()); 00350 } 00351 00352 void grab(const UmfpackMatrixRef &A) 00353 { 00354 if(&(A.derived()) != &mp_matrix) 00355 { 00356 mp_matrix.~UmfpackMatrixRef(); 00357 ::new (&mp_matrix) UmfpackMatrixRef(A); 00358 } 00359 } 00360 00361 // cached data to reduce reallocation, etc. 00362 mutable LUMatrixType m_l; 00363 int m_fact_errorCode; 00364 UmfpackControl m_control; 00365 00366 mutable LUMatrixType m_u; 00367 mutable IntColVectorType m_p; 00368 mutable IntRowVectorType m_q; 00369 00370 UmfpackMatrixType m_dummy; 00371 UmfpackMatrixRef mp_matrix; 00372 00373 void* m_numeric; 00374 void* m_symbolic; 00375 00376 mutable ComputationInfo m_info; 00377 int m_factorizationIsOk; 00378 int m_analysisIsOk; 00379 mutable bool m_extractedDataAreDirty; 00380 00381 private: 00382 UmfPackLU(const UmfPackLU& ) { } 00383 }; 00384 00385 00386 template<typename MatrixType> 00387 void UmfPackLU<MatrixType>::extractData() const 00388 { 00389 if (m_extractedDataAreDirty) 00390 { 00391 // get size of the data 00392 int lnz, unz, rows, cols, nz_udiag; 00393 umfpack_get_lunz(&lnz, &unz, &rows, &cols, &nz_udiag, m_numeric, Scalar()); 00394 00395 // allocate data 00396 m_l.resize(rows,(std::min)(rows,cols)); 00397 m_l.resizeNonZeros(lnz); 00398 00399 m_u.resize((std::min)(rows,cols),cols); 00400 m_u.resizeNonZeros(unz); 00401 00402 m_p.resize(rows); 00403 m_q.resize(cols); 00404 00405 // extract 00406 umfpack_get_numeric(m_l.outerIndexPtr(), m_l.innerIndexPtr(), m_l.valuePtr(), 00407 m_u.outerIndexPtr(), m_u.innerIndexPtr(), m_u.valuePtr(), 00408 m_p.data(), m_q.data(), 0, 0, 0, m_numeric); 00409 00410 m_extractedDataAreDirty = false; 00411 } 00412 } 00413 00414 template<typename MatrixType> 00415 typename UmfPackLU<MatrixType>::Scalar UmfPackLU<MatrixType>::determinant() const 00416 { 00417 Scalar det; 00418 umfpack_get_determinant(&det, 0, m_numeric, 0); 00419 return det; 00420 } 00421 00422 template<typename MatrixType> 00423 template<typename BDerived,typename XDerived> 00424 bool UmfPackLU<MatrixType>::_solve_impl(const MatrixBase<BDerived> &b, MatrixBase<XDerived> &x) const 00425 { 00426 Index rhsCols = b.cols(); 00427 eigen_assert((BDerived::Flags&RowMajorBit)==0 && "UmfPackLU backend does not support non col-major rhs yet"); 00428 eigen_assert((XDerived::Flags&RowMajorBit)==0 && "UmfPackLU backend does not support non col-major result yet"); 00429 eigen_assert(b.derived().data() != x.derived().data() && " Umfpack does not support inplace solve"); 00430 00431 int errorCode; 00432 Scalar* x_ptr = 0; 00433 Matrix<Scalar,Dynamic,1> x_tmp; 00434 if(x.innerStride()!=1) 00435 { 00436 x_tmp.resize(x.rows()); 00437 x_ptr = x_tmp.data(); 00438 } 00439 for (int j=0; j<rhsCols; ++j) 00440 { 00441 if(x.innerStride()==1) 00442 x_ptr = &x.col(j).coeffRef(0); 00443 errorCode = umfpack_solve(UMFPACK_A, 00444 mp_matrix.outerIndexPtr(), mp_matrix.innerIndexPtr(), mp_matrix.valuePtr(), 00445 x_ptr, &b.const_cast_derived().col(j).coeffRef(0), m_numeric, m_control.data(), 0); 00446 if(x.innerStride()!=1) 00447 x.col(j) = x_tmp; 00448 if (errorCode!=0) 00449 return false; 00450 } 00451 00452 return true; 00453 } 00454 00455 } // end namespace Eigen 00456 00457 #endif // EIGEN_UMFPACKSUPPORT_H