Решение простой разреженной линейной системы уравнений с использованием csparse: cs_cholsol

Я использую Microsoft Visual Studio 2008 на Windows 7 x64. Я пытаюсь решить следующую линейную систему Ax=b используя csparse, где A положительно определен

    | 1  0  0  1 |
A = | 0  3  1  0 |
| 0  1  2  1 |
| 1  0  1  2 |

| 1 |
b = | 1 |
| 1 |
| 1 |

Я использовал следующие коды

int Ncols = 4, Nrows = 4, nnz = 10;
int cols[]    = {0, 3, 1, 2, 1, 2, 3, 0, 2, 3};
int rows[]    = {0, 0, 1, 1, 2, 2, 2, 3, 3, 3};
double vals[] = {1, 1, 3, 1, 1, 2, 1, 1, 1, 2};

cs *Operator = cs_spalloc(Ncols,Nrows,nnz,1,1);

int j;
for(j = 0; j < nnz; j++)
{
Operator->i[j] = rows[j];
Operator->p[j] = cols[j];
Operator->x[j] = vals[j];
Operator->nz++;
}

for(j = 0; j < nnz; j++)
cout << Operator->i[j] << " " << Operator->p[j] << " " << Operator->x[j] << endl;

Operator = cs_compress(Operator);

for(j = 0; j < nnz; j++)
cout << Operator->i[j] << " " << Operator->p[j] << " " << Operator->x[j] << endl;

// Right hand side
double b[] = {1, 1, 1, 1};

// Solving Ax = b
int status = cs_cholsol(0, Operator, &b[0]); // status = 0 means error.

Чтобы убедиться, что я правильно создал разреженную переменную, я попытался распечатать индекс строк и столбцов, а также их значения в консоли до и после cs_compress, Ниже приведен результат этой распечатки.

До:

0 0 1
0 3 1
1 1 3
1 2 1
2 1 1
2 2 2
2 3 1
3 0 1
3 2 1
3 3 2

После:

0 0 1
3 2 1
1 4 3
2 7 1
1 10 1
2 -6076574517017313795 2
3 -6076574518398440533 1
0 -76843842582893653 1
2 0 1
3 0 2

Из-за значений мусора, которые можно наблюдать выше после вызова cs_compressрешение Ax=b не совпадает с тем, что я рассчитал с MATLAB. MATLAB приводит к следующему решению.

    | 2.0000 |
x = | 0.0000 |
| 1.0000 |
|-1.0000 |

Интересно, что у меня нет этой проблемы для следующих кодов, которые решают Ax=b, где A это 3×3 единичная матрица.

int Ncols = 3, Nrows = 3, nnz = Nrows;

cs *Operator = cs_spalloc(Ncols,Nrows,nnz,1,1);
int j;
for(j = 0; j < nnz; j++) {
Operator->i[j] = j;
Operator->p[j] = j;
Operator->x[j] = 1.0;
Operator->nz++;
}

Operator = cs_compress(Operator);

double b[] = {1, 2, 3};

int status = cs_cholsol(0, Operator, &b[0]); // status = 1 means no error.

Может кто-нибудь, пожалуйста, помогите мне решить проблему с cs_compress?