尖点CG收敛

Question

我使用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中设置正确？

Answer 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];