Method factorize

deps/Eigen/src/SparseLU/SparseLU.h:611–814

 */
template <typename MatrixType, typename OrderingType>
void SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix) {
  using internal::emptyIdxLU;
  eigen_assert(m_analysisIsOk && "analyzePattern() should be called first");
  eigen_assert((matrix.rows() == matrix.cols()) && "Only for squared matrices");

  m_isInitialized = true;

  // Apply the column permutation computed in analyzepattern()
  //   m_mat = matrix * m_perm_c.inverse();
  m_mat = matrix;
  if (m_perm_c.size()) {
    m_mat.uncompress();  // NOTE: The effect of this command is only to create the InnerNonzeros pointers.
    // Then, permute only the column pointers
    const StorageIndex* outerIndexPtr;
    if (matrix.isCompressed())
      outerIndexPtr = matrix.outerIndexPtr();
    else {
      StorageIndex* outerIndexPtr_t = new StorageIndex[matrix.cols() + 1];
      for (Index i = 0; i <= matrix.cols(); i++) outerIndexPtr_t[i] = m_mat.outerIndexPtr()[i];
      outerIndexPtr = outerIndexPtr_t;
    }
    for (Index i = 0; i < matrix.cols(); i++) {
      m_mat.outerIndexPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i];
      m_mat.innerNonZeroPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i + 1] - outerIndexPtr[i];
    }
    if (!matrix.isCompressed()) delete[] outerIndexPtr;
  } else {  // FIXME This should not be needed if the empty permutation is handled transparently
    m_perm_c.resize(matrix.cols());
    for (StorageIndex i = 0; i < matrix.cols(); ++i) m_perm_c.indices()(i) = i;
  }

  Index m = m_mat.rows();
  Index n = m_mat.cols();
  Index nnz = m_mat.nonZeros();
  Index maxpanel = m_perfv.panel_size * m;
  // Allocate working storage common to the factor routines
  Index lwork = 0;
  // Return the size of actually allocated memory when allocation failed,
  // and 0 on success.
  Index info = Base::memInit(m, n, nnz, lwork, m_perfv.fillfactor, m_perfv.panel_size, m_glu);
  if (info) {
    m_lastError = "UNABLE TO ALLOCATE WORKING MEMORY\n\n";
    m_factorizationIsOk = false;
    return;
  }

  // Set up pointers for integer working arrays
  IndexVector segrep(m);
  segrep.setZero();
  IndexVector parent(m);
  parent.setZero();
  IndexVector xplore(m);
  xplore.setZero();
  IndexVector repfnz(maxpanel);
  IndexVector panel_lsub(maxpanel);
  IndexVector xprune(n);
  xprune.setZero();
  IndexVector marker(m * internal::LUNoMarker);