SuiteSparse icon indicating copy to clipboard operation
SuiteSparse copied to clipboard

Out of memory (or invalid) with small matrix

Open kmfrick opened this issue 7 months ago • 0 comments

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

  1. Compile the code at the bottom.
  2. Run it, reading A and B from the attached CSV files. Yj.csv Aj.csv

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;
}

kmfrick avatar Jul 05 '24 18:07 kmfrick