opencv — не удалось использовать cs_qrsol из CXSparse для решения x = A \ b в C ++, когда матрица A большая

Я пытаюсь решить систему линейных уравнений x = A \ b, используя библиотеку CXSparse Тимом Дэвисом (http://faculty.cse.tamu.edu/davis/suitesparse.html). Я разрабатываю свою программу C ++ (с OpenCV), используя MS Visual Studio 2012 на Windows 7 x64.

У меня есть матрица A, которая является разреженной матрицей, и матрица b. Я хочу найти х = А \ б. Моя матрица А имеет размер ~ 700 000 х 30 000. Поскольку это не квадратная матрица, я вызываю функцию cs_qrsol (так как эта функция принимает неквадратную матрицу mxn.

Я попытался cs_qrsol решить небольшую проблему, используя мой код ниже, и это дает ожидаемый результат. Однако, когда я пытался применить его с моей реальной матрицей A, функция всегда возвращает 0. (Я пробовал order = от 0 до 3).

Чтобы доказать, что у моей задачи A, b есть решение x, я ввожу матрицу A, b в MATLAB и могу успешно получить результаты.

Итак, мне было интересно, что может быть что-то не так с моим кодом или есть какие-либо ограничения на функцию ??. Я был бы очень признателен, если бы кто-нибудь мог помочь указать, что я сделал не так. Большое спасибо.

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