I use CUSP conjugate gradient method to solve my symmetric sparse matrix. And I have no idea why it doesn't converge. Dimensions of matrices I use are not that large (1K to 100K). The same linear systems are easily solved by MKL, so the matrix is not ill-conditioned. However I tried adding preconditioner, but it gave no results:
diagonal preconditioner and AINV (incomplete Cholesky) gave unlimited growth of residual (as long as cg and bicgstab)
Here's my code:
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];
Why isn't it working correctly?
P.S. As I understand CUSP supposes zero-based index, that's why I decrease csrRowPtr and csrColInd. When I used nVidia cuSparse library there was an option to set other parameters like matrix type and fill mode. How can I be sure they are set correctly in CUSP?
Only elements from the upper triangle are stored in MKL's CSR format, but all the elements must be stored in CUSP's CSR format even if you are solving a symmetric linear system.
Also I think
is not a good idea; first transfer it back to host_memory
then copy it back to x_host or whatever your result array is