c++ - 使用 csparse : cs_cholsol 求解简单的稀疏线性方程组

Question

我在 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 的问题吗？

Answer 1

从未与 csparse 合作过之前，我浏览了 source code .

当您调用 cs_spalloc() 时要创建 Operator，您将创建一个三元组(通过将最后一个参数设置为 1 来表示)。但是，在调用 cs_copmress() 之后，结果不再是三元组(您可以通过检查结果检测到这一点，并看到压缩后 Operator->n 现在是 -1)。所以，像这样遍历矩阵是错误的。

您可以使用 cs_print()用于打印稀疏矩阵的 API。

顺便说一句，您的代码会泄漏内存，因为压缩矩阵是新分配的，而原始未压缩矩阵未被 cs_compress() 释放。