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 (at) 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 namespace Eigen { 14 15 /* TODO extract L, extract U, compute det, etc... */ 16 17 // generic double/complex<double> wrapper functions: 18 19 20 inline void umfpack_defaults(double control[UMFPACK_CONTROL], double) 21 { umfpack_di_defaults(control); } 22 23 inline void umfpack_defaults(double control[UMFPACK_CONTROL], std::complex<double>) 24 { umfpack_zi_defaults(control); } 25 26 inline void umfpack_report_info(double control[UMFPACK_CONTROL], double info[UMFPACK_INFO], double) 27 { umfpack_di_report_info(control, info);} 28 29 inline void umfpack_report_info(double control[UMFPACK_CONTROL], double info[UMFPACK_INFO], std::complex<double>) 30 { umfpack_zi_report_info(control, info);} 31 32 inline void umfpack_report_status(double control[UMFPACK_CONTROL], int status, double) 33 { umfpack_di_report_status(control, status);} 34 35 inline void umfpack_report_status(double control[UMFPACK_CONTROL], int status, std::complex<double>) 36 { umfpack_zi_report_status(control, status);} 37 38 inline void umfpack_report_control(double control[UMFPACK_CONTROL], double) 39 { umfpack_di_report_control(control);} 40 41 inline void umfpack_report_control(double control[UMFPACK_CONTROL], std::complex<double>) 42 { umfpack_zi_report_control(control);} 43 44 inline void umfpack_free_numeric(void **Numeric, double) 45 { umfpack_di_free_numeric(Numeric); *Numeric = 0; } 46 47 inline void umfpack_free_numeric(void **Numeric, std::complex<double>) 48 { umfpack_zi_free_numeric(Numeric); *Numeric = 0; } 49 50 inline void umfpack_free_symbolic(void **Symbolic, double) 51 { umfpack_di_free_symbolic(Symbolic); *Symbolic = 0; } 52 53 inline void umfpack_free_symbolic(void **Symbolic, std::complex<double>) 54 { umfpack_zi_free_symbolic(Symbolic); *Symbolic = 0; } 55 56 inline int umfpack_symbolic(int n_row,int n_col, 57 const int Ap[], const int Ai[], const double Ax[], void **Symbolic, 58 const double Control [UMFPACK_CONTROL], double Info [UMFPACK_INFO]) 59 { 60 return umfpack_di_symbolic(n_row,n_col,Ap,Ai,Ax,Symbolic,Control,Info); 61 } 62 63 inline int umfpack_symbolic(int n_row,int n_col, 64 const int Ap[], const int Ai[], const std::complex<double> Ax[], void **Symbolic, 65 const double Control [UMFPACK_CONTROL], double Info [UMFPACK_INFO]) 66 { 67 return umfpack_zi_symbolic(n_row,n_col,Ap,Ai,&numext::real_ref(Ax[0]),0,Symbolic,Control,Info); 68 } 69 70 inline int umfpack_numeric( const int Ap[], const int Ai[], const double Ax[], 71 void *Symbolic, void **Numeric, 72 const double Control[UMFPACK_CONTROL],double Info [UMFPACK_INFO]) 73 { 74 return umfpack_di_numeric(Ap,Ai,Ax,Symbolic,Numeric,Control,Info); 75 } 76 77 inline int umfpack_numeric( const int Ap[], const int Ai[], const std::complex<double> Ax[], 78 void *Symbolic, void **Numeric, 79 const double Control[UMFPACK_CONTROL],double Info [UMFPACK_INFO]) 80 { 81 return umfpack_zi_numeric(Ap,Ai,&numext::real_ref(Ax[0]),0,Symbolic,Numeric,Control,Info); 82 } 83 84 inline int umfpack_solve( int sys, const int Ap[], const int Ai[], const double Ax[], 85 double X[], const double B[], void *Numeric, 86 const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO]) 87 { 88 return umfpack_di_solve(sys,Ap,Ai,Ax,X,B,Numeric,Control,Info); 89 } 90 91 inline int umfpack_solve( int sys, const int Ap[], const int Ai[], const std::complex<double> Ax[], 92 std::complex<double> X[], const std::complex<double> B[], void *Numeric, 93 const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO]) 94 { 95 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); 96 } 97 98 inline int umfpack_get_lunz(int *lnz, int *unz, int *n_row, int *n_col, int *nz_udiag, void *Numeric, double) 99 { 100 return umfpack_di_get_lunz(lnz,unz,n_row,n_col,nz_udiag,Numeric); 101 } 102 103 inline int umfpack_get_lunz(int *lnz, int *unz, int *n_row, int *n_col, int *nz_udiag, void *Numeric, std::complex<double>) 104 { 105 return umfpack_zi_get_lunz(lnz,unz,n_row,n_col,nz_udiag,Numeric); 106 } 107 108 inline int umfpack_get_numeric(int Lp[], int Lj[], double Lx[], int Up[], int Ui[], double Ux[], 109 int P[], int Q[], double Dx[], int *do_recip, double Rs[], void *Numeric) 110 { 111 return umfpack_di_get_numeric(Lp,Lj,Lx,Up,Ui,Ux,P,Q,Dx,do_recip,Rs,Numeric); 112 } 113 114 inline int umfpack_get_numeric(int Lp[], int Lj[], std::complex<double> Lx[], int Up[], int Ui[], std::complex<double> Ux[], 115 int P[], int Q[], std::complex<double> Dx[], int *do_recip, double Rs[], void *Numeric) 116 { 117 double& lx0_real = numext::real_ref(Lx[0]); 118 double& ux0_real = numext::real_ref(Ux[0]); 119 double& dx0_real = numext::real_ref(Dx[0]); 120 return umfpack_zi_get_numeric(Lp,Lj,Lx?&lx0_real:0,0,Up,Ui,Ux?&ux0_real:0,0,P,Q, 121 Dx?&dx0_real:0,0,do_recip,Rs,Numeric); 122 } 123 124 inline int umfpack_get_determinant(double *Mx, double *Ex, void *NumericHandle, double User_Info [UMFPACK_INFO]) 125 { 126 return umfpack_di_get_determinant(Mx,Ex,NumericHandle,User_Info); 127 } 128 129 inline int umfpack_get_determinant(std::complex<double> *Mx, double *Ex, void *NumericHandle, double User_Info [UMFPACK_INFO]) 130 { 131 double& mx_real = numext::real_ref(*Mx); 132 return umfpack_zi_get_determinant(&mx_real,0,Ex,NumericHandle,User_Info); 133 } 134 135 136 /** \ingroup UmfPackSupport_Module 137 * \brief A sparse LU factorization and solver based on UmfPack 138 * 139 * This class allows to solve for A.X = B sparse linear problems via a LU factorization 140 * using the UmfPack library. The sparse matrix A must be squared and full rank. 141 * The vectors or matrices X and B can be either dense or sparse. 142 * 143 * \warning The input matrix A should be in a \b compressed and \b column-major form. 144 * Otherwise an expensive copy will be made. You can call the inexpensive makeCompressed() to get a compressed matrix. 145 * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> 146 * 147 * \implsparsesolverconcept 148 * 149 * \sa \ref TutorialSparseSolverConcept, class SparseLU 150 */ 151 template<typename _MatrixType> 152 class UmfPackLU : public SparseSolverBase<UmfPackLU<_MatrixType> > 153 { 154 protected: 155 typedef SparseSolverBase<UmfPackLU<_MatrixType> > Base; 156 using Base::m_isInitialized; 157 public: 158 using Base::_solve_impl; 159 typedef _MatrixType MatrixType; 160 typedef typename MatrixType::Scalar Scalar; 161 typedef typename MatrixType::RealScalar RealScalar; 162 typedef typename MatrixType::StorageIndex StorageIndex; 163 typedef Matrix<Scalar,Dynamic,1> Vector; 164 typedef Matrix<int, 1, MatrixType::ColsAtCompileTime> IntRowVectorType; 165 typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType; 166 typedef SparseMatrix<Scalar> LUMatrixType; 167 typedef SparseMatrix<Scalar,ColMajor,int> UmfpackMatrixType; 168 typedef Ref<const UmfpackMatrixType, StandardCompressedFormat> UmfpackMatrixRef; 169 enum { 170 ColsAtCompileTime = MatrixType::ColsAtCompileTime, 171 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime 172 }; 173 174 public: 175 176 typedef Array<double, UMFPACK_CONTROL, 1> UmfpackControl; 177 typedef Array<double, UMFPACK_INFO, 1> UmfpackInfo; 178 179 UmfPackLU() 180 : m_dummy(0,0), mp_matrix(m_dummy) 181 { 182 init(); 183 } 184 185 template<typename InputMatrixType> 186 explicit UmfPackLU(const InputMatrixType& matrix) 187 : mp_matrix(matrix) 188 { 189 init(); 190 compute(matrix); 191 } 192 193 ~UmfPackLU() 194 { 195 if(m_symbolic) umfpack_free_symbolic(&m_symbolic,Scalar()); 196 if(m_numeric) umfpack_free_numeric(&m_numeric,Scalar()); 197 } 198 199 inline Index rows() const { return mp_matrix.rows(); } 200 inline Index cols() const { return mp_matrix.cols(); } 201 202 /** \brief Reports whether previous computation was successful. 203 * 204 * \returns \c Success if computation was succesful, 205 * \c NumericalIssue if the matrix.appears to be negative. 206 */ 207 ComputationInfo info() const 208 { 209 eigen_assert(m_isInitialized && "Decomposition is not initialized."); 210 return m_info; 211 } 212 213 inline const LUMatrixType& matrixL() const 214 { 215 if (m_extractedDataAreDirty) extractData(); 216 return m_l; 217 } 218 219 inline const LUMatrixType& matrixU() const 220 { 221 if (m_extractedDataAreDirty) extractData(); 222 return m_u; 223 } 224 225 inline const IntColVectorType& permutationP() const 226 { 227 if (m_extractedDataAreDirty) extractData(); 228 return m_p; 229 } 230 231 inline const IntRowVectorType& permutationQ() const 232 { 233 if (m_extractedDataAreDirty) extractData(); 234 return m_q; 235 } 236 237 /** Computes the sparse Cholesky decomposition of \a matrix 238 * Note that the matrix should be column-major, and in compressed format for best performance. 239 * \sa SparseMatrix::makeCompressed(). 240 */ 241 template<typename InputMatrixType> 242 void compute(const InputMatrixType& matrix) 243 { 244 if(m_symbolic) umfpack_free_symbolic(&m_symbolic,Scalar()); 245 if(m_numeric) umfpack_free_numeric(&m_numeric,Scalar()); 246 grab(matrix.derived()); 247 analyzePattern_impl(); 248 factorize_impl(); 249 } 250 251 /** Performs a symbolic decomposition on the sparcity of \a matrix. 252 * 253 * This function is particularly useful when solving for several problems having the same structure. 254 * 255 * \sa factorize(), compute() 256 */ 257 template<typename InputMatrixType> 258 void analyzePattern(const InputMatrixType& matrix) 259 { 260 if(m_symbolic) umfpack_free_symbolic(&m_symbolic,Scalar()); 261 if(m_numeric) umfpack_free_numeric(&m_numeric,Scalar()); 262 263 grab(matrix.derived()); 264 265 analyzePattern_impl(); 266 } 267 268 /** Provides the return status code returned by UmfPack during the numeric 269 * factorization. 270 * 271 * \sa factorize(), compute() 272 */ 273 inline int umfpackFactorizeReturncode() const 274 { 275 eigen_assert(m_numeric && "UmfPackLU: you must first call factorize()"); 276 return m_fact_errorCode; 277 } 278 279 /** Provides access to the control settings array used by UmfPack. 280 * 281 * If this array contains NaN's, the default values are used. 282 * 283 * See UMFPACK documentation for details. 284 */ 285 inline const UmfpackControl& umfpackControl() const 286 { 287 return m_control; 288 } 289 290 /** Provides access to the control settings array used by UmfPack. 291 * 292 * If this array contains NaN's, the default values are used. 293 * 294 * See UMFPACK documentation for details. 295 */ 296 inline UmfpackControl& umfpackControl() 297 { 298 return m_control; 299 } 300 301 /** Performs a numeric decomposition of \a matrix 302 * 303 * The given matrix must has the same sparcity than the matrix on which the pattern anylysis has been performed. 304 * 305 * \sa analyzePattern(), compute() 306 */ 307 template<typename InputMatrixType> 308 void factorize(const InputMatrixType& matrix) 309 { 310 eigen_assert(m_analysisIsOk && "UmfPackLU: you must first call analyzePattern()"); 311 if(m_numeric) 312 umfpack_free_numeric(&m_numeric,Scalar()); 313 314 grab(matrix.derived()); 315 316 factorize_impl(); 317 } 318 319 /** Prints the current UmfPack control settings. 320 * 321 * \sa umfpackControl() 322 */ 323 void umfpackReportControl() 324 { 325 umfpack_report_control(m_control.data(), Scalar()); 326 } 327 328 /** Prints statistics collected by UmfPack. 329 * 330 * \sa analyzePattern(), compute() 331 */ 332 void umfpackReportInfo() 333 { 334 eigen_assert(m_analysisIsOk && "UmfPackLU: you must first call analyzePattern()"); 335 umfpack_report_info(m_control.data(), m_umfpackInfo.data(), Scalar()); 336 } 337 338 /** Prints the status of the previous factorization operation performed by UmfPack (symbolic or numerical factorization). 339 * 340 * \sa analyzePattern(), compute() 341 */ 342 void umfpackReportStatus() { 343 eigen_assert(m_analysisIsOk && "UmfPackLU: you must first call analyzePattern()"); 344 umfpack_report_status(m_control.data(), m_fact_errorCode, Scalar()); 345 } 346 347 /** \internal */ 348 template<typename BDerived,typename XDerived> 349 bool _solve_impl(const MatrixBase<BDerived> &b, MatrixBase<XDerived> &x) const; 350 351 Scalar determinant() const; 352 353 void extractData() const; 354 355 protected: 356 357 void init() 358 { 359 m_info = InvalidInput; 360 m_isInitialized = false; 361 m_numeric = 0; 362 m_symbolic = 0; 363 m_extractedDataAreDirty = true; 364 365 umfpack_defaults(m_control.data(), Scalar()); 366 } 367 368 void analyzePattern_impl() 369 { 370 m_fact_errorCode = umfpack_symbolic(internal::convert_index<int>(mp_matrix.rows()), 371 internal::convert_index<int>(mp_matrix.cols()), 372 mp_matrix.outerIndexPtr(), mp_matrix.innerIndexPtr(), mp_matrix.valuePtr(), 373 &m_symbolic, m_control.data(), m_umfpackInfo.data()); 374 375 m_isInitialized = true; 376 m_info = m_fact_errorCode ? InvalidInput : Success; 377 m_analysisIsOk = true; 378 m_factorizationIsOk = false; 379 m_extractedDataAreDirty = true; 380 } 381 382 void factorize_impl() 383 { 384 385 m_fact_errorCode = umfpack_numeric(mp_matrix.outerIndexPtr(), mp_matrix.innerIndexPtr(), mp_matrix.valuePtr(), 386 m_symbolic, &m_numeric, m_control.data(), m_umfpackInfo.data()); 387 388 m_info = m_fact_errorCode == UMFPACK_OK ? Success : NumericalIssue; 389 m_factorizationIsOk = true; 390 m_extractedDataAreDirty = true; 391 } 392 393 template<typename MatrixDerived> 394 void grab(const EigenBase<MatrixDerived> &A) 395 { 396 mp_matrix.~UmfpackMatrixRef(); 397 ::new (&mp_matrix) UmfpackMatrixRef(A.derived()); 398 } 399 400 void grab(const UmfpackMatrixRef &A) 401 { 402 if(&(A.derived()) != &mp_matrix) 403 { 404 mp_matrix.~UmfpackMatrixRef(); 405 ::new (&mp_matrix) UmfpackMatrixRef(A); 406 } 407 } 408 409 // cached data to reduce reallocation, etc. 410 mutable LUMatrixType m_l; 411 int m_fact_errorCode; 412 UmfpackControl m_control; 413 mutable UmfpackInfo m_umfpackInfo; 414 415 mutable LUMatrixType m_u; 416 mutable IntColVectorType m_p; 417 mutable IntRowVectorType m_q; 418 419 UmfpackMatrixType m_dummy; 420 UmfpackMatrixRef mp_matrix; 421 422 void* m_numeric; 423 void* m_symbolic; 424 425 mutable ComputationInfo m_info; 426 int m_factorizationIsOk; 427 int m_analysisIsOk; 428 mutable bool m_extractedDataAreDirty; 429 430 private: 431 UmfPackLU(const UmfPackLU& ) { } 432 }; 433 434 435 template<typename MatrixType> 436 void UmfPackLU<MatrixType>::extractData() const 437 { 438 if (m_extractedDataAreDirty) 439 { 440 // get size of the data 441 int lnz, unz, rows, cols, nz_udiag; 442 umfpack_get_lunz(&lnz, &unz, &rows, &cols, &nz_udiag, m_numeric, Scalar()); 443 444 // allocate data 445 m_l.resize(rows,(std::min)(rows,cols)); 446 m_l.resizeNonZeros(lnz); 447 448 m_u.resize((std::min)(rows,cols),cols); 449 m_u.resizeNonZeros(unz); 450 451 m_p.resize(rows); 452 m_q.resize(cols); 453 454 // extract 455 umfpack_get_numeric(m_l.outerIndexPtr(), m_l.innerIndexPtr(), m_l.valuePtr(), 456 m_u.outerIndexPtr(), m_u.innerIndexPtr(), m_u.valuePtr(), 457 m_p.data(), m_q.data(), 0, 0, 0, m_numeric); 458 459 m_extractedDataAreDirty = false; 460 } 461 } 462 463 template<typename MatrixType> 464 typename UmfPackLU<MatrixType>::Scalar UmfPackLU<MatrixType>::determinant() const 465 { 466 Scalar det; 467 umfpack_get_determinant(&det, 0, m_numeric, 0); 468 return det; 469 } 470 471 template<typename MatrixType> 472 template<typename BDerived,typename XDerived> 473 bool UmfPackLU<MatrixType>::_solve_impl(const MatrixBase<BDerived> &b, MatrixBase<XDerived> &x) const 474 { 475 Index rhsCols = b.cols(); 476 eigen_assert((BDerived::Flags&RowMajorBit)==0 && "UmfPackLU backend does not support non col-major rhs yet"); 477 eigen_assert((XDerived::Flags&RowMajorBit)==0 && "UmfPackLU backend does not support non col-major result yet"); 478 eigen_assert(b.derived().data() != x.derived().data() && " Umfpack does not support inplace solve"); 479 480 int errorCode; 481 Scalar* x_ptr = 0; 482 Matrix<Scalar,Dynamic,1> x_tmp; 483 if(x.innerStride()!=1) 484 { 485 x_tmp.resize(x.rows()); 486 x_ptr = x_tmp.data(); 487 } 488 for (int j=0; j<rhsCols; ++j) 489 { 490 if(x.innerStride()==1) 491 x_ptr = &x.col(j).coeffRef(0); 492 errorCode = umfpack_solve(UMFPACK_A, 493 mp_matrix.outerIndexPtr(), mp_matrix.innerIndexPtr(), mp_matrix.valuePtr(), 494 x_ptr, &b.const_cast_derived().col(j).coeffRef(0), m_numeric, m_control.data(), m_umfpackInfo.data()); 495 if(x.innerStride()!=1) 496 x.col(j) = x_tmp; 497 if (errorCode!=0) 498 return false; 499 } 500 501 return true; 502 } 503 504 } // end namespace Eigen 505 506 #endif // EIGEN_UMFPACKSUPPORT_H 507