KLUSupport.h (11555B)
1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2017 Kyle Macfarlan <kyle.macfarlan@gmail.com> 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_KLUSUPPORT_H 11 #define EIGEN_KLUSUPPORT_H 12 13 namespace Eigen { 14 15 /* TODO extract L, extract U, compute det, etc... */ 16 17 /** \ingroup KLUSupport_Module 18 * \brief A sparse LU factorization and solver based on KLU 19 * 20 * This class allows to solve for A.X = B sparse linear problems via a LU factorization 21 * using the KLU library. The sparse matrix A must be squared and full rank. 22 * The vectors or matrices X and B can be either dense or sparse. 23 * 24 * \warning The input matrix A should be in a \b compressed and \b column-major form. 25 * Otherwise an expensive copy will be made. You can call the inexpensive makeCompressed() to get a compressed matrix. 26 * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> 27 * 28 * \implsparsesolverconcept 29 * 30 * \sa \ref TutorialSparseSolverConcept, class UmfPackLU, class SparseLU 31 */ 32 33 34 inline int klu_solve(klu_symbolic *Symbolic, klu_numeric *Numeric, Index ldim, Index nrhs, double B [ ], klu_common *Common, double) { 35 return klu_solve(Symbolic, Numeric, internal::convert_index<int>(ldim), internal::convert_index<int>(nrhs), B, Common); 36 } 37 38 inline int klu_solve(klu_symbolic *Symbolic, klu_numeric *Numeric, Index ldim, Index nrhs, std::complex<double>B[], klu_common *Common, std::complex<double>) { 39 return klu_z_solve(Symbolic, Numeric, internal::convert_index<int>(ldim), internal::convert_index<int>(nrhs), &numext::real_ref(B[0]), Common); 40 } 41 42 inline int klu_tsolve(klu_symbolic *Symbolic, klu_numeric *Numeric, Index ldim, Index nrhs, double B[], klu_common *Common, double) { 43 return klu_tsolve(Symbolic, Numeric, internal::convert_index<int>(ldim), internal::convert_index<int>(nrhs), B, Common); 44 } 45 46 inline int klu_tsolve(klu_symbolic *Symbolic, klu_numeric *Numeric, Index ldim, Index nrhs, std::complex<double>B[], klu_common *Common, std::complex<double>) { 47 return klu_z_tsolve(Symbolic, Numeric, internal::convert_index<int>(ldim), internal::convert_index<int>(nrhs), &numext::real_ref(B[0]), 0, Common); 48 } 49 50 inline klu_numeric* klu_factor(int Ap [ ], int Ai [ ], double Ax [ ], klu_symbolic *Symbolic, klu_common *Common, double) { 51 return klu_factor(Ap, Ai, Ax, Symbolic, Common); 52 } 53 54 inline klu_numeric* klu_factor(int Ap[], int Ai[], std::complex<double> Ax[], klu_symbolic *Symbolic, klu_common *Common, std::complex<double>) { 55 return klu_z_factor(Ap, Ai, &numext::real_ref(Ax[0]), Symbolic, Common); 56 } 57 58 59 template<typename _MatrixType> 60 class KLU : public SparseSolverBase<KLU<_MatrixType> > 61 { 62 protected: 63 typedef SparseSolverBase<KLU<_MatrixType> > Base; 64 using Base::m_isInitialized; 65 public: 66 using Base::_solve_impl; 67 typedef _MatrixType MatrixType; 68 typedef typename MatrixType::Scalar Scalar; 69 typedef typename MatrixType::RealScalar RealScalar; 70 typedef typename MatrixType::StorageIndex StorageIndex; 71 typedef Matrix<Scalar,Dynamic,1> Vector; 72 typedef Matrix<int, 1, MatrixType::ColsAtCompileTime> IntRowVectorType; 73 typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType; 74 typedef SparseMatrix<Scalar> LUMatrixType; 75 typedef SparseMatrix<Scalar,ColMajor,int> KLUMatrixType; 76 typedef Ref<const KLUMatrixType, StandardCompressedFormat> KLUMatrixRef; 77 enum { 78 ColsAtCompileTime = MatrixType::ColsAtCompileTime, 79 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime 80 }; 81 82 public: 83 84 KLU() 85 : m_dummy(0,0), mp_matrix(m_dummy) 86 { 87 init(); 88 } 89 90 template<typename InputMatrixType> 91 explicit KLU(const InputMatrixType& matrix) 92 : mp_matrix(matrix) 93 { 94 init(); 95 compute(matrix); 96 } 97 98 ~KLU() 99 { 100 if(m_symbolic) klu_free_symbolic(&m_symbolic,&m_common); 101 if(m_numeric) klu_free_numeric(&m_numeric,&m_common); 102 } 103 104 EIGEN_CONSTEXPR inline Index rows() const EIGEN_NOEXCEPT { return mp_matrix.rows(); } 105 EIGEN_CONSTEXPR inline Index cols() const EIGEN_NOEXCEPT { return mp_matrix.cols(); } 106 107 /** \brief Reports whether previous computation was successful. 108 * 109 * \returns \c Success if computation was successful, 110 * \c NumericalIssue if the matrix.appears to be negative. 111 */ 112 ComputationInfo info() const 113 { 114 eigen_assert(m_isInitialized && "Decomposition is not initialized."); 115 return m_info; 116 } 117 #if 0 // not implemented yet 118 inline const LUMatrixType& matrixL() const 119 { 120 if (m_extractedDataAreDirty) extractData(); 121 return m_l; 122 } 123 124 inline const LUMatrixType& matrixU() const 125 { 126 if (m_extractedDataAreDirty) extractData(); 127 return m_u; 128 } 129 130 inline const IntColVectorType& permutationP() const 131 { 132 if (m_extractedDataAreDirty) extractData(); 133 return m_p; 134 } 135 136 inline const IntRowVectorType& permutationQ() const 137 { 138 if (m_extractedDataAreDirty) extractData(); 139 return m_q; 140 } 141 #endif 142 /** Computes the sparse Cholesky decomposition of \a matrix 143 * Note that the matrix should be column-major, and in compressed format for best performance. 144 * \sa SparseMatrix::makeCompressed(). 145 */ 146 template<typename InputMatrixType> 147 void compute(const InputMatrixType& matrix) 148 { 149 if(m_symbolic) klu_free_symbolic(&m_symbolic, &m_common); 150 if(m_numeric) klu_free_numeric(&m_numeric, &m_common); 151 grab(matrix.derived()); 152 analyzePattern_impl(); 153 factorize_impl(); 154 } 155 156 /** Performs a symbolic decomposition on the sparcity of \a matrix. 157 * 158 * This function is particularly useful when solving for several problems having the same structure. 159 * 160 * \sa factorize(), compute() 161 */ 162 template<typename InputMatrixType> 163 void analyzePattern(const InputMatrixType& matrix) 164 { 165 if(m_symbolic) klu_free_symbolic(&m_symbolic, &m_common); 166 if(m_numeric) klu_free_numeric(&m_numeric, &m_common); 167 168 grab(matrix.derived()); 169 170 analyzePattern_impl(); 171 } 172 173 174 /** Provides access to the control settings array used by KLU. 175 * 176 * See KLU documentation for details. 177 */ 178 inline const klu_common& kluCommon() const 179 { 180 return m_common; 181 } 182 183 /** Provides access to the control settings array used by UmfPack. 184 * 185 * If this array contains NaN's, the default values are used. 186 * 187 * See KLU documentation for details. 188 */ 189 inline klu_common& kluCommon() 190 { 191 return m_common; 192 } 193 194 /** Performs a numeric decomposition of \a matrix 195 * 196 * The given matrix must has the same sparcity than the matrix on which the pattern anylysis has been performed. 197 * 198 * \sa analyzePattern(), compute() 199 */ 200 template<typename InputMatrixType> 201 void factorize(const InputMatrixType& matrix) 202 { 203 eigen_assert(m_analysisIsOk && "KLU: you must first call analyzePattern()"); 204 if(m_numeric) 205 klu_free_numeric(&m_numeric,&m_common); 206 207 grab(matrix.derived()); 208 209 factorize_impl(); 210 } 211 212 /** \internal */ 213 template<typename BDerived,typename XDerived> 214 bool _solve_impl(const MatrixBase<BDerived> &b, MatrixBase<XDerived> &x) const; 215 216 #if 0 // not implemented yet 217 Scalar determinant() const; 218 219 void extractData() const; 220 #endif 221 222 protected: 223 224 void init() 225 { 226 m_info = InvalidInput; 227 m_isInitialized = false; 228 m_numeric = 0; 229 m_symbolic = 0; 230 m_extractedDataAreDirty = true; 231 232 klu_defaults(&m_common); 233 } 234 235 void analyzePattern_impl() 236 { 237 m_info = InvalidInput; 238 m_analysisIsOk = false; 239 m_factorizationIsOk = false; 240 m_symbolic = klu_analyze(internal::convert_index<int>(mp_matrix.rows()), 241 const_cast<StorageIndex*>(mp_matrix.outerIndexPtr()), const_cast<StorageIndex*>(mp_matrix.innerIndexPtr()), 242 &m_common); 243 if (m_symbolic) { 244 m_isInitialized = true; 245 m_info = Success; 246 m_analysisIsOk = true; 247 m_extractedDataAreDirty = true; 248 } 249 } 250 251 void factorize_impl() 252 { 253 254 m_numeric = klu_factor(const_cast<StorageIndex*>(mp_matrix.outerIndexPtr()), const_cast<StorageIndex*>(mp_matrix.innerIndexPtr()), const_cast<Scalar*>(mp_matrix.valuePtr()), 255 m_symbolic, &m_common, Scalar()); 256 257 258 m_info = m_numeric ? Success : NumericalIssue; 259 m_factorizationIsOk = m_numeric ? 1 : 0; 260 m_extractedDataAreDirty = true; 261 } 262 263 template<typename MatrixDerived> 264 void grab(const EigenBase<MatrixDerived> &A) 265 { 266 mp_matrix.~KLUMatrixRef(); 267 ::new (&mp_matrix) KLUMatrixRef(A.derived()); 268 } 269 270 void grab(const KLUMatrixRef &A) 271 { 272 if(&(A.derived()) != &mp_matrix) 273 { 274 mp_matrix.~KLUMatrixRef(); 275 ::new (&mp_matrix) KLUMatrixRef(A); 276 } 277 } 278 279 // cached data to reduce reallocation, etc. 280 #if 0 // not implemented yet 281 mutable LUMatrixType m_l; 282 mutable LUMatrixType m_u; 283 mutable IntColVectorType m_p; 284 mutable IntRowVectorType m_q; 285 #endif 286 287 KLUMatrixType m_dummy; 288 KLUMatrixRef mp_matrix; 289 290 klu_numeric* m_numeric; 291 klu_symbolic* m_symbolic; 292 klu_common m_common; 293 mutable ComputationInfo m_info; 294 int m_factorizationIsOk; 295 int m_analysisIsOk; 296 mutable bool m_extractedDataAreDirty; 297 298 private: 299 KLU(const KLU& ) { } 300 }; 301 302 #if 0 // not implemented yet 303 template<typename MatrixType> 304 void KLU<MatrixType>::extractData() const 305 { 306 if (m_extractedDataAreDirty) 307 { 308 eigen_assert(false && "KLU: extractData Not Yet Implemented"); 309 310 // get size of the data 311 int lnz, unz, rows, cols, nz_udiag; 312 umfpack_get_lunz(&lnz, &unz, &rows, &cols, &nz_udiag, m_numeric, Scalar()); 313 314 // allocate data 315 m_l.resize(rows,(std::min)(rows,cols)); 316 m_l.resizeNonZeros(lnz); 317 318 m_u.resize((std::min)(rows,cols),cols); 319 m_u.resizeNonZeros(unz); 320 321 m_p.resize(rows); 322 m_q.resize(cols); 323 324 // extract 325 umfpack_get_numeric(m_l.outerIndexPtr(), m_l.innerIndexPtr(), m_l.valuePtr(), 326 m_u.outerIndexPtr(), m_u.innerIndexPtr(), m_u.valuePtr(), 327 m_p.data(), m_q.data(), 0, 0, 0, m_numeric); 328 329 m_extractedDataAreDirty = false; 330 } 331 } 332 333 template<typename MatrixType> 334 typename KLU<MatrixType>::Scalar KLU<MatrixType>::determinant() const 335 { 336 eigen_assert(false && "KLU: extractData Not Yet Implemented"); 337 return Scalar(); 338 } 339 #endif 340 341 template<typename MatrixType> 342 template<typename BDerived,typename XDerived> 343 bool KLU<MatrixType>::_solve_impl(const MatrixBase<BDerived> &b, MatrixBase<XDerived> &x) const 344 { 345 Index rhsCols = b.cols(); 346 EIGEN_STATIC_ASSERT((XDerived::Flags&RowMajorBit)==0, THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES); 347 eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or analyzePattern()/factorize()"); 348 349 x = b; 350 int info = klu_solve(m_symbolic, m_numeric, b.rows(), rhsCols, x.const_cast_derived().data(), const_cast<klu_common*>(&m_common), Scalar()); 351 352 m_info = info!=0 ? Success : NumericalIssue; 353 return true; 354 } 355 356 } // end namespace Eigen 357 358 #endif // EIGEN_KLUSUPPORT_H