cuSolver 未返回正确的解

cuSolver doesn't return the correct solution

本文关键字:返回 cuSolver      更新时间:2023-10-16

我正在尝试在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: %dn", (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 等......

更新 按照建议更正了csrIndColA(相应地更改了代码(数组。现在求解器工作了(即我不再得到我之前得到的错误(,尽管我得到的结果仍然是 0。

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

示例中

csrColIndA数组太短,因此 cuSOLVER 会尝试读取它的末尾。

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

所以你的示例输出应该有

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