2

我在 Windows 7 x64 上使用 Microsoft Visual Studio 2008。我试图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， 的解与我用 MATLAB 计算的解Ax=b不匹配。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吗？

4

1 回答 1

2

之前从未使用过csparse，我浏览了源代码。

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

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

顺便说一句，您的代码会泄漏内存，因为压缩矩阵是一个新分配，并且原始未压缩矩阵没有被cs_compress().

于 2013-07-30T03:10:58.420 回答