使用csparse求解一个简单的稀疏线性方程组:cs_cholsol

我在Windows 7 x64上使用Microsoft Visual Studio 2008。我试图用csparse求解下面的线性系统Ax=b，其中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的问题吗？