1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2012 Dsir Nuentsa-Wakam <desire.nuentsa_wakam (at) inria.fr> 5 // Copyright (C) 2012 Gael Guennebaud <gael.guennebaud (at) inria.fr> 6 // 7 // This Source Code Form is subject to the terms of the Mozilla 8 // Public License v. 2.0. If a copy of the MPL was not distributed 9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 10 11 12 #ifndef EIGEN_SPARSE_LU_H 13 #define EIGEN_SPARSE_LU_H 14 15 namespace Eigen { 16 17 template <typename _MatrixType, typename _OrderingType = COLAMDOrdering<typename _MatrixType::Index> > class SparseLU; 18 template <typename MappedSparseMatrixType> struct SparseLUMatrixLReturnType; 19 template <typename MatrixLType, typename MatrixUType> struct SparseLUMatrixUReturnType; 20 21 /** \ingroup SparseLU_Module 22 * \class SparseLU 23 * 24 * \brief Sparse supernodal LU factorization for general matrices 25 * 26 * This class implements the supernodal LU factorization for general matrices. 27 * It uses the main techniques from the sequential SuperLU package 28 * (http://crd-legacy.lbl.gov/~xiaoye/SuperLU/). It handles transparently real 29 * and complex arithmetics with single and double precision, depending on the 30 * scalar type of your input matrix. 31 * The code has been optimized to provide BLAS-3 operations during supernode-panel updates. 32 * It benefits directly from the built-in high-performant Eigen BLAS routines. 33 * Moreover, when the size of a supernode is very small, the BLAS calls are avoided to 34 * enable a better optimization from the compiler. For best performance, 35 * you should compile it with NDEBUG flag to avoid the numerous bounds checking on vectors. 36 * 37 * An important parameter of this class is the ordering method. It is used to reorder the columns 38 * (and eventually the rows) of the matrix to reduce the number of new elements that are created during 39 * numerical factorization. The cheapest method available is COLAMD. 40 * See \link OrderingMethods_Module the OrderingMethods module \endlink for the list of 41 * built-in and external ordering methods. 42 * 43 * Simple example with key steps 44 * \code 45 * VectorXd x(n), b(n); 46 * SparseMatrix<double, ColMajor> A; 47 * SparseLU<SparseMatrix<scalar, ColMajor>, COLAMDOrdering<Index> > solver; 48 * // fill A and b; 49 * // Compute the ordering permutation vector from the structural pattern of A 50 * solver.analyzePattern(A); 51 * // Compute the numerical factorization 52 * solver.factorize(A); 53 * //Use the factors to solve the linear system 54 * x = solver.solve(b); 55 * \endcode 56 * 57 * \warning The input matrix A should be in a \b compressed and \b column-major form. 58 * Otherwise an expensive copy will be made. You can call the inexpensive makeCompressed() to get a compressed matrix. 59 * 60 * \note Unlike the initial SuperLU implementation, there is no step to equilibrate the matrix. 61 * For badly scaled matrices, this step can be useful to reduce the pivoting during factorization. 62 * If this is the case for your matrices, you can try the basic scaling method at 63 * "unsupported/Eigen/src/IterativeSolvers/Scaling.h" 64 * 65 * \tparam _MatrixType The type of the sparse matrix. It must be a column-major SparseMatrix<> 66 * \tparam _OrderingType The ordering method to use, either AMD, COLAMD or METIS. Default is COLMAD 67 * 68 * 69 * \sa \ref TutorialSparseDirectSolvers 70 * \sa \ref OrderingMethods_Module 71 */ 72 template <typename _MatrixType, typename _OrderingType> 73 class SparseLU : public internal::SparseLUImpl<typename _MatrixType::Scalar, typename _MatrixType::Index> 74 { 75 public: 76 typedef _MatrixType MatrixType; 77 typedef _OrderingType OrderingType; 78 typedef typename MatrixType::Scalar Scalar; 79 typedef typename MatrixType::RealScalar RealScalar; 80 typedef typename MatrixType::Index Index; 81 typedef SparseMatrix<Scalar,ColMajor,Index> NCMatrix; 82 typedef internal::MappedSuperNodalMatrix<Scalar, Index> SCMatrix; 83 typedef Matrix<Scalar,Dynamic,1> ScalarVector; 84 typedef Matrix<Index,Dynamic,1> IndexVector; 85 typedef PermutationMatrix<Dynamic, Dynamic, Index> PermutationType; 86 typedef internal::SparseLUImpl<Scalar, Index> Base; 87 88 public: 89 SparseLU():m_isInitialized(true),m_lastError(""),m_Ustore(0,0,0,0,0,0),m_symmetricmode(false),m_diagpivotthresh(1.0),m_detPermR(1) 90 { 91 initperfvalues(); 92 } 93 SparseLU(const MatrixType& matrix):m_isInitialized(true),m_lastError(""),m_Ustore(0,0,0,0,0,0),m_symmetricmode(false),m_diagpivotthresh(1.0),m_detPermR(1) 94 { 95 initperfvalues(); 96 compute(matrix); 97 } 98 99 ~SparseLU() 100 { 101 // Free all explicit dynamic pointers 102 } 103 104 void analyzePattern (const MatrixType& matrix); 105 void factorize (const MatrixType& matrix); 106 void simplicialfactorize(const MatrixType& matrix); 107 108 /** 109 * Compute the symbolic and numeric factorization of the input sparse matrix. 110 * The input matrix should be in column-major storage. 111 */ 112 void compute (const MatrixType& matrix) 113 { 114 // Analyze 115 analyzePattern(matrix); 116 //Factorize 117 factorize(matrix); 118 } 119 120 inline Index rows() const { return m_mat.rows(); } 121 inline Index cols() const { return m_mat.cols(); } 122 /** Indicate that the pattern of the input matrix is symmetric */ 123 void isSymmetric(bool sym) 124 { 125 m_symmetricmode = sym; 126 } 127 128 /** \returns an expression of the matrix L, internally stored as supernodes 129 * The only operation available with this expression is the triangular solve 130 * \code 131 * y = b; matrixL().solveInPlace(y); 132 * \endcode 133 */ 134 SparseLUMatrixLReturnType<SCMatrix> matrixL() const 135 { 136 return SparseLUMatrixLReturnType<SCMatrix>(m_Lstore); 137 } 138 /** \returns an expression of the matrix U, 139 * The only operation available with this expression is the triangular solve 140 * \code 141 * y = b; matrixU().solveInPlace(y); 142 * \endcode 143 */ 144 SparseLUMatrixUReturnType<SCMatrix,MappedSparseMatrix<Scalar,ColMajor,Index> > matrixU() const 145 { 146 return SparseLUMatrixUReturnType<SCMatrix, MappedSparseMatrix<Scalar,ColMajor,Index> >(m_Lstore, m_Ustore); 147 } 148 149 /** 150 * \returns a reference to the row matrix permutation \f$ P_r \f$ such that \f$P_r A P_c^T = L U\f$ 151 * \sa colsPermutation() 152 */ 153 inline const PermutationType& rowsPermutation() const 154 { 155 return m_perm_r; 156 } 157 /** 158 * \returns a reference to the column matrix permutation\f$ P_c^T \f$ such that \f$P_r A P_c^T = L U\f$ 159 * \sa rowsPermutation() 160 */ 161 inline const PermutationType& colsPermutation() const 162 { 163 return m_perm_c; 164 } 165 /** Set the threshold used for a diagonal entry to be an acceptable pivot. */ 166 void setPivotThreshold(const RealScalar& thresh) 167 { 168 m_diagpivotthresh = thresh; 169 } 170 171 /** \returns the solution X of \f$ A X = B \f$ using the current decomposition of A. 172 * 173 * \warning the destination matrix X in X = this->solve(B) must be colmun-major. 174 * 175 * \sa compute() 176 */ 177 template<typename Rhs> 178 inline const internal::solve_retval<SparseLU, Rhs> solve(const MatrixBase<Rhs>& B) const 179 { 180 eigen_assert(m_factorizationIsOk && "SparseLU is not initialized."); 181 eigen_assert(rows()==B.rows() 182 && "SparseLU::solve(): invalid number of rows of the right hand side matrix B"); 183 return internal::solve_retval<SparseLU, Rhs>(*this, B.derived()); 184 } 185 186 /** \returns the solution X of \f$ A X = B \f$ using the current decomposition of A. 187 * 188 * \sa compute() 189 */ 190 template<typename Rhs> 191 inline const internal::sparse_solve_retval<SparseLU, Rhs> solve(const SparseMatrixBase<Rhs>& B) const 192 { 193 eigen_assert(m_factorizationIsOk && "SparseLU is not initialized."); 194 eigen_assert(rows()==B.rows() 195 && "SparseLU::solve(): invalid number of rows of the right hand side matrix B"); 196 return internal::sparse_solve_retval<SparseLU, Rhs>(*this, B.derived()); 197 } 198 199 /** \brief Reports whether previous computation was successful. 200 * 201 * \returns \c Success if computation was succesful, 202 * \c NumericalIssue if the LU factorization reports a problem, zero diagonal for instance 203 * \c InvalidInput if the input matrix is invalid 204 * 205 * \sa iparm() 206 */ 207 ComputationInfo info() const 208 { 209 eigen_assert(m_isInitialized && "Decomposition is not initialized."); 210 return m_info; 211 } 212 213 /** 214 * \returns A string describing the type of error 215 */ 216 std::string lastErrorMessage() const 217 { 218 return m_lastError; 219 } 220 221 template<typename Rhs, typename Dest> 222 bool _solve(const MatrixBase<Rhs> &B, MatrixBase<Dest> &X_base) const 223 { 224 Dest& X(X_base.derived()); 225 eigen_assert(m_factorizationIsOk && "The matrix should be factorized first"); 226 EIGEN_STATIC_ASSERT((Dest::Flags&RowMajorBit)==0, 227 THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES); 228 229 // Permute the right hand side to form X = Pr*B 230 // on return, X is overwritten by the computed solution 231 X.resize(B.rows(),B.cols()); 232 233 // this ugly const_cast_derived() helps to detect aliasing when applying the permutations 234 for(Index j = 0; j < B.cols(); ++j) 235 X.col(j) = rowsPermutation() * B.const_cast_derived().col(j); 236 237 //Forward substitution with L 238 this->matrixL().solveInPlace(X); 239 this->matrixU().solveInPlace(X); 240 241 // Permute back the solution 242 for (Index j = 0; j < B.cols(); ++j) 243 X.col(j) = colsPermutation().inverse() * X.col(j); 244 245 return true; 246 } 247 248 /** 249 * \returns the absolute value of the determinant of the matrix of which 250 * *this is the QR decomposition. 251 * 252 * \warning a determinant can be very big or small, so for matrices 253 * of large enough dimension, there is a risk of overflow/underflow. 254 * One way to work around that is to use logAbsDeterminant() instead. 255 * 256 * \sa logAbsDeterminant(), signDeterminant() 257 */ 258 Scalar absDeterminant() 259 { 260 eigen_assert(m_factorizationIsOk && "The matrix should be factorized first."); 261 // Initialize with the determinant of the row matrix 262 Scalar det = Scalar(1.); 263 //Note that the diagonal blocks of U are stored in supernodes, 264 // which are available in the L part :) 265 for (Index j = 0; j < this->cols(); ++j) 266 { 267 for (typename SCMatrix::InnerIterator it(m_Lstore, j); it; ++it) 268 { 269 if(it.row() < j) continue; 270 if(it.row() == j) 271 { 272 det *= (std::abs)(it.value()); 273 break; 274 } 275 } 276 } 277 return det; 278 } 279 280 /** \returns the natural log of the absolute value of the determinant of the matrix 281 * of which **this is the QR decomposition 282 * 283 * \note This method is useful to work around the risk of overflow/underflow that's 284 * inherent to the determinant computation. 285 * 286 * \sa absDeterminant(), signDeterminant() 287 */ 288 Scalar logAbsDeterminant() const 289 { 290 eigen_assert(m_factorizationIsOk && "The matrix should be factorized first."); 291 Scalar det = Scalar(0.); 292 for (Index j = 0; j < this->cols(); ++j) 293 { 294 for (typename SCMatrix::InnerIterator it(m_Lstore, j); it; ++it) 295 { 296 if(it.row() < j) continue; 297 if(it.row() == j) 298 { 299 det += (std::log)((std::abs)(it.value())); 300 break; 301 } 302 } 303 } 304 return det; 305 } 306 307 /** \returns A number representing the sign of the determinant 308 * 309 * \sa absDeterminant(), logAbsDeterminant() 310 */ 311 Scalar signDeterminant() 312 { 313 eigen_assert(m_factorizationIsOk && "The matrix should be factorized first."); 314 return Scalar(m_detPermR); 315 } 316 317 protected: 318 // Functions 319 void initperfvalues() 320 { 321 m_perfv.panel_size = 1; 322 m_perfv.relax = 1; 323 m_perfv.maxsuper = 128; 324 m_perfv.rowblk = 16; 325 m_perfv.colblk = 8; 326 m_perfv.fillfactor = 20; 327 } 328 329 // Variables 330 mutable ComputationInfo m_info; 331 bool m_isInitialized; 332 bool m_factorizationIsOk; 333 bool m_analysisIsOk; 334 std::string m_lastError; 335 NCMatrix m_mat; // The input (permuted ) matrix 336 SCMatrix m_Lstore; // The lower triangular matrix (supernodal) 337 MappedSparseMatrix<Scalar,ColMajor,Index> m_Ustore; // The upper triangular matrix 338 PermutationType m_perm_c; // Column permutation 339 PermutationType m_perm_r ; // Row permutation 340 IndexVector m_etree; // Column elimination tree 341 342 typename Base::GlobalLU_t m_glu; 343 344 // SparseLU options 345 bool m_symmetricmode; 346 // values for performance 347 internal::perfvalues<Index> m_perfv; 348 RealScalar m_diagpivotthresh; // Specifies the threshold used for a diagonal entry to be an acceptable pivot 349 Index m_nnzL, m_nnzU; // Nonzeros in L and U factors 350 Index m_detPermR; // Determinant of the coefficient matrix 351 private: 352 // Disable copy constructor 353 SparseLU (const SparseLU& ); 354 355 }; // End class SparseLU 356 357 358 359 // Functions needed by the anaysis phase 360 /** 361 * Compute the column permutation to minimize the fill-in 362 * 363 * - Apply this permutation to the input matrix - 364 * 365 * - Compute the column elimination tree on the permuted matrix 366 * 367 * - Postorder the elimination tree and the column permutation 368 * 369 */ 370 template <typename MatrixType, typename OrderingType> 371 void SparseLU<MatrixType, OrderingType>::analyzePattern(const MatrixType& mat) 372 { 373 374 //TODO It is possible as in SuperLU to compute row and columns scaling vectors to equilibrate the matrix mat. 375 376 OrderingType ord; 377 ord(mat,m_perm_c); 378 379 // Apply the permutation to the column of the input matrix 380 //First copy the whole input matrix. 381 m_mat = mat; 382 if (m_perm_c.size()) { 383 m_mat.uncompress(); //NOTE: The effect of this command is only to create the InnerNonzeros pointers. FIXME : This vector is filled but not subsequently used. 384 //Then, permute only the column pointers 385 const Index * outerIndexPtr; 386 if (mat.isCompressed()) outerIndexPtr = mat.outerIndexPtr(); 387 else 388 { 389 Index *outerIndexPtr_t = new Index[mat.cols()+1]; 390 for(Index i = 0; i <= mat.cols(); i++) outerIndexPtr_t[i] = m_mat.outerIndexPtr()[i]; 391 outerIndexPtr = outerIndexPtr_t; 392 } 393 for (Index i = 0; i < mat.cols(); i++) 394 { 395 m_mat.outerIndexPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i]; 396 m_mat.innerNonZeroPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i+1] - outerIndexPtr[i]; 397 } 398 if(!mat.isCompressed()) delete[] outerIndexPtr; 399 } 400 // Compute the column elimination tree of the permuted matrix 401 IndexVector firstRowElt; 402 internal::coletree(m_mat, m_etree,firstRowElt); 403 404 // In symmetric mode, do not do postorder here 405 if (!m_symmetricmode) { 406 IndexVector post, iwork; 407 // Post order etree 408 internal::treePostorder(m_mat.cols(), m_etree, post); 409 410 411 // Renumber etree in postorder 412 Index m = m_mat.cols(); 413 iwork.resize(m+1); 414 for (Index i = 0; i < m; ++i) iwork(post(i)) = post(m_etree(i)); 415 m_etree = iwork; 416 417 // Postmultiply A*Pc by post, i.e reorder the matrix according to the postorder of the etree 418 PermutationType post_perm(m); 419 for (Index i = 0; i < m; i++) 420 post_perm.indices()(i) = post(i); 421 422 // Combine the two permutations : postorder the permutation for future use 423 if(m_perm_c.size()) { 424 m_perm_c = post_perm * m_perm_c; 425 } 426 427 } // end postordering 428 429 m_analysisIsOk = true; 430 } 431 432 // Functions needed by the numerical factorization phase 433 434 435 /** 436 * - Numerical factorization 437 * - Interleaved with the symbolic factorization 438 * On exit, info is 439 * 440 * = 0: successful factorization 441 * 442 * > 0: if info = i, and i is 443 * 444 * <= A->ncol: U(i,i) is exactly zero. The factorization has 445 * been completed, but the factor U is exactly singular, 446 * and division by zero will occur if it is used to solve a 447 * system of equations. 448 * 449 * > A->ncol: number of bytes allocated when memory allocation 450 * failure occurred, plus A->ncol. If lwork = -1, it is 451 * the estimated amount of space needed, plus A->ncol. 452 */ 453 template <typename MatrixType, typename OrderingType> 454 void SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix) 455 { 456 using internal::emptyIdxLU; 457 eigen_assert(m_analysisIsOk && "analyzePattern() should be called first"); 458 eigen_assert((matrix.rows() == matrix.cols()) && "Only for squared matrices"); 459 460 typedef typename IndexVector::Scalar Index; 461 462 463 // Apply the column permutation computed in analyzepattern() 464 // m_mat = matrix * m_perm_c.inverse(); 465 m_mat = matrix; 466 if (m_perm_c.size()) 467 { 468 m_mat.uncompress(); //NOTE: The effect of this command is only to create the InnerNonzeros pointers. 469 //Then, permute only the column pointers 470 const Index * outerIndexPtr; 471 if (matrix.isCompressed()) outerIndexPtr = matrix.outerIndexPtr(); 472 else 473 { 474 Index* outerIndexPtr_t = new Index[matrix.cols()+1]; 475 for(Index i = 0; i <= matrix.cols(); i++) outerIndexPtr_t[i] = m_mat.outerIndexPtr()[i]; 476 outerIndexPtr = outerIndexPtr_t; 477 } 478 for (Index i = 0; i < matrix.cols(); i++) 479 { 480 m_mat.outerIndexPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i]; 481 m_mat.innerNonZeroPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i+1] - outerIndexPtr[i]; 482 } 483 if(!matrix.isCompressed()) delete[] outerIndexPtr; 484 } 485 else 486 { //FIXME This should not be needed if the empty permutation is handled transparently 487 m_perm_c.resize(matrix.cols()); 488 for(Index i = 0; i < matrix.cols(); ++i) m_perm_c.indices()(i) = i; 489 } 490 491 Index m = m_mat.rows(); 492 Index n = m_mat.cols(); 493 Index nnz = m_mat.nonZeros(); 494 Index maxpanel = m_perfv.panel_size * m; 495 // Allocate working storage common to the factor routines 496 Index lwork = 0; 497 Index info = Base::memInit(m, n, nnz, lwork, m_perfv.fillfactor, m_perfv.panel_size, m_glu); 498 if (info) 499 { 500 m_lastError = "UNABLE TO ALLOCATE WORKING MEMORY\n\n" ; 501 m_factorizationIsOk = false; 502 return ; 503 } 504 505 // Set up pointers for integer working arrays 506 IndexVector segrep(m); segrep.setZero(); 507 IndexVector parent(m); parent.setZero(); 508 IndexVector xplore(m); xplore.setZero(); 509 IndexVector repfnz(maxpanel); 510 IndexVector panel_lsub(maxpanel); 511 IndexVector xprune(n); xprune.setZero(); 512 IndexVector marker(m*internal::LUNoMarker); marker.setZero(); 513 514 repfnz.setConstant(-1); 515 panel_lsub.setConstant(-1); 516 517 // Set up pointers for scalar working arrays 518 ScalarVector dense; 519 dense.setZero(maxpanel); 520 ScalarVector tempv; 521 tempv.setZero(internal::LUnumTempV(m, m_perfv.panel_size, m_perfv.maxsuper, /*m_perfv.rowblk*/m) ); 522 523 // Compute the inverse of perm_c 524 PermutationType iperm_c(m_perm_c.inverse()); 525 526 // Identify initial relaxed snodes 527 IndexVector relax_end(n); 528 if ( m_symmetricmode == true ) 529 Base::heap_relax_snode(n, m_etree, m_perfv.relax, marker, relax_end); 530 else 531 Base::relax_snode(n, m_etree, m_perfv.relax, marker, relax_end); 532 533 534 m_perm_r.resize(m); 535 m_perm_r.indices().setConstant(-1); 536 marker.setConstant(-1); 537 m_detPermR = 1; // Record the determinant of the row permutation 538 539 m_glu.supno(0) = emptyIdxLU; m_glu.xsup.setConstant(0); 540 m_glu.xsup(0) = m_glu.xlsub(0) = m_glu.xusub(0) = m_glu.xlusup(0) = Index(0); 541 542 // Work on one 'panel' at a time. A panel is one of the following : 543 // (a) a relaxed supernode at the bottom of the etree, or 544 // (b) panel_size contiguous columns, <panel_size> defined by the user 545 Index jcol; 546 IndexVector panel_histo(n); 547 Index pivrow; // Pivotal row number in the original row matrix 548 Index nseg1; // Number of segments in U-column above panel row jcol 549 Index nseg; // Number of segments in each U-column 550 Index irep; 551 Index i, k, jj; 552 for (jcol = 0; jcol < n; ) 553 { 554 // Adjust panel size so that a panel won't overlap with the next relaxed snode. 555 Index panel_size = m_perfv.panel_size; // upper bound on panel width 556 for (k = jcol + 1; k < (std::min)(jcol+panel_size, n); k++) 557 { 558 if (relax_end(k) != emptyIdxLU) 559 { 560 panel_size = k - jcol; 561 break; 562 } 563 } 564 if (k == n) 565 panel_size = n - jcol; 566 567 // Symbolic outer factorization on a panel of columns 568 Base::panel_dfs(m, panel_size, jcol, m_mat, m_perm_r.indices(), nseg1, dense, panel_lsub, segrep, repfnz, xprune, marker, parent, xplore, m_glu); 569 570 // Numeric sup-panel updates in topological order 571 Base::panel_bmod(m, panel_size, jcol, nseg1, dense, tempv, segrep, repfnz, m_glu); 572 573 // Sparse LU within the panel, and below the panel diagonal 574 for ( jj = jcol; jj< jcol + panel_size; jj++) 575 { 576 k = (jj - jcol) * m; // Column index for w-wide arrays 577 578 nseg = nseg1; // begin after all the panel segments 579 //Depth-first-search for the current column 580 VectorBlock<IndexVector> panel_lsubk(panel_lsub, k, m); 581 VectorBlock<IndexVector> repfnz_k(repfnz, k, m); 582 info = Base::column_dfs(m, jj, m_perm_r.indices(), m_perfv.maxsuper, nseg, panel_lsubk, segrep, repfnz_k, xprune, marker, parent, xplore, m_glu); 583 if ( info ) 584 { 585 m_lastError = "UNABLE TO EXPAND MEMORY IN COLUMN_DFS() "; 586 m_info = NumericalIssue; 587 m_factorizationIsOk = false; 588 return; 589 } 590 // Numeric updates to this column 591 VectorBlock<ScalarVector> dense_k(dense, k, m); 592 VectorBlock<IndexVector> segrep_k(segrep, nseg1, m-nseg1); 593 info = Base::column_bmod(jj, (nseg - nseg1), dense_k, tempv, segrep_k, repfnz_k, jcol, m_glu); 594 if ( info ) 595 { 596 m_lastError = "UNABLE TO EXPAND MEMORY IN COLUMN_BMOD() "; 597 m_info = NumericalIssue; 598 m_factorizationIsOk = false; 599 return; 600 } 601 602 // Copy the U-segments to ucol(*) 603 info = Base::copy_to_ucol(jj, nseg, segrep, repfnz_k ,m_perm_r.indices(), dense_k, m_glu); 604 if ( info ) 605 { 606 m_lastError = "UNABLE TO EXPAND MEMORY IN COPY_TO_UCOL() "; 607 m_info = NumericalIssue; 608 m_factorizationIsOk = false; 609 return; 610 } 611 612 // Form the L-segment 613 info = Base::pivotL(jj, m_diagpivotthresh, m_perm_r.indices(), iperm_c.indices(), pivrow, m_glu); 614 if ( info ) 615 { 616 m_lastError = "THE MATRIX IS STRUCTURALLY SINGULAR ... ZERO COLUMN AT "; 617 std::ostringstream returnInfo; 618 returnInfo << info; 619 m_lastError += returnInfo.str(); 620 m_info = NumericalIssue; 621 m_factorizationIsOk = false; 622 return; 623 } 624 625 // Update the determinant of the row permutation matrix 626 if (pivrow != jj) m_detPermR *= -1; 627 628 // Prune columns (0:jj-1) using column jj 629 Base::pruneL(jj, m_perm_r.indices(), pivrow, nseg, segrep, repfnz_k, xprune, m_glu); 630 631 // Reset repfnz for this column 632 for (i = 0; i < nseg; i++) 633 { 634 irep = segrep(i); 635 repfnz_k(irep) = emptyIdxLU; 636 } 637 } // end SparseLU within the panel 638 jcol += panel_size; // Move to the next panel 639 } // end for -- end elimination 640 641 // Count the number of nonzeros in factors 642 Base::countnz(n, m_nnzL, m_nnzU, m_glu); 643 // Apply permutation to the L subscripts 644 Base::fixupL(n, m_perm_r.indices(), m_glu); 645 646 // Create supernode matrix L 647 m_Lstore.setInfos(m, n, m_glu.lusup, m_glu.xlusup, m_glu.lsub, m_glu.xlsub, m_glu.supno, m_glu.xsup); 648 // Create the column major upper sparse matrix U; 649 new (&m_Ustore) MappedSparseMatrix<Scalar, ColMajor, Index> ( m, n, m_nnzU, m_glu.xusub.data(), m_glu.usub.data(), m_glu.ucol.data() ); 650 651 m_info = Success; 652 m_factorizationIsOk = true; 653 } 654 655 template<typename MappedSupernodalType> 656 struct SparseLUMatrixLReturnType : internal::no_assignment_operator 657 { 658 typedef typename MappedSupernodalType::Index Index; 659 typedef typename MappedSupernodalType::Scalar Scalar; 660 SparseLUMatrixLReturnType(const MappedSupernodalType& mapL) : m_mapL(mapL) 661 { } 662 Index rows() { return m_mapL.rows(); } 663 Index cols() { return m_mapL.cols(); } 664 template<typename Dest> 665 void solveInPlace( MatrixBase<Dest> &X) const 666 { 667 m_mapL.solveInPlace(X); 668 } 669 const MappedSupernodalType& m_mapL; 670 }; 671 672 template<typename MatrixLType, typename MatrixUType> 673 struct SparseLUMatrixUReturnType : internal::no_assignment_operator 674 { 675 typedef typename MatrixLType::Index Index; 676 typedef typename MatrixLType::Scalar Scalar; 677 SparseLUMatrixUReturnType(const MatrixLType& mapL, const MatrixUType& mapU) 678 : m_mapL(mapL),m_mapU(mapU) 679 { } 680 Index rows() { return m_mapL.rows(); } 681 Index cols() { return m_mapL.cols(); } 682 683 template<typename Dest> void solveInPlace(MatrixBase<Dest> &X) const 684 { 685 Index nrhs = X.cols(); 686 Index n = X.rows(); 687 // Backward solve with U 688 for (Index k = m_mapL.nsuper(); k >= 0; k--) 689 { 690 Index fsupc = m_mapL.supToCol()[k]; 691 Index lda = m_mapL.colIndexPtr()[fsupc+1] - m_mapL.colIndexPtr()[fsupc]; // leading dimension 692 Index nsupc = m_mapL.supToCol()[k+1] - fsupc; 693 Index luptr = m_mapL.colIndexPtr()[fsupc]; 694 695 if (nsupc == 1) 696 { 697 for (Index j = 0; j < nrhs; j++) 698 { 699 X(fsupc, j) /= m_mapL.valuePtr()[luptr]; 700 } 701 } 702 else 703 { 704 Map<const Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride<> > A( &(m_mapL.valuePtr()[luptr]), nsupc, nsupc, OuterStride<>(lda) ); 705 Map< Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride<> > U (&(X(fsupc,0)), nsupc, nrhs, OuterStride<>(n) ); 706 U = A.template triangularView<Upper>().solve(U); 707 } 708 709 for (Index j = 0; j < nrhs; ++j) 710 { 711 for (Index jcol = fsupc; jcol < fsupc + nsupc; jcol++) 712 { 713 typename MatrixUType::InnerIterator it(m_mapU, jcol); 714 for ( ; it; ++it) 715 { 716 Index irow = it.index(); 717 X(irow, j) -= X(jcol, j) * it.value(); 718 } 719 } 720 } 721 } // End For U-solve 722 } 723 const MatrixLType& m_mapL; 724 const MatrixUType& m_mapU; 725 }; 726 727 namespace internal { 728 729 template<typename _MatrixType, typename Derived, typename Rhs> 730 struct solve_retval<SparseLU<_MatrixType,Derived>, Rhs> 731 : solve_retval_base<SparseLU<_MatrixType,Derived>, Rhs> 732 { 733 typedef SparseLU<_MatrixType,Derived> Dec; 734 EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) 735 736 template<typename Dest> void evalTo(Dest& dst) const 737 { 738 dec()._solve(rhs(),dst); 739 } 740 }; 741 742 template<typename _MatrixType, typename Derived, typename Rhs> 743 struct sparse_solve_retval<SparseLU<_MatrixType,Derived>, Rhs> 744 : sparse_solve_retval_base<SparseLU<_MatrixType,Derived>, Rhs> 745 { 746 typedef SparseLU<_MatrixType,Derived> Dec; 747 EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs) 748 749 template<typename Dest> void evalTo(Dest& dst) const 750 { 751 this->defaultEvalTo(dst); 752 } 753 }; 754 } // end namespace internal 755 756 } // End namespace Eigen 757 758 #endif 759