SparseLU.h

// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr>
// Copyright (C) 2012-2014 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
// SPDX-License-Identifier: MPL-2.0
 
#ifndef EIGEN_SPARSE_LU_H
#define EIGEN_SPARSE_LU_H
 
// IWYU pragma: private
#include "./InternalHeaderCheck.h"
 
namespace Eigen {
 
template <typename MatrixType_, typename OrderingType_ = COLAMDOrdering<typename MatrixType_::StorageIndex>>
class SparseLU;
template <typename MappedSparseMatrixType>
struct SparseLUMatrixLReturnType;
template <typename MatrixLType, typename MatrixUType>
struct SparseLUMatrixUReturnType;
 
template <bool Conjugate, class SparseLUType>
class SparseLUTransposeView : public SparseSolverBase<SparseLUTransposeView<Conjugate, SparseLUType>> {
 protected:
  using APIBase = SparseSolverBase<SparseLUTransposeView<Conjugate, SparseLUType>>;
  using APIBase::m_isInitialized;
 
 public:
  using Scalar = typename SparseLUType::Scalar;
  using StorageIndex = typename SparseLUType::StorageIndex;
  using MatrixType = typename SparseLUType::MatrixType;
  using OrderingType = typename SparseLUType::OrderingType;
 
  enum { ColsAtCompileTime = MatrixType::ColsAtCompileTime, MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime };
 
  SparseLUTransposeView() = default;
  SparseLUTransposeView(const SparseLUTransposeView& view) : APIBase() {
    this->m_sparseLU = view.m_sparseLU;
    this->m_isInitialized = view.m_isInitialized;
  }
  void setIsInitialized(const bool isInitialized) { this->m_isInitialized = isInitialized; }
  void setSparseLU(SparseLUType* sparseLU) { m_sparseLU = sparseLU; }
  using APIBase::_solve_impl;
  template <typename Rhs, typename Dest>
  bool _solve_impl(const MatrixBase<Rhs>& B, MatrixBase<Dest>& X_base) const {
    Dest& X(X_base.derived());
    eigen_assert(m_sparseLU->info() == Success && "The matrix should be factorized first");
    EIGEN_STATIC_ASSERT((Dest::Flags & RowMajorBit) == 0, THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
 
    // const_cast_derived() is needed to enable aliasing detection when applying the permutations.
    for (Index j = 0; j < B.cols(); ++j) {
      X.col(j) = m_sparseLU->colsPermutation() * B.const_cast_derived().col(j);
    }
    // Forward substitution with transposed or adjoint of U
    m_sparseLU->matrixU().template solveTransposedInPlace<Conjugate>(X);
 
    // Backward substitution with transposed or adjoint of L
    m_sparseLU->matrixL().template solveTransposedInPlace<Conjugate>(X);
 
    // Permute back the solution
    for (Index j = 0; j < B.cols(); ++j) X.col(j) = m_sparseLU->rowsPermutation().transpose() * X.col(j);
    return true;
  }
  inline Index rows() const { return m_sparseLU->rows(); }
  inline Index cols() const { return m_sparseLU->cols(); }
 
 private:
  SparseLUType* m_sparseLU = nullptr;
  SparseLUTransposeView& operator=(const SparseLUTransposeView&) = delete;
};

template <typename MatrixType_, typename OrderingType_>
class SparseLU : public SparseSolverBase<SparseLU<MatrixType_, OrderingType_>>,
                 public internal::SparseLUImpl<typename MatrixType_::Scalar, typename MatrixType_::StorageIndex> {
 protected:
  using APIBase = SparseSolverBase<SparseLU<MatrixType_, OrderingType_>>;
  using APIBase::m_isInitialized;
 
 public:
  using APIBase::_solve_impl;
 
  using MatrixType = MatrixType_;
  using OrderingType = OrderingType_;
  using Scalar = typename MatrixType::Scalar;
  using RealScalar = typename MatrixType::RealScalar;
  using StorageIndex = typename MatrixType::StorageIndex;
  using NCMatrix = SparseMatrix<Scalar, ColMajor, StorageIndex>;
  using SCMatrix = internal::MappedSuperNodalMatrix<Scalar, StorageIndex>;
  using ScalarVector = Matrix<Scalar, Dynamic, 1>;
  using IndexVector = Matrix<StorageIndex, Dynamic, 1>;
  using PermutationType = PermutationMatrix<Dynamic, Dynamic, StorageIndex>;
  using Base = internal::SparseLUImpl<Scalar, StorageIndex>;
 
  enum { ColsAtCompileTime = MatrixType::ColsAtCompileTime, MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime };
 
 public:
  SparseLU()
      : m_lastError(""), m_Ustore(0, 0, 0, 0, 0, 0), m_symmetricmode(false), m_diagpivotthresh(1.0), m_detPermR(1) {
    initperfvalues();
  }

  explicit SparseLU(const MatrixType& matrix)
      : m_lastError(""), m_Ustore(0, 0, 0, 0, 0, 0), m_symmetricmode(false), m_diagpivotthresh(1.0), m_detPermR(1) {
    initperfvalues();
    compute(matrix);
  }
 
  ~SparseLU() {
    // Free all explicit dynamic pointers
  }
 
  void analyzePattern(const MatrixType& matrix);
  void factorize(const MatrixType& matrix);
  void simplicialfactorize(const MatrixType& matrix);

  void compute(const MatrixType& matrix) {
    // Analyze
    analyzePattern(matrix);
    // Factorize
    factorize(matrix);
  }

  const SparseLUTransposeView<false, SparseLU<MatrixType_, OrderingType_>> transpose() {
    SparseLUTransposeView<false, SparseLU<MatrixType_, OrderingType_>> transposeView;
    transposeView.setSparseLU(this);
    transposeView.setIsInitialized(this->m_isInitialized);
    return transposeView;
  }

  const SparseLUTransposeView<true, SparseLU<MatrixType_, OrderingType_>> adjoint() {
    SparseLUTransposeView<true, SparseLU<MatrixType_, OrderingType_>> adjointView;
    adjointView.setSparseLU(this);
    adjointView.setIsInitialized(this->m_isInitialized);
    return adjointView;
  }

  inline Index rows() const { return m_mat.rows(); }
  inline Index cols() const { return m_mat.cols(); }
  void isSymmetric(bool sym) { m_symmetricmode = sym; }

  SparseLUMatrixLReturnType<SCMatrix> matrixL() const { return SparseLUMatrixLReturnType<SCMatrix>(m_Lstore); }
  SparseLUMatrixUReturnType<SCMatrix, Map<SparseMatrix<Scalar, ColMajor, StorageIndex>>> matrixU() const {
    return SparseLUMatrixUReturnType<SCMatrix, Map<SparseMatrix<Scalar, ColMajor, StorageIndex>>>(m_Lstore, m_Ustore);
  }

  inline const PermutationType& rowsPermutation() const { return m_perm_r; }
  inline const PermutationType& colsPermutation() const { return m_perm_c; }
  void setPivotThreshold(const RealScalar& thresh) { m_diagpivotthresh = thresh; }
 
#ifdef EIGEN_PARSED_BY_DOXYGEN
  template <typename Rhs>
  inline Solve<SparseLU, Rhs> solve(const MatrixBase<Rhs>& B) const;
#endif  // EIGEN_PARSED_BY_DOXYGEN

  ComputationInfo info() const {
    eigen_assert(m_isInitialized && "Decomposition is not initialized.");
    return m_info;
  }

  std::string lastErrorMessage() const { return m_lastError; }
 
  template <typename Rhs, typename Dest>
  bool _solve_impl(const MatrixBase<Rhs>& B, MatrixBase<Dest>& X_base) const {
    Dest& X(X_base.derived());
    eigen_assert(m_factorizationIsOk && "The matrix should be factorized first");
    EIGEN_STATIC_ASSERT((Dest::Flags & RowMajorBit) == 0, THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
 
    // Permute the right hand side to form X = Pr*B
    // on return, X is overwritten by the computed solution
    X.resize(B.rows(), B.cols());
 
    // const_cast_derived() is needed to enable aliasing detection when applying the permutations.
    for (Index j = 0; j < B.cols(); ++j) X.col(j) = rowsPermutation() * B.const_cast_derived().col(j);
 
    // Forward substitution with L
    this->matrixL().solveInPlace(X);
    this->matrixU().solveInPlace(X);
 
    // Permute back the solution
    for (Index j = 0; j < B.cols(); ++j) X.col(j) = colsPermutation().inverse() * X.col(j);
 
    return true;
  }

  Scalar absDeterminant() const {
    using std::abs;
    eigen_assert(m_factorizationIsOk && "The matrix should be factorized first.");
    // Initialize with the determinant of the row matrix
    Scalar det = Scalar(1.);
    // Note that the diagonal blocks of U are stored in supernodes,
    // which are available in the  L part :)
    for (Index j = 0; j < this->cols(); ++j) {
      for (typename SCMatrix::InnerIterator it(m_Lstore, j); it; ++it) {
        if (it.index() == j) {
          det *= abs(it.value());
          break;
        }
      }
    }
    return det;
  }

  Scalar logAbsDeterminant() const {
    using std::abs;
    using std::log;
 
    eigen_assert(m_factorizationIsOk && "The matrix should be factorized first.");
    Scalar det = Scalar(0.);
    for (Index j = 0; j < this->cols(); ++j) {
      for (typename SCMatrix::InnerIterator it(m_Lstore, j); it; ++it) {
        if (it.row() < j) continue;
        if (it.row() == j) {
          det += log(abs(it.value()));
          break;
        }
      }
    }
    return det;
  }

  Scalar signDeterminant() const {
    eigen_assert(m_factorizationIsOk && "The matrix should be factorized first.");
    // Initialize with the determinant of the row matrix
    Index det = 1;
    // Note that the diagonal blocks of U are stored in supernodes,
    // which are available in the  L part :)
    for (Index j = 0; j < this->cols(); ++j) {
      for (typename SCMatrix::InnerIterator it(m_Lstore, j); it; ++it) {
        if (it.index() == j) {
          if (it.value() < 0)
            det = -det;
          else if (it.value() == 0)
            return 0;
          break;
        }
      }
    }
    return det * m_detPermR * m_detPermC;
  }

  Scalar determinant() const {
    eigen_assert(m_factorizationIsOk && "The matrix should be factorized first.");
    // Initialize with the determinant of the row matrix
    Scalar det = Scalar(1.);
    // Note that the diagonal blocks of U are stored in supernodes,
    // which are available in the  L part :)
    for (Index j = 0; j < this->cols(); ++j) {
      for (typename SCMatrix::InnerIterator it(m_Lstore, j); it; ++it) {
        if (it.index() == j) {
          det *= it.value();
          break;
        }
      }
    }
    return (m_detPermR * m_detPermC) > 0 ? det : -det;
  }

  Index nnzL() const { return m_nnzL; }
  Index nnzU() const { return m_nnzU; }
 
 protected:
  // Functions
  void initperfvalues() {
    m_perfv.panel_size = 16;
    m_perfv.relax = 1;
    m_perfv.maxsuper = 128;
    m_perfv.rowblk = 16;
    m_perfv.colblk = 8;
    m_perfv.fillfactor = 20;
  }
 
  // Variables
  mutable ComputationInfo m_info;
  bool m_factorizationIsOk;
  bool m_analysisIsOk;
  std::string m_lastError;
  NCMatrix m_mat;                                              // The input (permuted ) matrix
  SCMatrix m_Lstore;                                           // The lower triangular matrix (supernodal)
  Map<SparseMatrix<Scalar, ColMajor, StorageIndex>> m_Ustore;  // The upper triangular matrix
  PermutationType m_perm_c;                                    // Column permutation
  PermutationType m_perm_r;                                    // Row permutation
  IndexVector m_etree;                                         // Column elimination tree
 
  typename Base::GlobalLU_t m_glu;
 
  // SparseLU options
  bool m_symmetricmode;
  // values for performance
  internal::perfvalues m_perfv;
  RealScalar m_diagpivotthresh;  // Specifies the threshold used for a diagonal entry to be an acceptable pivot
  Index m_nnzL, m_nnzU;          // Nonzeros in L and U factors
  Index m_detPermR, m_detPermC;  // Determinants of the permutation matrices
 private:
  SparseLU(const SparseLU&) = delete;
};  // End class SparseLU
 
// Functions needed by the analysis phase
template <typename MatrixType, typename OrderingType>
void SparseLU<MatrixType, OrderingType>::analyzePattern(const MatrixType& mat) {
  // TODO  It is possible as in SuperLU to compute row and column scaling vectors to equilibrate the matrix mat.
 
  // Firstly, copy the whole input matrix.
  m_mat = mat;
 
  // Compute fill-in ordering
  OrderingType ord;
  ord(m_mat, m_perm_c);
 
  // Apply the permutation to the column of the input matrix
  if (m_perm_c.size()) {
    // Switch to uncompressed mode so innerNonZeroPtr() exists and can be
    // permuted consistently with outerIndexPtr().
    // Downstream sparse traversals may also rely on these per-column counts
    // while m_mat remains uncompressed.
    m_mat.uncompress();
    // A compressed column-major input already exposes valid column pointers.
    // Otherwise snapshot the internal column-major structure before permuting in place.
    const bool useInputOuterIndex = !MatrixType::IsRowMajor && mat.isCompressed();
    ei_declare_aligned_stack_constructed_variable(
        StorageIndex, outerIndexPtr, m_mat.cols() + 1,
        useInputOuterIndex ? const_cast<StorageIndex*>(mat.outerIndexPtr()) : 0);
    if (!useInputOuterIndex)
      IndexVector::Map(outerIndexPtr, m_mat.cols() + 1) = IndexVector::Map(m_mat.outerIndexPtr(), m_mat.cols() + 1);
 
    // Apply the permutation and compute the nnz per column.
    for (Index i = 0; i < mat.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];
    }
  }
 
  // Compute the column elimination tree of the permuted matrix
  IndexVector firstRowElt;
  internal::coletree(m_mat, m_etree, firstRowElt);
 
  // In symmetric mode, do not do postorder here
  if (!m_symmetricmode) {
    IndexVector post, iwork;
    // Post order etree
    internal::treePostorder(StorageIndex(m_mat.cols()), m_etree, post);
 
    // Renumber etree in postorder
    Index m = m_mat.cols();
    iwork.resize(m + 1);
    for (Index i = 0; i < m; ++i) iwork(post(i)) = post(m_etree(i));
    m_etree = iwork;
 
    // Postmultiply A*Pc by post, i.e reorder the matrix according to the postorder of the etree
    PermutationType post_perm(m);
    for (Index i = 0; i < m; i++) post_perm.indices()(i) = post(i);
 
    // Combine the two permutations : postorder the permutation for future use
    if (m_perm_c.size()) {
      m_perm_c = post_perm * m_perm_c;
    }
 
  }  // end postordering
 
  m_analysisIsOk = true;
}
 
// Functions needed by the numerical factorization phase

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;
 
  // Reset state from any prior factorize() so info() and lastErrorMessage()
  // describe this call's outcome, not the previous matrix's.
  m_info = Success;
  m_lastError.clear();
 
  // Apply the column permutation computed in analyzepattern()
  m_mat = matrix;
  if (m_perm_c.size()) {
    // Switch to uncompressed mode so innerNonZeroPtr() exists and can be
    // permuted consistently with outerIndexPtr().
    m_mat.uncompress();
    const bool useInputOuterIndex = !MatrixType::IsRowMajor && matrix.isCompressed();
    ei_declare_aligned_stack_constructed_variable(
        StorageIndex, outerIndexPtr, m_mat.cols() + 1,
        useInputOuterIndex ? const_cast<StorageIndex*>(matrix.outerIndexPtr()) : 0);
    if (!useInputOuterIndex)
      IndexVector::Map(outerIndexPtr, m_mat.cols() + 1) = IndexVector::Map(m_mat.outerIndexPtr(), m_mat.cols() + 1);
    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];
    }
  } 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);
  marker.setZero();
 
  repfnz.setConstant(-1);
  panel_lsub.setConstant(-1);
 
  // Set up pointers for scalar working arrays
  ScalarVector dense;
  dense.setZero(maxpanel);
  ScalarVector tempv;
  tempv.setZero(internal::LUnumTempV(m, m_perfv.panel_size, m_perfv.maxsuper, /*m_perfv.rowblk*/ m));
 
  // Compute the inverse of perm_c
  PermutationType iperm_c(m_perm_c.inverse());
 
  // Identify initial relaxed snodes
  IndexVector relax_end(n);
  if (m_symmetricmode == true)
    Base::heap_relax_snode(n, m_etree, m_perfv.relax, marker, relax_end);
  else
    Base::relax_snode(n, m_etree, m_perfv.relax, marker, relax_end);
 
  m_perm_r.resize(m);
  m_perm_r.indices().setConstant(-1);
  marker.setConstant(-1);
  m_detPermR = 1;  // Record the determinant of the row permutation
 
  m_glu.supno(0) = emptyIdxLU;
  m_glu.xsup.setConstant(0);
  m_glu.xsup(0) = m_glu.xlsub(0) = m_glu.xusub(0) = m_glu.xlusup(0) = Index(0);
 
  // Work on one 'panel' at a time. A panel is one of the following :
  //  (a) a relaxed supernode at the bottom of the etree, or
  //  (b) panel_size contiguous columns, <panel_size> defined by the user
  Index jcol;
  Index pivrow;  // Pivotal row number in the original row matrix
  Index nseg1;   // Number of segments in U-column above panel row jcol
  Index nseg;    // Number of segments in each U-column
  Index irep;
  Index i, k, jj;
  for (jcol = 0; jcol < n;) {
    // Adjust panel size so that a panel won't overlap with the next relaxed snode.
    Index panel_size = m_perfv.panel_size;  // upper bound on panel width
    for (k = jcol + 1; k < (std::min)(jcol + panel_size, n); k++) {
      if (relax_end(k) != emptyIdxLU) {
        panel_size = k - jcol;
        break;
      }
    }
    if (k == n) panel_size = n - jcol;
 
    // Symbolic outer factorization on a panel of columns
    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);
 
    // Numeric sup-panel updates in topological order
    Base::panel_bmod(m, panel_size, jcol, nseg1, dense, tempv, segrep, repfnz, m_glu);
 
    // Sparse LU within the panel, and below the panel diagonal
    for (jj = jcol; jj < jcol + panel_size; jj++) {
      k = (jj - jcol) * m;  // Column index for w-wide arrays
 
      nseg = nseg1;  // begin after all the panel segments
      // Depth-first-search for the current column
      VectorBlock<IndexVector> panel_lsubk(panel_lsub, k, m);
      VectorBlock<IndexVector> repfnz_k(repfnz, k, m);
      // Return 0 on success and > 0 number of bytes allocated when run out of space.
      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);
      if (info) {
        m_lastError = "UNABLE TO EXPAND MEMORY IN COLUMN_DFS() ";
        m_info = NumericalIssue;
        m_factorizationIsOk = false;
        return;
      }
      // Numeric updates to this column
      VectorBlock<ScalarVector> dense_k(dense, k, m);
      VectorBlock<IndexVector> segrep_k(segrep, nseg1, m - nseg1);
      // Return 0 on success and > 0 number of bytes allocated when run out of space.
      info = Base::column_bmod(jj, (nseg - nseg1), dense_k, tempv, segrep_k, repfnz_k, jcol, m_glu);
      if (info) {
        m_lastError = "UNABLE TO EXPAND MEMORY IN COLUMN_BMOD() ";
        m_info = NumericalIssue;
        m_factorizationIsOk = false;
        return;
      }
 
      // Copy the U-segments to ucol(*)
      // Return 0 on success and > 0 number of bytes allocated when run out of space.
      info = Base::copy_to_ucol(jj, nseg, segrep, repfnz_k, m_perm_r.indices(), dense_k, m_glu);
      if (info) {
        m_lastError = "UNABLE TO EXPAND MEMORY IN COPY_TO_UCOL() ";
        m_info = NumericalIssue;
        m_factorizationIsOk = false;
        return;
      }
 
      // Form the L-segment
      // Return 0 if success, i > 0 if U(i, i) is exactly zero.
      info = Base::pivotL(jj, m_diagpivotthresh, m_perm_r.indices(), iperm_c.indices(), pivrow, m_glu);
      if (info) {
        m_lastError = "THE MATRIX IS STRUCTURALLY SINGULAR";
#ifndef EIGEN_NO_IO
        std::ostringstream returnInfo;
        returnInfo << " ... ZERO COLUMN AT ";
        returnInfo << info;
        m_lastError += returnInfo.str();
#endif
        m_info = NumericalIssue;
        m_factorizationIsOk = false;
        return;
      }
 
      // Prune columns (0:jj-1) using column jj
      Base::pruneL(jj, m_perm_r.indices(), pivrow, nseg, segrep, repfnz_k, xprune, m_glu);
 
      // Reset repfnz for this column
      for (i = 0; i < nseg; i++) {
        irep = segrep(i);
        repfnz_k(irep) = emptyIdxLU;
      }
    }                    // end SparseLU within the panel
    jcol += panel_size;  // Move to the next panel
  }                      // end for -- end elimination
 
  m_detPermR = m_perm_r.determinant();
  m_detPermC = m_perm_c.determinant();
 
  // Count the number of nonzeros in factors
  Base::countnz(n, m_nnzL, m_nnzU, m_glu);
  // Apply permutation  to the L subscripts
  Base::fixupL(n, m_perm_r.indices(), m_glu);
 
  // Create supernode matrix L
  m_Lstore.setInfos(m, n, m_glu.lusup, m_glu.xlusup, m_glu.lsub, m_glu.xlsub, m_glu.supno, m_glu.xsup);
  // Create the column major upper sparse matrix  U;
  new (&m_Ustore) Map<SparseMatrix<Scalar, ColMajor, StorageIndex>>(m, n, m_nnzU, m_glu.xusub.data(), m_glu.usub.data(),
                                                                    m_glu.ucol.data());
 
  m_info = Success;
  m_factorizationIsOk = true;
}
 
template <typename MappedSupernodalType>
struct SparseLUMatrixLReturnType : internal::no_assignment_operator {
  using Scalar = typename MappedSupernodalType::Scalar;
  explicit SparseLUMatrixLReturnType(const MappedSupernodalType& mapL) : m_mapL(mapL) {}
  Index rows() const { return m_mapL.rows(); }
  Index cols() const { return m_mapL.cols(); }
  template <typename Dest>
  void solveInPlace(MatrixBase<Dest>& X) const {
    m_mapL.solveInPlace(X);
  }
  template <bool Conjugate, typename Dest>
  void solveTransposedInPlace(MatrixBase<Dest>& X) const {
    m_mapL.template solveTransposedInPlace<Conjugate>(X);
  }
 
  SparseMatrix<Scalar, ColMajor, Index> toSparse() const {
    ArrayXi colCount = ArrayXi::Ones(cols());
    for (Index i = 0; i < cols(); i++) {
      typename MappedSupernodalType::InnerIterator iter(m_mapL, i);
      for (; iter; ++iter) {
        if (iter.row() > iter.col()) {
          colCount(iter.col())++;
        }
      }
    }
    SparseMatrix<Scalar, ColMajor, Index> sL(rows(), cols());
    sL.reserve(colCount);
    for (Index i = 0; i < cols(); i++) {
      sL.insert(i, i) = 1.0;
      typename MappedSupernodalType::InnerIterator iter(m_mapL, i);
      for (; iter; ++iter) {
        if (iter.row() > iter.col()) {
          sL.insert(iter.row(), iter.col()) = iter.value();
        }
      }
    }
    sL.makeCompressed();
    return sL;
  }
 
  const MappedSupernodalType& m_mapL;
};
 
template <typename MatrixLType, typename MatrixUType>
struct SparseLUMatrixUReturnType : internal::no_assignment_operator {
  using Scalar = typename MatrixLType::Scalar;
  SparseLUMatrixUReturnType(const MatrixLType& mapL, const MatrixUType& mapU) : m_mapL(mapL), m_mapU(mapU) {}
  Index rows() const { return m_mapL.rows(); }
  Index cols() const { return m_mapL.cols(); }
 
  template <typename Dest>
  void solveInPlace(MatrixBase<Dest>& X) const {
    Index nrhs = X.cols();
    // Backward solve with U
    for (Index k = m_mapL.nsuper(); k >= 0; k--) {
      Index fsupc = m_mapL.supToCol()[k];
      Index lda = m_mapL.colIndexPtr()[fsupc + 1] - m_mapL.colIndexPtr()[fsupc];  // leading dimension
      Index nsupc = m_mapL.supToCol()[k + 1] - fsupc;
      Index luptr = m_mapL.colIndexPtr()[fsupc];
 
      if (nsupc == 1) {
        for (Index j = 0; j < nrhs; j++) {
          X(fsupc, j) /= m_mapL.valuePtr()[luptr];
        }
      } else {
        // FIXME: the following lines should use Block expressions and not Map!
        Map<const Matrix<Scalar, Dynamic, Dynamic, ColMajor>, 0, OuterStride<>> A(&(m_mapL.valuePtr()[luptr]), nsupc,
                                                                                  nsupc, OuterStride<>(lda));
        typename Dest::RowsBlockXpr U = X.derived().middleRows(fsupc, nsupc);
        U = A.template triangularView<Upper>().solve(U);
      }
 
      for (Index j = 0; j < nrhs; ++j) {
        for (Index jcol = fsupc; jcol < fsupc + nsupc; jcol++) {
          typename MatrixUType::InnerIterator it(m_mapU, jcol);
          for (; it; ++it) {
            Index irow = it.index();
            X(irow, j) -= X(jcol, j) * it.value();
          }
        }
      }
    }  // End For U-solve
  }
 
  template <bool Conjugate, typename Dest>
  void solveTransposedInPlace(MatrixBase<Dest>& X) const {
    using numext::conj;
    Index nrhs = X.cols();
    // Forward solve with U
    for (Index k = 0; k <= m_mapL.nsuper(); k++) {
      Index fsupc = m_mapL.supToCol()[k];
      Index lda = m_mapL.colIndexPtr()[fsupc + 1] - m_mapL.colIndexPtr()[fsupc];  // leading dimension
      Index nsupc = m_mapL.supToCol()[k + 1] - fsupc;
      Index luptr = m_mapL.colIndexPtr()[fsupc];
 
      for (Index j = 0; j < nrhs; ++j) {
        for (Index jcol = fsupc; jcol < fsupc + nsupc; jcol++) {
          typename MatrixUType::InnerIterator it(m_mapU, jcol);
          for (; it; ++it) {
            Index irow = it.index();
            X(jcol, j) -= X(irow, j) * (Conjugate ? conj(it.value()) : it.value());
          }
        }
      }
      if (nsupc == 1) {
        for (Index j = 0; j < nrhs; j++) {
          X(fsupc, j) /= (Conjugate ? conj(m_mapL.valuePtr()[luptr]) : m_mapL.valuePtr()[luptr]);
        }
      } else {
        Map<const Matrix<Scalar, Dynamic, Dynamic, ColMajor>, 0, OuterStride<>> A(&(m_mapL.valuePtr()[luptr]), nsupc,
                                                                                  nsupc, OuterStride<>(lda));
        typename Dest::RowsBlockXpr U = X.derived().middleRows(fsupc, nsupc);
        EIGEN_IF_CONSTEXPR (Conjugate)
          U = A.adjoint().template triangularView<Lower>().solve(U);
        else
          U = A.transpose().template triangularView<Lower>().solve(U);
      }
    }  // End For U-solve
  }
 
  SparseMatrix<Scalar, RowMajor, Index> toSparse() {
    ArrayXi rowCount = ArrayXi::Zero(rows());
    for (Index i = 0; i < cols(); i++) {
      typename MatrixLType::InnerIterator iter(m_mapL, i);
      for (; iter; ++iter) {
        if (iter.row() <= iter.col()) {
          rowCount(iter.row())++;
        }
      }
    }
 
    SparseMatrix<Scalar, RowMajor, Index> sU(rows(), cols());
    sU.reserve(rowCount);
    for (Index i = 0; i < cols(); i++) {
      typename MatrixLType::InnerIterator iter(m_mapL, i);
      for (; iter; ++iter) {
        if (iter.row() <= iter.col()) {
          sU.insert(iter.row(), iter.col()) = iter.value();
        }
      }
    }
    sU.makeCompressed();
    const SparseMatrix<Scalar, RowMajor, Index> u = m_mapU;  // convert to RowMajor
    sU += u;
    return sU;
  }
 
  const MatrixLType& m_mapL;
  const MatrixUType& m_mapU;
};
 
}  // End namespace Eigen
 
#endif