Я пытаюсь использовать SuiteSparse SPQR для решения системы линейных уравнений x = A \ b; Моя матрица A разрежена, и это прямоугольная матрица, поэтому я выбрал SPQR для решения этой проблемы.
Я построил SuiteSparse, используя MS Visual Studio 2012 на Windows 7 x64, используя те, которые предоставляются https://github.com/jlblancoc/suitesparse-metis-for-windows.
Чтобы протестировать функцию, я изменил проект spqr_example, чтобы выделить трипеты перед преобразованием в разреженную матрицу, вместо первоначального чтения входных данных из stdin для создания разреженной матрицы. Я ввел небольшую матрицу А и В для тестирования. Программа успешно скомпилирована. Я отладил программу и обнаружил, что мой вызов cholmod_allocate_triplet () не удался, потому что в объявлении этой функции он имеет следующий код:
RETURN_IF_NULL_COMMON (NULL) ;
Это всегда возвращает ложь (даже если моя общая начинается успешно).
Я не хочу явно вносить изменения в эту строку, поскольку я мог где-то ошибиться или я забыл сделать что-то, что мне нужно сделать, потому что я новичок в использовании библиотеки.
Кто-нибудь может помочь дать мне несколько советов о том, как правильно запустить мою программу? Мой код ниже изменен из предоставленного spqr_example. Большое спасибо.
#include <iostream>
#include "SuiteSparseQR.hpp"
int main (int argc, char **argv)
{
cholmod_common Common, *cc ;
cholmod_sparse *A ;
cholmod_dense *X, *B, *Residual ;
double rnorm, one [2] = {1,0}, minusone [2] = {-1,0} ;
int mtype ;
// start CHOLMOD
cc = &Common ;
cholmod_l_start (cc) ;
// load A
//A = (cholmod_sparse *) cholmod_l_read_matrix (stdin, 1, &mtype, cc) ;
// A = [ 1 0 0 0;
// -1 1 0 0; ...
// 0 -1 1 0; ...
// 0 0 -1 1; ...
// 0 0 0 -1];
int row[] = {0, 1, 1, 2, 2, 3, 3, 4};
int col[] = {0, 0, 1, 1, 2, 2, 3, 3};
double val[] = {1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0};
int numEq = 5;
int numElement = 8;
int numSol = 4;
double b[] = {5.0, -5.0, 2.0, 1.0, 0.0};
cholmod_triplet* triplet = cholmod_allocate_triplet(5,4,5*4,0,CHOLMOD_REAL,cc);
int * triplet_i = (int *)(triplet->i);
int * triplet_j = (int *)(triplet->j);
double * triplet_x = (double *)(triplet->x);
for (int ne=0; ne<numElement; ne++)
{
triplet_i[triplet->nnz] = row[ne];
triplet_j[triplet->nnz] = col[ne];
triplet_x[triplet->nnz] = val[ne];
triplet->nnz++;
}
// Convert triplet to sparse matrix
A = cholmod_triplet_to_sparse(triplet, numElement, cc);
cholmod_free_triplet(&triplet, cc);
// B = ones (size (A,1),1)
//B = cholmod_l_ones (A->nrow, 1, A->xtype, cc) ;
B = cholmod_l_zeros(numEq, 1, CHOLMOD_REAL, cc);
for (int ne=0; ne<numEq; ne++)
{
((double *)(B->x))[ne] = val[ne];
}
// X = A\B
X = SuiteSparseQR<double>(A,B,cc);
//X = SuiteSparseQR <double> (A, B, cc) ;
// Print out the result
double *sol = static_cast<double *>(malloc(sizeof(X->x)));
sol = (double *)(X->x);
for (int r=0; r<numSol; r++)
{
std::cout << "x[" << r << "] = " << sol << std::endl;
sol++;
}
///// END HERE
// rnorm = norm (B-A*X)
Residual = cholmod_l_copy_dense (B, cc) ;
cholmod_l_sdmult (A, 0, minusone, one, X, Residual, cc) ;
rnorm = cholmod_l_norm_dense (Residual, 2, cc) ;
printf ("2-norm of residual: %8.1e\n", rnorm) ;
printf ("rank %ld\n", cc->SPQR_istat [4]) ;
// free everything and finish CHOLMOD
cholmod_l_free_dense (&Residual, cc) ;
cholmod_l_free_sparse (&A, cc) ;
cholmod_l_free_dense (&X, cc) ;
cholmod_l_free_dense (&B, cc) ;
cholmod_l_finish (cc) ;
return (0) ;
}
Я наконец узнал, почему моя программа сломалась после строки ниже
cholmod_triplet * triplet = cholmod_allocate_triplet (5,4,5 * 4,0, CHOLMOD_REAL, cc);
в результате внутренний вызов cholmod_allocate_triplet () вызывает RETURN_IF_NULL_COMMON (NULL), которые возвращают false.
Причина в том, что я запускаю процесс вызова
cholmod_l_start (cc);
которая является длинной int-версией cholmod_start ().
Чтобы решить эту проблему, я должен вызвать cholmod_l_allocate_triplet () вместо cholmod_allocate_triplet (), а также изменить все другие функции, чтобы использовать cholmod_l вместо вызова только cholmod_
Других решений пока нет …