Search code examples
c++opencvsparse-matrixleast-squaressuitesparse

Failed to use cs_qrsol from CXSparse to solve x=A\b in C++ when A matrix is large


I am trying to solve a linear equation system x = A\b, using CXSparse library by Tim Davis (http://faculty.cse.tamu.edu/davis/suitesparse.html). I develop my C++ program (with OpenCV) using MS Visual Studio 2012 on Windows 7 x64.

I have an A matrix, which is a sparse matrix and b matrix. I want to find x = A\b. My A matrix is of size ~ 700,000 x 30,000. Since this is not a square matrix, I call the function cs_qrsol (as this function accept a non-square mxn matrix.

I have tried cs_qrsol to solve a small problem using my code below, and it gives an expected result. However, when I tried to apply it with my real A matrix, the function always returns 0. (I have tried order = 0 to 3).

To proof that my problem A, b has a solution x, I input the matrix A,b in MATLAB and I can get the results successfully.

So, I was wondering there might be something wrong with my code or there is any limitation on the function??. I would really appreciate if anybody could help point out what I did wrong. Thank you very much.

bool sparseSolve(cv::Mat& i, cv::Mat& j, cv::Mat& s, cv::Mat& d, int k, int n)
{
// Solve linear equations: x = A\b
// i,j,s are (k x 1) matrices - (i,j) are row and col, s are values so the trippet is (i,j,s)
// d is an (n x 1) matrix - this is the b matrix

cs *T, *A;
double * b;
int y, m, status;

// Assign the sparse matrix A
T = cs_spalloc (0, 0, 1, 1, 1) ; 
for (y = 0; y < k; y++)
{
    if (!cs_entry (T, i.at<int>(y,0), j.at<int>(y,0), s.at<double>(y,0)))
    {
        cs_spfree(T);
        std::cout << "Failed to add entry to the matrix A at " << y << std::endl;
        return false;
    }
}

A = cs_compress(T);
cs_spfree(T);
if (!A)
{
    std::cout << "Failed to create A as the compressed-column form of T" << std::endl;
    return false;   
}
m = A->m;   // # rows of A

if (n != m)
    std::cout << "# rows of A (" << m << ") not equal # equations (" << n << ")" << std::endl;

// Allocate the b matrix
b = static_cast<double *>(malloc(n*sizeof(double)));
if (!b)
{
    std::cout << "Failed to allocate the b matrix" << std::cout;
    cs_spfree(A);
    return false;
}

// Assign the b matrix
for (y = 0; y < n; y++)
{
    b[y] = d.at<double>(y,0);
}

// Obtain the results x=A\b in the b matrix
status = cs_cholsol(0, A, b); // This returns 0 as failed.

cs_spfree(A);
free(b);

return (1==status);
}

Solution

  • Now, I can solve my own problem. I made change to my program to call the version cs_dl_* instead of cs_* functions. Also, the memory in my machine is really relevant. To make it work, I have to close all the opened application to make sure that I have enough space of memory for cs_spalloc or cs_dl_spalloc.