CUSP CG收敛答案

【问题标题】：CUSP CG convergenceCUSP CG收敛
【发布时间】：2014-03-06 19:42:22
【问题描述】：

我使用 CUSP 共轭梯度法来求解我的对称稀疏矩阵。而且我不知道为什么它不收敛。我使用的矩阵尺寸并不大（1K 到 100K）。相同的线性系统很容易用 MKL 求解，所以矩阵不是病态的。但是我尝试添加预处理器，但没有结果：

对角线预处理器和 AINV（不完全 Cholesky）给出了残差的无限增长（只要 cg 和 bicgstab）

这是我的代码：

cusp::csr_matrix <int, float, cusp::device_memory> A (n, n, nnz);

for (i = 0; i < n + 1; i++)
    A.row_offsets[i] = csrRowPtr[i] - 1;
for (i = 0; i < nnz; i++)
    A.values[i] = csrVal[i];
for (i = 0; i < nnz; i++)
    A.column_indices[i] = csrColInd[i] - 1;

cusp::array1d <float, cusp::device_memory> x (A.num_rows, 0);
cusp::array1d <float, cusp::device_memory> b (A.num_rows, 1);

for (i = 0; i < n; i++)
    b[i] = b_host[i];

cusp::verbose_monitor<float> monitor(b, 100, 1e-3);
cusp::identity_operator<float, MemorySpace> M(A.num_rows, A.num_rows);
    /*
    cusp::precond::diagonal<float, MemorySpace> M(A);
    cusp::precond::scaled_bridson_ainv<float, MemorySpace> M(A, .1);
    */
cusp::krylov::cg(A, x, b, monitor, M);

for (i = 0; i < n; i++)
    x_host[i] = x[i];

为什么不能正常工作？

附：据我了解，CUSP 假设索引从零开始，这就是我减少 csrRowPtr 和 csrColInd 的原因。当我使用 nVidia cuSparse 库时，有一个选项可以设置其他参数，例如矩阵类型和填充模式。我如何确定它们在 CUSP 中设置正确？

【问题讨论】：

MKL中什么方法成功解决了系统？
预处理CG。但我也用来自 CUDA SDK 的样本共轭梯度解决了这些系统（40K 以下）。
你确定你没有做错什么吗？由于您使用的是 MKL 和 CUDA SDK 并且两者都收敛，因此我假设您知道解决方案。如果您从已有的解决方案点开始使用 cusp 的 cg 会发生什么（我不是 cusp 用户，所以我不知道这是否可能）？尖点是否卡在起点，这是人们应该期待的？这是我通常用来检查优化算法中的错误的一致性测试。
这段代码似乎很大程度上来自cusp sample code。是否有某些原因您不能提供显示问题的完整可编译示例代码？您的建议似乎是 cg 求解器在您尝试的各种示例问题上都失败了，那么 cusp 示例或 CG SDK 示例中给出的合成案例怎么样？您在某处有 MemorySpace 的 typedef 吗？请提供完整的示例，例如 cusp 或 SDK CG 示例。
我怀疑您的A 矩阵组件存在问题。尝试使用cusp/verify.h 中提供的 CSR 矩阵验证。这是一个complete worked example，显示了使用cuda SDK 示例中的例程组装A 矩阵。对我来说，它似乎很快收敛。

标签： matrix cuda linear-algebra cusp-library

【解决方案1】：

只有上三角形的元素以 MKL 的 CSR 格式存储，但所有元素都必须以 CUSP 的 CSR 格式存储，即使您正在求解对称线性系统。

我也觉得

for (i = 0; i < n; i++)
    x_host[i] = x[i];

不是一个好主意；先传回host_memory

cusp::array1d<float, cusp::host_memory> _x = x;

然后将其复制回 x_host 或您的结果数组是什么

for (i = 0; i < n; i++)
    x_host[i] = _x[i];

【讨论】：