| 494 | */ |
| 495 | template <typename MatrixType, typename OrderingType> |
| 496 | void SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix) |
| 497 | { |
| 498 | using internal::emptyIdxLU; |
| 499 | eigen_assert(m_analysisIsOk && "analyzePattern() should be called first"); |
| 500 | eigen_assert((matrix.rows() == matrix.cols()) && "Only for squared matrices"); |
| 501 | |
| 502 | typedef typename IndexVector::Scalar StorageIndex; |
| 503 | |
| 504 | m_isInitialized = true; |
| 505 | |
| 506 | |
| 507 | // Apply the column permutation computed in analyzepattern() |
| 508 | // m_mat = matrix * m_perm_c.inverse(); |
| 509 | m_mat = matrix; |
| 510 | if (m_perm_c.size()) |
| 511 | { |
| 512 | m_mat.uncompress(); //NOTE: The effect of this command is only to create the InnerNonzeros pointers. |
| 513 | //Then, permute only the column pointers |
| 514 | const StorageIndex * outerIndexPtr; |
| 515 | if (matrix.isCompressed()) outerIndexPtr = matrix.outerIndexPtr(); |
| 516 | else |
| 517 | { |
| 518 | StorageIndex* outerIndexPtr_t = new StorageIndex[matrix.cols()+1]; |
| 519 | for(Index i = 0; i <= matrix.cols(); i++) outerIndexPtr_t[i] = m_mat.outerIndexPtr()[i]; |
| 520 | outerIndexPtr = outerIndexPtr_t; |
| 521 | } |
| 522 | for (Index i = 0; i < matrix.cols(); i++) |
| 523 | { |
| 524 | m_mat.outerIndexPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i]; |
| 525 | m_mat.innerNonZeroPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i+1] - outerIndexPtr[i]; |
| 526 | } |
| 527 | if(!matrix.isCompressed()) delete[] outerIndexPtr; |
| 528 | } |
| 529 | else |
| 530 | { //FIXME This should not be needed if the empty permutation is handled transparently |
| 531 | m_perm_c.resize(matrix.cols()); |
| 532 | for(StorageIndex i = 0; i < matrix.cols(); ++i) m_perm_c.indices()(i) = i; |
| 533 | } |
| 534 | |
| 535 | Index m = m_mat.rows(); |
| 536 | Index n = m_mat.cols(); |
| 537 | Index nnz = m_mat.nonZeros(); |
| 538 | Index maxpanel = m_perfv.panel_size * m; |
| 539 | // Allocate working storage common to the factor routines |
| 540 | Index lwork = 0; |
| 541 | Index info = Base::memInit(m, n, nnz, lwork, m_perfv.fillfactor, m_perfv.panel_size, m_glu); |
| 542 | if (info) |
| 543 | { |
| 544 | m_lastError = "UNABLE TO ALLOCATE WORKING MEMORY\n\n" ; |
| 545 | m_factorizationIsOk = false; |
| 546 | return ; |
| 547 | } |
| 548 | |
| 549 | // Set up pointers for integer working arrays |
| 550 | IndexVector segrep(m); segrep.setZero(); |
| 551 | IndexVector parent(m); parent.setZero(); |
| 552 | IndexVector xplore(m); xplore.setZero(); |
| 553 | IndexVector repfnz(maxpanel); |