CUSP CG конвергенция

Question

Я использую метод сопряженных градиентов CUSP для решения моей симметричной разреженной матрицы. И я понятия не имею, почему это не сходится. Размеры используемых мной матриц не так велики (от 1К до 100К). Такие же линейные системы легко решаются с помощью 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];

Почему это не работает правильно?

PS Как я понимаю, CUSP предполагает нулевой индекс, поэтому я уменьшаю csrRowPtr и csrColInd. Когда я использовал библиотеку nVidia cuSparse, была возможность установить другие параметры, такие как тип матрицы и режим заполнения. Как я могу быть уверен, что они правильно установлены в CUSP?

Answer 1

Только элементы из верхнего треугольника хранятся в формате CSR MKL, но все элементы должны храниться в формате CSR CUSP, даже если вы решаете симметричную линейную систему.

Также я думаю

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