| 450 | */ |
| 451 | template <typename MatrixType, typename OrderingType> |
| 452 | void SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix) |
| 453 | { |
| 454 | using internal::emptyIdxLU; |
| 455 | eigen_assert(m_analysisIsOk && "analyzePattern() should be called first"); |
| 456 | eigen_assert((matrix.rows() == matrix.cols()) && "Only for squared matrices"); |
| 457 | |
| 458 | typedef typename IndexVector::Scalar Index; |
| 459 | |
| 460 | |
| 461 | // Apply the column permutation computed in analyzepattern() |
| 462 | // m_mat = matrix * m_perm_c.inverse(); |
| 463 | m_mat = matrix; |
| 464 | if (m_perm_c.size()) |
| 465 | { |
| 466 | m_mat.uncompress(); //NOTE: The effect of this command is only to create the InnerNonzeros pointers. |
| 467 | //Then, permute only the column pointers |
| 468 | const Index * outerIndexPtr; |
| 469 | if (matrix.isCompressed()) outerIndexPtr = matrix.outerIndexPtr(); |
| 470 | else |
| 471 | { |
| 472 | Index* outerIndexPtr_t = new Index[matrix.cols()+1]; |
| 473 | for(Index i = 0; i <= matrix.cols(); i++) outerIndexPtr_t[i] = m_mat.outerIndexPtr()[i]; |
| 474 | outerIndexPtr = outerIndexPtr_t; |
| 475 | } |
| 476 | for (Index i = 0; i < matrix.cols(); i++) |
| 477 | { |
| 478 | m_mat.outerIndexPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i]; |
| 479 | m_mat.innerNonZeroPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i+1] - outerIndexPtr[i]; |
| 480 | } |
| 481 | if(!matrix.isCompressed()) delete[] outerIndexPtr; |
| 482 | } |
| 483 | else |
| 484 | { //FIXME This should not be needed if the empty permutation is handled transparently |
| 485 | m_perm_c.resize(matrix.cols()); |
| 486 | for(Index i = 0; i < matrix.cols(); ++i) m_perm_c.indices()(i) = i; |
| 487 | } |
| 488 | |
| 489 | Index m = m_mat.rows(); |
| 490 | Index n = m_mat.cols(); |
| 491 | Index nnz = m_mat.nonZeros(); |
| 492 | Index maxpanel = m_perfv.panel_size * m; |
| 493 | // Allocate working storage common to the factor routines |
| 494 | Index lwork = 0; |
| 495 | Index info = Base::memInit(m, n, nnz, lwork, m_perfv.fillfactor, m_perfv.panel_size, m_glu); |
| 496 | if (info) |
| 497 | { |
| 498 | m_lastError = "UNABLE TO ALLOCATE WORKING MEMORY\n\n" ; |
| 499 | m_factorizationIsOk = false; |
| 500 | return ; |
| 501 | } |
| 502 | |
| 503 | // Set up pointers for integer working arrays |
| 504 | IndexVector segrep(m); segrep.setZero(); |
| 505 | IndexVector parent(m); parent.setZero(); |
| 506 | IndexVector xplore(m); xplore.setZero(); |
| 507 | IndexVector repfnz(maxpanel); |
| 508 | IndexVector panel_lsub(maxpanel); |
| 509 | IndexVector xprune(n); xprune.setZero(); |