UmfPackSupport.h (24456B)
1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2008-2011 Gael Guennebaud <gael.guennebaud@inria.fr> 5 // 6 // This Source Code Form is subject to the terms of the Mozilla 7 // Public License v. 2.0. If a copy of the MPL was not distributed 8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9 10 #ifndef EIGEN_UMFPACKSUPPORT_H 11 #define EIGEN_UMFPACKSUPPORT_H 12 13 // for compatibility with super old version of umfpack, 14 // not sure this is really needed, but this is harmless. 15 #ifndef SuiteSparse_long 16 #ifdef UF_long 17 #define SuiteSparse_long UF_long 18 #else 19 #error neither SuiteSparse_long nor UF_long are defined 20 #endif 21 #endif 22 23 namespace Eigen { 24 25 /* TODO extract L, extract U, compute det, etc... */ 26 27 // generic double/complex<double> wrapper functions: 28 29 30 // Defaults 31 inline void umfpack_defaults(double control[UMFPACK_CONTROL], double, int) 32 { umfpack_di_defaults(control); } 33 34 inline void umfpack_defaults(double control[UMFPACK_CONTROL], std::complex<double>, int) 35 { umfpack_zi_defaults(control); } 36 37 inline void umfpack_defaults(double control[UMFPACK_CONTROL], double, SuiteSparse_long) 38 { umfpack_dl_defaults(control); } 39 40 inline void umfpack_defaults(double control[UMFPACK_CONTROL], std::complex<double>, SuiteSparse_long) 41 { umfpack_zl_defaults(control); } 42 43 // Report info 44 inline void umfpack_report_info(double control[UMFPACK_CONTROL], double info[UMFPACK_INFO], double, int) 45 { umfpack_di_report_info(control, info);} 46 47 inline void umfpack_report_info(double control[UMFPACK_CONTROL], double info[UMFPACK_INFO], std::complex<double>, int) 48 { umfpack_zi_report_info(control, info);} 49 50 inline void umfpack_report_info(double control[UMFPACK_CONTROL], double info[UMFPACK_INFO], double, SuiteSparse_long) 51 { umfpack_dl_report_info(control, info);} 52 53 inline void umfpack_report_info(double control[UMFPACK_CONTROL], double info[UMFPACK_INFO], std::complex<double>, SuiteSparse_long) 54 { umfpack_zl_report_info(control, info);} 55 56 // Report status 57 inline void umfpack_report_status(double control[UMFPACK_CONTROL], int status, double, int) 58 { umfpack_di_report_status(control, status);} 59 60 inline void umfpack_report_status(double control[UMFPACK_CONTROL], int status, std::complex<double>, int) 61 { umfpack_zi_report_status(control, status);} 62 63 inline void umfpack_report_status(double control[UMFPACK_CONTROL], int status, double, SuiteSparse_long) 64 { umfpack_dl_report_status(control, status);} 65 66 inline void umfpack_report_status(double control[UMFPACK_CONTROL], int status, std::complex<double>, SuiteSparse_long) 67 { umfpack_zl_report_status(control, status);} 68 69 // report control 70 inline void umfpack_report_control(double control[UMFPACK_CONTROL], double, int) 71 { umfpack_di_report_control(control);} 72 73 inline void umfpack_report_control(double control[UMFPACK_CONTROL], std::complex<double>, int) 74 { umfpack_zi_report_control(control);} 75 76 inline void umfpack_report_control(double control[UMFPACK_CONTROL], double, SuiteSparse_long) 77 { umfpack_dl_report_control(control);} 78 79 inline void umfpack_report_control(double control[UMFPACK_CONTROL], std::complex<double>, SuiteSparse_long) 80 { umfpack_zl_report_control(control);} 81 82 // Free numeric 83 inline void umfpack_free_numeric(void **Numeric, double, int) 84 { umfpack_di_free_numeric(Numeric); *Numeric = 0; } 85 86 inline void umfpack_free_numeric(void **Numeric, std::complex<double>, int) 87 { umfpack_zi_free_numeric(Numeric); *Numeric = 0; } 88 89 inline void umfpack_free_numeric(void **Numeric, double, SuiteSparse_long) 90 { umfpack_dl_free_numeric(Numeric); *Numeric = 0; } 91 92 inline void umfpack_free_numeric(void **Numeric, std::complex<double>, SuiteSparse_long) 93 { umfpack_zl_free_numeric(Numeric); *Numeric = 0; } 94 95 // Free symbolic 96 inline void umfpack_free_symbolic(void **Symbolic, double, int) 97 { umfpack_di_free_symbolic(Symbolic); *Symbolic = 0; } 98 99 inline void umfpack_free_symbolic(void **Symbolic, std::complex<double>, int) 100 { umfpack_zi_free_symbolic(Symbolic); *Symbolic = 0; } 101 102 inline void umfpack_free_symbolic(void **Symbolic, double, SuiteSparse_long) 103 { umfpack_dl_free_symbolic(Symbolic); *Symbolic = 0; } 104 105 inline void umfpack_free_symbolic(void **Symbolic, std::complex<double>, SuiteSparse_long) 106 { umfpack_zl_free_symbolic(Symbolic); *Symbolic = 0; } 107 108 // Symbolic 109 inline int umfpack_symbolic(int n_row,int n_col, 110 const int Ap[], const int Ai[], const double Ax[], void **Symbolic, 111 const double Control [UMFPACK_CONTROL], double Info [UMFPACK_INFO]) 112 { 113 return umfpack_di_symbolic(n_row,n_col,Ap,Ai,Ax,Symbolic,Control,Info); 114 } 115 116 inline int umfpack_symbolic(int n_row,int n_col, 117 const int Ap[], const int Ai[], const std::complex<double> Ax[], void **Symbolic, 118 const double Control [UMFPACK_CONTROL], double Info [UMFPACK_INFO]) 119 { 120 return umfpack_zi_symbolic(n_row,n_col,Ap,Ai,&numext::real_ref(Ax[0]),0,Symbolic,Control,Info); 121 } 122 inline SuiteSparse_long umfpack_symbolic( SuiteSparse_long n_row,SuiteSparse_long n_col, 123 const SuiteSparse_long Ap[], const SuiteSparse_long Ai[], const double Ax[], void **Symbolic, 124 const double Control [UMFPACK_CONTROL], double Info [UMFPACK_INFO]) 125 { 126 return umfpack_dl_symbolic(n_row,n_col,Ap,Ai,Ax,Symbolic,Control,Info); 127 } 128 129 inline SuiteSparse_long umfpack_symbolic( SuiteSparse_long n_row,SuiteSparse_long n_col, 130 const SuiteSparse_long Ap[], const SuiteSparse_long Ai[], const std::complex<double> Ax[], void **Symbolic, 131 const double Control [UMFPACK_CONTROL], double Info [UMFPACK_INFO]) 132 { 133 return umfpack_zl_symbolic(n_row,n_col,Ap,Ai,&numext::real_ref(Ax[0]),0,Symbolic,Control,Info); 134 } 135 136 // Numeric 137 inline int umfpack_numeric( const int Ap[], const int Ai[], const double Ax[], 138 void *Symbolic, void **Numeric, 139 const double Control[UMFPACK_CONTROL],double Info [UMFPACK_INFO]) 140 { 141 return umfpack_di_numeric(Ap,Ai,Ax,Symbolic,Numeric,Control,Info); 142 } 143 144 inline int umfpack_numeric( const int Ap[], const int Ai[], const std::complex<double> Ax[], 145 void *Symbolic, void **Numeric, 146 const double Control[UMFPACK_CONTROL],double Info [UMFPACK_INFO]) 147 { 148 return umfpack_zi_numeric(Ap,Ai,&numext::real_ref(Ax[0]),0,Symbolic,Numeric,Control,Info); 149 } 150 inline SuiteSparse_long umfpack_numeric(const SuiteSparse_long Ap[], const SuiteSparse_long Ai[], const double Ax[], 151 void *Symbolic, void **Numeric, 152 const double Control[UMFPACK_CONTROL],double Info [UMFPACK_INFO]) 153 { 154 return umfpack_dl_numeric(Ap,Ai,Ax,Symbolic,Numeric,Control,Info); 155 } 156 157 inline SuiteSparse_long umfpack_numeric(const SuiteSparse_long Ap[], const SuiteSparse_long Ai[], const std::complex<double> Ax[], 158 void *Symbolic, void **Numeric, 159 const double Control[UMFPACK_CONTROL],double Info [UMFPACK_INFO]) 160 { 161 return umfpack_zl_numeric(Ap,Ai,&numext::real_ref(Ax[0]),0,Symbolic,Numeric,Control,Info); 162 } 163 164 // solve 165 inline int umfpack_solve( int sys, const int Ap[], const int Ai[], const double Ax[], 166 double X[], const double B[], void *Numeric, 167 const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO]) 168 { 169 return umfpack_di_solve(sys,Ap,Ai,Ax,X,B,Numeric,Control,Info); 170 } 171 172 inline int umfpack_solve( int sys, const int Ap[], const int Ai[], const std::complex<double> Ax[], 173 std::complex<double> X[], const std::complex<double> B[], void *Numeric, 174 const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO]) 175 { 176 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); 177 } 178 179 inline SuiteSparse_long umfpack_solve(int sys, const SuiteSparse_long Ap[], const SuiteSparse_long Ai[], const double Ax[], 180 double X[], const double B[], void *Numeric, 181 const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO]) 182 { 183 return umfpack_dl_solve(sys,Ap,Ai,Ax,X,B,Numeric,Control,Info); 184 } 185 186 inline SuiteSparse_long umfpack_solve(int sys, const SuiteSparse_long Ap[], const SuiteSparse_long Ai[], const std::complex<double> Ax[], 187 std::complex<double> X[], const std::complex<double> B[], void *Numeric, 188 const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO]) 189 { 190 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); 191 } 192 193 // Get Lunz 194 inline int umfpack_get_lunz(int *lnz, int *unz, int *n_row, int *n_col, int *nz_udiag, void *Numeric, double) 195 { 196 return umfpack_di_get_lunz(lnz,unz,n_row,n_col,nz_udiag,Numeric); 197 } 198 199 inline int umfpack_get_lunz(int *lnz, int *unz, int *n_row, int *n_col, int *nz_udiag, void *Numeric, std::complex<double>) 200 { 201 return umfpack_zi_get_lunz(lnz,unz,n_row,n_col,nz_udiag,Numeric); 202 } 203 204 inline SuiteSparse_long umfpack_get_lunz( SuiteSparse_long *lnz, SuiteSparse_long *unz, SuiteSparse_long *n_row, SuiteSparse_long *n_col, 205 SuiteSparse_long *nz_udiag, void *Numeric, double) 206 { 207 return umfpack_dl_get_lunz(lnz,unz,n_row,n_col,nz_udiag,Numeric); 208 } 209 210 inline SuiteSparse_long umfpack_get_lunz( SuiteSparse_long *lnz, SuiteSparse_long *unz, SuiteSparse_long *n_row, SuiteSparse_long *n_col, 211 SuiteSparse_long *nz_udiag, void *Numeric, std::complex<double>) 212 { 213 return umfpack_zl_get_lunz(lnz,unz,n_row,n_col,nz_udiag,Numeric); 214 } 215 216 // Get Numeric 217 inline int umfpack_get_numeric(int Lp[], int Lj[], double Lx[], int Up[], int Ui[], double Ux[], 218 int P[], int Q[], double Dx[], int *do_recip, double Rs[], void *Numeric) 219 { 220 return umfpack_di_get_numeric(Lp,Lj,Lx,Up,Ui,Ux,P,Q,Dx,do_recip,Rs,Numeric); 221 } 222 223 inline int umfpack_get_numeric(int Lp[], int Lj[], std::complex<double> Lx[], int Up[], int Ui[], std::complex<double> Ux[], 224 int P[], int Q[], std::complex<double> Dx[], int *do_recip, double Rs[], void *Numeric) 225 { 226 double& lx0_real = numext::real_ref(Lx[0]); 227 double& ux0_real = numext::real_ref(Ux[0]); 228 double& dx0_real = numext::real_ref(Dx[0]); 229 return umfpack_zi_get_numeric(Lp,Lj,Lx?&lx0_real:0,0,Up,Ui,Ux?&ux0_real:0,0,P,Q, 230 Dx?&dx0_real:0,0,do_recip,Rs,Numeric); 231 } 232 inline SuiteSparse_long umfpack_get_numeric(SuiteSparse_long Lp[], SuiteSparse_long Lj[], double Lx[], SuiteSparse_long Up[], SuiteSparse_long Ui[], double Ux[], 233 SuiteSparse_long P[], SuiteSparse_long Q[], double Dx[], SuiteSparse_long *do_recip, double Rs[], void *Numeric) 234 { 235 return umfpack_dl_get_numeric(Lp,Lj,Lx,Up,Ui,Ux,P,Q,Dx,do_recip,Rs,Numeric); 236 } 237 238 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[], 239 SuiteSparse_long P[], SuiteSparse_long Q[], std::complex<double> Dx[], SuiteSparse_long *do_recip, double Rs[], void *Numeric) 240 { 241 double& lx0_real = numext::real_ref(Lx[0]); 242 double& ux0_real = numext::real_ref(Ux[0]); 243 double& dx0_real = numext::real_ref(Dx[0]); 244 return umfpack_zl_get_numeric(Lp,Lj,Lx?&lx0_real:0,0,Up,Ui,Ux?&ux0_real:0,0,P,Q, 245 Dx?&dx0_real:0,0,do_recip,Rs,Numeric); 246 } 247 248 // Get Determinant 249 inline int umfpack_get_determinant(double *Mx, double *Ex, void *NumericHandle, double User_Info [UMFPACK_INFO], int) 250 { 251 return umfpack_di_get_determinant(Mx,Ex,NumericHandle,User_Info); 252 } 253 254 inline int umfpack_get_determinant(std::complex<double> *Mx, double *Ex, void *NumericHandle, double User_Info [UMFPACK_INFO], int) 255 { 256 double& mx_real = numext::real_ref(*Mx); 257 return umfpack_zi_get_determinant(&mx_real,0,Ex,NumericHandle,User_Info); 258 } 259 260 inline SuiteSparse_long umfpack_get_determinant(double *Mx, double *Ex, void *NumericHandle, double User_Info [UMFPACK_INFO], SuiteSparse_long) 261 { 262 return umfpack_dl_get_determinant(Mx,Ex,NumericHandle,User_Info); 263 } 264 265 inline SuiteSparse_long umfpack_get_determinant(std::complex<double> *Mx, double *Ex, void *NumericHandle, double User_Info [UMFPACK_INFO], SuiteSparse_long) 266 { 267 double& mx_real = numext::real_ref(*Mx); 268 return umfpack_zl_get_determinant(&mx_real,0,Ex,NumericHandle,User_Info); 269 } 270 271 272 /** \ingroup UmfPackSupport_Module 273 * \brief A sparse LU factorization and solver based on UmfPack 274 * 275 * This class allows to solve for A.X = B sparse linear problems via a LU factorization 276 * using the UmfPack library. The sparse matrix A must be squared and full rank. 277 * The vectors or matrices X and B can be either dense or sparse. 278 * 279 * \warning The input matrix A should be in a \b compressed and \b column-major form. 280 * Otherwise an expensive copy will be made. You can call the inexpensive makeCompressed() to get a compressed matrix. 281 * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> 282 * 283 * \implsparsesolverconcept 284 * 285 * \sa \ref TutorialSparseSolverConcept, class SparseLU 286 */ 287 template<typename _MatrixType> 288 class UmfPackLU : public SparseSolverBase<UmfPackLU<_MatrixType> > 289 { 290 protected: 291 typedef SparseSolverBase<UmfPackLU<_MatrixType> > Base; 292 using Base::m_isInitialized; 293 public: 294 using Base::_solve_impl; 295 typedef _MatrixType MatrixType; 296 typedef typename MatrixType::Scalar Scalar; 297 typedef typename MatrixType::RealScalar RealScalar; 298 typedef typename MatrixType::StorageIndex StorageIndex; 299 typedef Matrix<Scalar,Dynamic,1> Vector; 300 typedef Matrix<int, 1, MatrixType::ColsAtCompileTime> IntRowVectorType; 301 typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType; 302 typedef SparseMatrix<Scalar> LUMatrixType; 303 typedef SparseMatrix<Scalar,ColMajor,StorageIndex> UmfpackMatrixType; 304 typedef Ref<const UmfpackMatrixType, StandardCompressedFormat> UmfpackMatrixRef; 305 enum { 306 ColsAtCompileTime = MatrixType::ColsAtCompileTime, 307 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime 308 }; 309 310 public: 311 312 typedef Array<double, UMFPACK_CONTROL, 1> UmfpackControl; 313 typedef Array<double, UMFPACK_INFO, 1> UmfpackInfo; 314 315 UmfPackLU() 316 : m_dummy(0,0), mp_matrix(m_dummy) 317 { 318 init(); 319 } 320 321 template<typename InputMatrixType> 322 explicit UmfPackLU(const InputMatrixType& matrix) 323 : mp_matrix(matrix) 324 { 325 init(); 326 compute(matrix); 327 } 328 329 ~UmfPackLU() 330 { 331 if(m_symbolic) umfpack_free_symbolic(&m_symbolic,Scalar(), StorageIndex()); 332 if(m_numeric) umfpack_free_numeric(&m_numeric,Scalar(), StorageIndex()); 333 } 334 335 inline Index rows() const { return mp_matrix.rows(); } 336 inline Index cols() const { return mp_matrix.cols(); } 337 338 /** \brief Reports whether previous computation was successful. 339 * 340 * \returns \c Success if computation was successful, 341 * \c NumericalIssue if the matrix.appears to be negative. 342 */ 343 ComputationInfo info() const 344 { 345 eigen_assert(m_isInitialized && "Decomposition is not initialized."); 346 return m_info; 347 } 348 349 inline const LUMatrixType& matrixL() const 350 { 351 if (m_extractedDataAreDirty) extractData(); 352 return m_l; 353 } 354 355 inline const LUMatrixType& matrixU() const 356 { 357 if (m_extractedDataAreDirty) extractData(); 358 return m_u; 359 } 360 361 inline const IntColVectorType& permutationP() const 362 { 363 if (m_extractedDataAreDirty) extractData(); 364 return m_p; 365 } 366 367 inline const IntRowVectorType& permutationQ() const 368 { 369 if (m_extractedDataAreDirty) extractData(); 370 return m_q; 371 } 372 373 /** Computes the sparse Cholesky decomposition of \a matrix 374 * Note that the matrix should be column-major, and in compressed format for best performance. 375 * \sa SparseMatrix::makeCompressed(). 376 */ 377 template<typename InputMatrixType> 378 void compute(const InputMatrixType& matrix) 379 { 380 if(m_symbolic) umfpack_free_symbolic(&m_symbolic,Scalar(),StorageIndex()); 381 if(m_numeric) umfpack_free_numeric(&m_numeric,Scalar(),StorageIndex()); 382 grab(matrix.derived()); 383 analyzePattern_impl(); 384 factorize_impl(); 385 } 386 387 /** Performs a symbolic decomposition on the sparcity of \a matrix. 388 * 389 * This function is particularly useful when solving for several problems having the same structure. 390 * 391 * \sa factorize(), compute() 392 */ 393 template<typename InputMatrixType> 394 void analyzePattern(const InputMatrixType& matrix) 395 { 396 if(m_symbolic) umfpack_free_symbolic(&m_symbolic,Scalar(),StorageIndex()); 397 if(m_numeric) umfpack_free_numeric(&m_numeric,Scalar(),StorageIndex()); 398 399 grab(matrix.derived()); 400 401 analyzePattern_impl(); 402 } 403 404 /** Provides the return status code returned by UmfPack during the numeric 405 * factorization. 406 * 407 * \sa factorize(), compute() 408 */ 409 inline int umfpackFactorizeReturncode() const 410 { 411 eigen_assert(m_numeric && "UmfPackLU: you must first call factorize()"); 412 return m_fact_errorCode; 413 } 414 415 /** Provides access to the control settings array used by UmfPack. 416 * 417 * If this array contains NaN's, the default values are used. 418 * 419 * See UMFPACK documentation for details. 420 */ 421 inline const UmfpackControl& umfpackControl() const 422 { 423 return m_control; 424 } 425 426 /** Provides access to the control settings array used by UmfPack. 427 * 428 * If this array contains NaN's, the default values are used. 429 * 430 * See UMFPACK documentation for details. 431 */ 432 inline UmfpackControl& umfpackControl() 433 { 434 return m_control; 435 } 436 437 /** Performs a numeric decomposition of \a matrix 438 * 439 * The given matrix must has the same sparcity than the matrix on which the pattern anylysis has been performed. 440 * 441 * \sa analyzePattern(), compute() 442 */ 443 template<typename InputMatrixType> 444 void factorize(const InputMatrixType& matrix) 445 { 446 eigen_assert(m_analysisIsOk && "UmfPackLU: you must first call analyzePattern()"); 447 if(m_numeric) 448 umfpack_free_numeric(&m_numeric,Scalar(),StorageIndex()); 449 450 grab(matrix.derived()); 451 452 factorize_impl(); 453 } 454 455 /** Prints the current UmfPack control settings. 456 * 457 * \sa umfpackControl() 458 */ 459 void printUmfpackControl() 460 { 461 umfpack_report_control(m_control.data(), Scalar(),StorageIndex()); 462 } 463 464 /** Prints statistics collected by UmfPack. 465 * 466 * \sa analyzePattern(), compute() 467 */ 468 void printUmfpackInfo() 469 { 470 eigen_assert(m_analysisIsOk && "UmfPackLU: you must first call analyzePattern()"); 471 umfpack_report_info(m_control.data(), m_umfpackInfo.data(), Scalar(),StorageIndex()); 472 } 473 474 /** Prints the status of the previous factorization operation performed by UmfPack (symbolic or numerical factorization). 475 * 476 * \sa analyzePattern(), compute() 477 */ 478 void printUmfpackStatus() { 479 eigen_assert(m_analysisIsOk && "UmfPackLU: you must first call analyzePattern()"); 480 umfpack_report_status(m_control.data(), m_fact_errorCode, Scalar(),StorageIndex()); 481 } 482 483 /** \internal */ 484 template<typename BDerived,typename XDerived> 485 bool _solve_impl(const MatrixBase<BDerived> &b, MatrixBase<XDerived> &x) const; 486 487 Scalar determinant() const; 488 489 void extractData() const; 490 491 protected: 492 493 void init() 494 { 495 m_info = InvalidInput; 496 m_isInitialized = false; 497 m_numeric = 0; 498 m_symbolic = 0; 499 m_extractedDataAreDirty = true; 500 501 umfpack_defaults(m_control.data(), Scalar(),StorageIndex()); 502 } 503 504 void analyzePattern_impl() 505 { 506 m_fact_errorCode = umfpack_symbolic(internal::convert_index<StorageIndex>(mp_matrix.rows()), 507 internal::convert_index<StorageIndex>(mp_matrix.cols()), 508 mp_matrix.outerIndexPtr(), mp_matrix.innerIndexPtr(), mp_matrix.valuePtr(), 509 &m_symbolic, m_control.data(), m_umfpackInfo.data()); 510 511 m_isInitialized = true; 512 m_info = m_fact_errorCode ? InvalidInput : Success; 513 m_analysisIsOk = true; 514 m_factorizationIsOk = false; 515 m_extractedDataAreDirty = true; 516 } 517 518 void factorize_impl() 519 { 520 521 m_fact_errorCode = umfpack_numeric(mp_matrix.outerIndexPtr(), mp_matrix.innerIndexPtr(), mp_matrix.valuePtr(), 522 m_symbolic, &m_numeric, m_control.data(), m_umfpackInfo.data()); 523 524 m_info = m_fact_errorCode == UMFPACK_OK ? Success : NumericalIssue; 525 m_factorizationIsOk = true; 526 m_extractedDataAreDirty = true; 527 } 528 529 template<typename MatrixDerived> 530 void grab(const EigenBase<MatrixDerived> &A) 531 { 532 mp_matrix.~UmfpackMatrixRef(); 533 ::new (&mp_matrix) UmfpackMatrixRef(A.derived()); 534 } 535 536 void grab(const UmfpackMatrixRef &A) 537 { 538 if(&(A.derived()) != &mp_matrix) 539 { 540 mp_matrix.~UmfpackMatrixRef(); 541 ::new (&mp_matrix) UmfpackMatrixRef(A); 542 } 543 } 544 545 // cached data to reduce reallocation, etc. 546 mutable LUMatrixType m_l; 547 StorageIndex m_fact_errorCode; 548 UmfpackControl m_control; 549 mutable UmfpackInfo m_umfpackInfo; 550 551 mutable LUMatrixType m_u; 552 mutable IntColVectorType m_p; 553 mutable IntRowVectorType m_q; 554 555 UmfpackMatrixType m_dummy; 556 UmfpackMatrixRef mp_matrix; 557 558 void* m_numeric; 559 void* m_symbolic; 560 561 mutable ComputationInfo m_info; 562 int m_factorizationIsOk; 563 int m_analysisIsOk; 564 mutable bool m_extractedDataAreDirty; 565 566 private: 567 UmfPackLU(const UmfPackLU& ) { } 568 }; 569 570 571 template<typename MatrixType> 572 void UmfPackLU<MatrixType>::extractData() const 573 { 574 if (m_extractedDataAreDirty) 575 { 576 // get size of the data 577 StorageIndex lnz, unz, rows, cols, nz_udiag; 578 umfpack_get_lunz(&lnz, &unz, &rows, &cols, &nz_udiag, m_numeric, Scalar()); 579 580 // allocate data 581 m_l.resize(rows,(std::min)(rows,cols)); 582 m_l.resizeNonZeros(lnz); 583 584 m_u.resize((std::min)(rows,cols),cols); 585 m_u.resizeNonZeros(unz); 586 587 m_p.resize(rows); 588 m_q.resize(cols); 589 590 // extract 591 umfpack_get_numeric(m_l.outerIndexPtr(), m_l.innerIndexPtr(), m_l.valuePtr(), 592 m_u.outerIndexPtr(), m_u.innerIndexPtr(), m_u.valuePtr(), 593 m_p.data(), m_q.data(), 0, 0, 0, m_numeric); 594 595 m_extractedDataAreDirty = false; 596 } 597 } 598 599 template<typename MatrixType> 600 typename UmfPackLU<MatrixType>::Scalar UmfPackLU<MatrixType>::determinant() const 601 { 602 Scalar det; 603 umfpack_get_determinant(&det, 0, m_numeric, 0, StorageIndex()); 604 return det; 605 } 606 607 template<typename MatrixType> 608 template<typename BDerived,typename XDerived> 609 bool UmfPackLU<MatrixType>::_solve_impl(const MatrixBase<BDerived> &b, MatrixBase<XDerived> &x) const 610 { 611 Index rhsCols = b.cols(); 612 eigen_assert((BDerived::Flags&RowMajorBit)==0 && "UmfPackLU backend does not support non col-major rhs yet"); 613 eigen_assert((XDerived::Flags&RowMajorBit)==0 && "UmfPackLU backend does not support non col-major result yet"); 614 eigen_assert(b.derived().data() != x.derived().data() && " Umfpack does not support inplace solve"); 615 616 Scalar* x_ptr = 0; 617 Matrix<Scalar,Dynamic,1> x_tmp; 618 if(x.innerStride()!=1) 619 { 620 x_tmp.resize(x.rows()); 621 x_ptr = x_tmp.data(); 622 } 623 for (int j=0; j<rhsCols; ++j) 624 { 625 if(x.innerStride()==1) 626 x_ptr = &x.col(j).coeffRef(0); 627 StorageIndex errorCode = umfpack_solve(UMFPACK_A, 628 mp_matrix.outerIndexPtr(), mp_matrix.innerIndexPtr(), mp_matrix.valuePtr(), 629 x_ptr, &b.const_cast_derived().col(j).coeffRef(0), 630 m_numeric, m_control.data(), m_umfpackInfo.data()); 631 if(x.innerStride()!=1) 632 x.col(j) = x_tmp; 633 if (errorCode!=0) 634 return false; 635 } 636 637 return true; 638 } 639 640 } // end namespace Eigen 641 642 #endif // EIGEN_UMFPACKSUPPORT_H