LU-разложение на разреженную прямоугольную матрицу в R

Решено. Следующий код использует UMFPACK, ту же библиотеку, которую использует MATLAB. Чтобы получить нужные элементы, нам нужно добавить Control[UMFPACK_SCALE] = 0;. Затем нам нужно зафиксировать размер полученной матрицы U в R.

#include <suitesparse/umfpack.h>
#include <Rcpp.h>
#include <algorithm>
#include <vector>
#include <iostream>

// [[Rcpp::export]]
Rcpp::List sparseLU(const std::vector<int> &Ap, const std::vector<int> &Ai, const std::vector<double> &Ax) {
  int n = Ap.size() - 1;
  int m = *std::max_element(Ai.begin(), Ai.end()) + 1;
  int nnz = Ax.size();
  double Control [UMFPACK_CONTROL], Info [UMFPACK_INFO];
  umfpack_di_defaults(Control);
  Control[UMFPACK_PRL] = 2;
  Control[UMFPACK_SCALE] = 0;
  Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_SYMMETRIC;

  void *symbolic, *numeric;
  int status;
  status = umfpack_di_symbolic(m, n, Ap.data(), Ai.data(), Ax.data(), &symbolic, Control, Info);
  if (status != UMFPACK_OK) {
    umfpack_di_report_status(Control, status);
    if (status < 0) {
      Rcpp::stop("umfpack_di_symbolic failed with status %d", status);
    }
  }

  status = umfpack_di_numeric(Ap.data(), Ai.data(), Ax.data(), symbolic, &numeric, Control, Info);
  if (status != UMFPACK_OK) {
    umfpack_di_report_status(Control, status);
    if (status < 0) {
      Rcpp::stop("umfpack_di_numeric failed with status %d", status);
      umfpack_di_free_symbolic(&symbolic);
    }
  }

  umfpack_di_free_symbolic(&symbolic);

  int lnz, unz, nz_udiag;
  status = umfpack_di_get_lunz(&lnz, &unz, &m, &n, &nz_udiag, numeric);
  if (status != UMFPACK_OK) {
    umfpack_di_report_status(Control, status);
    if (status < 0) {
      Rcpp::stop("umfpack_di_get_lunz failed with status %d", status);
      umfpack_di_free_numeric(&numeric);
    }
  }


  std::vector<int> Lp(m+1), Lj(lnz), Up(n+1), Ui(unz), P(m), Q(n);
  std::vector<double> Lx(lnz), Ux(unz), D(std::min(m, n)), Rs(m);
  int do_recip;

  status = umfpack_di_get_numeric(Lp.data(), Lj.data(), Lx.data(), Up.data(), Ui.data(), Ux.data(), P.data(), Q.data(), D.data(), &do_recip, Rs.data(), numeric);
  if (status != UMFPACK_OK) {
    umfpack_di_report_status(Control, status);
    if (status < 0) {
      umfpack_di_free_numeric(&numeric);
      Rcpp::stop("umfpack_di_get_numeric failed with status %d", status);
    }
  }

  umfpack_di_free_numeric(&numeric);

  Rcpp::List ret = Rcpp::List::create(
      Rcpp::Named("L") = Rcpp::List::create(Rcpp::Named("j") = Lj, Rcpp::Named("p") = Lp, Rcpp::Named("x") = Lx),
      Rcpp::Named("U") = Rcpp::List::create(Rcpp::Named("i") = Ui, Rcpp::Named("p") = Up, Rcpp::Named("x") = Ux),
      Rcpp::Named("P") = P,
      Rcpp::Named("Q") = Q
      );


  return ret;
}

Затем мы можем загрузить его в R вот так

# A is a sparseMatrix() from package Matrix
LUPQ <- sparseLU(A@p, A@i, A@x) 
n <- ncol(A)
L <- sparseMatrix(j = LUPQ$L$j, p = LUPQ$L$p, x = LUPQ$L$x, index1 = FALSE) 
U <- sparseMatrix(i = LUPQ$U$i, p = LUPQ$U$p, x = LUPQ$U$x, index1 = FALSE, dims = c(n, n))
Q <- sparseMatrix(i = LUPQ$Q + 1, j = 1:length(LUPQ$Q), x = 1)