cuSolver没有返回正确的解决方案

Question

我试图在cuSOLVER中使用QR线性系统求解器，如下

#include <cusparse_v2.h>
#include <stdio.h>
#include <cuda.h>
#include <cuda_runtime.h>
#include "device_launch_parameters.h"
#include <iostream>
#include <cassert>

#include <cublas_v2.h>
#include <cusolverDn.h>
#include <cusolverSp.h>

#include <thrust/host_vector.h>
#include <thrust/device_vector.h>

using namespace thrust;

#define CUSPARSE_CHECK(x) {cusparseStatus_t _c=x; if (_c != CUSPARSE_STATUS_SUCCESS) {printf("cusparse fail: %d, line: %d\n", (int)_c, __LINE__); exit(-1);}}

void init_handlers_and_matr_descriptor(cusolverSpHandle_t& cusolverH,cusparseHandle_t& cusparseH,cusparseMatDescr_t& descrA) {
    assert(CUSOLVER_STATUS_SUCCESS == cusolverSpCreate(&cusolverH));
    assert(CUSPARSE_STATUS_SUCCESS == cusparseCreate(&cusparseH));
    assert(CUSPARSE_STATUS_SUCCESS == cusparseCreateMatDescr(&descrA)); //CUSPARSE_INDEX_ZERO,CUSPARSE_MATRIX_TYPE_GENERAL
    assert(CUSPARSE_STATUS_SUCCESS == cusparseSetMatIndexBase(descrA,CUSPARSE_INDEX_BASE_ZERO));
    assert(CUSPARSE_STATUS_SUCCESS == cusparseSetMatType(descrA, CUSPARSE_MATRIX_TYPE_GENERAL));
}

int sparse_solver_test() {

    //Init csr format A and b for solving Ax = b
    /*
    A =
    [  1.0  2.0  0.0  0.0  0.0  0.0 ]
    [  3.0  4.0  5.0  0.0  0.0  0.0 ]
    [  0.0  6.0  7.0  8.0  0.0  0.0 ]
    [  0.0  0.0  9.0 10.0 11.0  0.0 ]
    [  0.0  0.0  0.0 12.0 13.0 14.0 ]
    [  0.0  0.0  0.0  0.0 15.0 16.0 ]

    b = [0.0,2.0,4.0,6.0,8.0,10.0]^T
    */

    int nnz_A = 16, m = 6;

    cusolverSpHandle_t cusolverH = NULL;
    cusparseHandle_t cusparseH = NULL; //cuBLAS or cuSPARSE?
    cusolverStatus_t cusolver_status = CUSOLVER_STATUS_SUCCESS;
    cusparseMatDescr_t descrA;
    cudaError_t cudaStat;

    init_handlers_and_matr_descriptor(cusolverH, cusparseH, descrA);

    host_vector<double> h_csrValA(nnz_A), h_dataB(m), h_dataX(m);
    host_vector<int> h_csrRowPtrA(m + 1), h_csrColIndA(nnz_A);
    device_vector<double> d_csrValA(nnz_A), d_dataB(m), d_dataX(m);
    device_vector<int> d_csrRowPtrA(m + 1), d_csrColIndA(nnz_A);

    //Init matrix
    for (auto i = 0; i < nnz_A; ++i) h_csrValA[i] = static_cast<double>(i + 1);
    h_csrRowPtrA[0] = 0;
    h_csrRowPtrA[1] = 2;
    for (auto i = 2; i < m; ++i) h_csrRowPtrA[i] = h_csrRowPtrA[i - 1] + 3;
    h_csrRowPtrA[m] = nnz_A;

    h_csrColIndA[0] = 0;
    int v[] = {0, 1, 0, 1, 2, 1, 2, 3, 2, 3, 4, 3, 4, 5, 4, 5 };
    for (auto i = 0; i < nnz_A; ++i) h_csrColIndA[i] = v[i];

    for (auto i = 0; i < m; ++i) h_dataB[i] = static_cast<double>(2 * i);

    //device memory and descriptor A init
    d_csrValA = h_csrValA;
    d_csrRowPtrA = h_csrRowPtrA;
    d_csrColIndA = h_csrColIndA;
    d_dataB = h_dataB;

    //step4, solve the linear system?
    int singularity;
    cusolver_status = cusolverSpDcsrlsvqr(
        cusolverH, m, nnz_A, descrA,
        d_csrValA.data().get(), d_csrRowPtrA.data().get(), d_csrColIndA.data().get(), d_dataB.data().get(),
        0.0, 0, d_dataX.data().get(), &singularity);

    std::cout << "singularity = " << singularity << std::endl;
    assert(CUSOLVER_STATUS_SUCCESS == cusolver_status);

    h_dataX = d_dataX;

    std::cout << "x = (";
    for (auto i = 0; i < m; ++i) {
        std::cout << h_dataX[i];
        if (i < m - 1) std::cout << ", ";
    }
    std::cout << std::endl;

    if (cusparseH)
        cusparseDestroy(cusparseH);
    if (cusolverH)
        cusolverSpDestroy(cusolverH);
    if (descrA)
        cusparseDestroyMatDescr(descrA);
    cudaDeviceReset();
    return 0;
}

int main(int argc, char** argv) {
    sparse_solver_test();
    return 0;
}

不确定我的函数设置是否错误，有人能帮忙吗？

更新我使用推力库对代码进行了一些简化，尽管错误仍然是一样的，但是至少我去掉了所有的malloc等等.

Update按照建议更正了csrIndColA (相应地更改了代码)数组。现在，求解器可以工作了(也就是说，我不再得到之前得到的错误了)，尽管结果是0。

更新--我所做的所有更改--也忘记了初始化h_dataB，以及解决问题的csrIndColA索引，完整的代码在上面供将来参考。

Answer 1

在您的示例中，csrColIndA数组太短了，因此cuSOLVER试图从它的末尾读取。

根据cuSOLVER文档和一般约定，列索引数组的长度与非零矩阵条目数组的长度相同，并存储每个非零元素的列索引(而不仅仅是每个列中的第一个非零元素，这将限制格式为所有非零元素垂直连续的稀疏模式)。

所以您的示例输出应该有

csrColIndA = {0, 1, 0, 1, 2, 1, 2, 3, 2, 3, 4, 3, 4, 5, 4, 5}