SuiteSparse
SuiteSparse copied to clipboard
Out of memory (or invalid) with small matrix
Describe the bug
I'm trying to reproduce what A\b
does in MATLAB using SQPR. It goes out of memory but the matrix is relatively small ( 17261-by-5547, nz 17261) and MATLAB can process it without issues. In the comments for SuiteSparseQR_min2norm
I find that
// After A and B are valid, the only way any of the above functions can
// fail is by running out of memory. However, if (say) AT is NULL
// because of insufficient memory, SuiteSparseQR_factorize will report
// CHOLMOD_INVALID since its input is NULL. In that case, restore the
// error back to CHOLMOD_OUT_OF_MEMORY.
Indeed, if I use SuiteSparseQR, I get status -4 (INVALID).
Please let me know if there's any more information you may need.
To Reproduce
Expected behavior I get a vector with 4328 non-zero elements whose first 10 elements are: 29.6637 9.8166 0 25.4857 46.3880 17.5530 5.8954 0 4.9303 26.3932
Observed behavior Returns status -2.
Desktop (please complete the following information):
- OS: macOS 14.5 (23F79)
- compiler: Apple clang version 15.0.0 (clang-1500.3.9.4) Target: x86_64-apple-darwin23.5.0 Thread model: posix
- Version: 7.7.0 installed with Homebrew
#include <suitesparse/SuiteSparseQR.hpp>
#include <vector>
#include <iostream>
#include <fstream>
int main() {
int m = 17261, n = 5547, nz = 17262;
std::vector<int> Ai;
std::vector<int> Aj;
std::vector<double> Ax;
std::ifstream fileA("Aj.csv");
std::vector<double> b;
std::ifstream fileB("Yj.csv");
for (int i = 0; i < nz; i++) {
int a, b;
double c;
fileA >> a;
fileA >> b;
fileA >> c;
Ai.emplace_back(a - 1);
Aj.emplace_back(b - 1);
Ax.emplace_back(c);
}
for (int i = 0; i < nz; i++) {
std::cout << Ai[i] << " " << Aj[i] << " " << Ax[i] << std::endl;
}
for (int i = 0; i < m; i++) {
double a;
fileB >> a;
b.emplace_back(a);
}
cholmod_common c;
cholmod_start(&c);
c.print = 4;
cholmod_triplet *T = cholmod_allocate_triplet(m, n, nz, 0, CHOLMOD_REAL, &c);
for (int i = 0; i < nz; i++) {
static_cast<int*>(T->i)[i] = Ai[i];
static_cast<int*>(T->j)[i] = Aj[i];
static_cast<double*>(T->x)[i] = Ax[i];
}
T->nnz = nz;
cholmod_sparse *A = cholmod_triplet_to_sparse(T, nz, &c);
cholmod_dense *B, *X;
int status_check = cholmod_print_sparse(A, "Aj", &c);
// Copy vector b to cholmod_dense format
B = cholmod_allocate_dense(b.size(), 1, b.size(), CHOLMOD_REAL, &c);
for (int i = 0; i < b.size(); i++) {
static_cast<double*>(B->x)[i] = b[i];
}
status_check = cholmod_print_dense(B, "Yj", &c);
// Solve Ax = b using SPQR
X = SuiteSparseQR_min2norm<double>(SPQR_ORDERING_DEFAULT, SPQR_DEFAULT_TOL, A, B, &c);
if (c.status < 0) {
std::cerr << "SPQR solve failed with status" << c.status << std::endl;
}
cholmod_print_dense(X, "X", &c);
// Extract solution
std::vector<double> x(m);
for (int i = 0; i < m; i++) {
x[i] = static_cast<double*>(X->x)[i];
}
cholmod_free_dense(&B, &c);
cholmod_free_dense(&X, &c);
cholmod_free_sparse(&A, &c);
cholmod_finish(&c);
for (auto x_val : x) {
std::cout << x_val << std::endl;
}
return 0;
}