Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Out of memory (or invalid) with small matrix #850

Open
kmfrick opened this issue Jul 5, 2024 · 0 comments
Open

Out of memory (or invalid) with small matrix #850

kmfrick opened this issue Jul 5, 2024 · 0 comments
Assignees

Comments

@kmfrick
Copy link

kmfrick commented Jul 5, 2024

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;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants