测试CUDA 11 cusolverDnDSgels

问题描述

试图理解cusolverDnDSgels功能。如果我像在文档中那样使用简单的3x3示例运行它，那么它可以工作，但是当我用我的数据运行它时，d_info返回-1，正如文档所说，如果d_info = -i则第i个参数无效。

下面我发布了3×3和4×3矩阵的代码，前者起作用而第二个不起作用。

作为参考，我使用了该网站计算器https://adrianstoll.com/linear-algebra/least-squares.html

#include <stdio.h>
#include <stdlib.h>
#include <assert.h>

#include <cuda_runtime.h>
#include <cusolverDn.h>


void printMatrix(int m,int n,const double* A,int lda,const char* name)
{
    for (int row = 0; row < m; row++) {
        for (int col = 0; col < n; coL++) {
            double Areg = A[row + col * lda];
            printf("%s(%d,%d) = %f\n",name,row + 1,col + 1,Areg);
        }
    }
}

int main(int argc,char*argv[])
{
    // 3x3 example works fine
    int m = 3;
    int n = 3;
    double A[9] = { 1.0,4.0,2.0,5.0,1.0,3.0,6.0,1.0 };
    double B[3] = { 6.0,15.0,4.0 };
    
    // 4x3 example d_info/info_gpu returns -1
    //int m = 4;
    //int n = 3;
    //double A[12] = { 1.0,2.0 };
    //double B[4] = { 6.0,5.0 };
    
    double X[3];
    
    int lda = m;
    int ldb = m;
    int ldx = n;
    int nrhs = 1;
    int niter = 0;
    int info_gpu = 0;
    size_t lwork = 0;
    
    double *d_A = NULL;
    double *d_B = NULL;
    double *d_X = NULL;
    double *d_work = NULL;
    int* d_info = NULL;
    
    cusolverDnHandle_t cusolverH = NULL;
    cudaError_t cudaStat = cudaSuccess;
    cusolverStatus_t cusolver_status = CUSOLVER_STATUS_SUCCESS;
    
    cusolver_status = cusolverDnCreate(&cusolverH);
    assert(CUSOLVER_STATUS_SUCCESS == cusolver_status);
    
    // Allocate space in the GPU
    cudaStat = cudamalloc((void**)&d_A,sizeof(double) * m * n);
    assert(cudaSuccess == cudaStat);
    
    cudaStat = cudamalloc((void**)&d_B,sizeof(double) * m * nrhs);
    assert(cudaSuccess == cudaStat);
    
    cudaStat = cudamalloc((void**)&d_X,sizeof(double) * n * nrhs);
    assert(cudaSuccess == cudaStat);
    
    cudaStat = cudamalloc((void**)&d_info,sizeof(int));
    assert(cudaSuccess == cudaStat);
    
    // copy matrices into GPU space
    cudaStat = cudamemcpy(d_A,A,sizeof(double) * m * n,cudamemcpyHostToDevice);
    assert(cudaSuccess == cudaStat);
    cudaStat = cudamemcpy(d_B,B,sizeof(double) * m * nrhs,cudamemcpyHostToDevice);
    assert(cudaSuccess == cudaStat);
    
    // Get work buffer size
    cusolver_status = cusolverDnDSgels_bufferSize(cusolverH,m,n,nrhs,d_A,lda,d_B,ldb,d_X,ldx,d_work,&lwork);
    assert(CUSOLVER_STATUS_SUCCESS == cusolver_status);
    
    // Allocate workspace
    cudaStat = cudamalloc((void**)&d_work,sizeof(float) * lwork);
    assert(cudaSuccess == cudaStat);
    
    // Run solver
    cusolver_status = cusolverDnDSgels(cusolverH,lwork,&niter,d_info);
    
    // Sync threads
    cudaStat = cudaDeviceSynchronize();
    assert(cudaSuccess == cudaStat);
    
    // copy GPU info
    cudaStat = cudamemcpy(&info_gpu,d_info,sizeof(int),cudamemcpyDevicetoHost);
    assert(cudaSuccess == cudaStat);
    
    // Get solved data
    cudaStat = cudamemcpy(X,sizeof(double) * n * nrhs,cudamemcpyDevicetoHost);
    assert(cudaSuccess == cudaStat);
    
    printf("after DDgels: info_gpu = %d\n",info_gpu);
    printMatrix(n,X,"X");
    
    assert(CUSOLVER_STATUS_SUCCESS == cusolver_status);
    
    if (d_A) cudaFree(d_A);
    if (d_B) cudaFree(d_B);
    if (d_X) cudaFree(d_X);
    if (d_info) cudaFree(d_info);
    if (d_work) cudaFree(d_work);
    if (cusolverH) cusolverDnDestroy(cusolverH);
    cudaDeviceReset();
    return 0;
}

解决方法

不幸的是，cuSolver设置存在不一致之处，从而造成了此问题。 有一种方法可以避免这种问题，方法是调用专家API“ cusolverDnIRSXgels”“ cusolverDnIRSXgels_bufferSize”，从而为用户提供更多控制权。

因此在您的代码中替换

    cusolver_status = cusolverDnDDgels_bufferSize(cusolverH,m,n,nrhs,d_A,lda,d_B,ldb,d_X,ldx,d_work,&lwork);
    assert(CUSOLVER_STATUS_SUCCESS == cusolver_status);

    // Allocate workspace
    cudaStat = cudaMalloc((void**)&d_work,lwork);
    assert(cudaSuccess == cudaStat);

    // Run solver
    cusolver_status = cusolverDnDDgels(cusolverH,lwork,&niter,d_info);
    printf("gels status: %d\n",int(cusolver_status));

作者

    // create the params and info structure for the expert interface
    cusolverDnIRSParams_t gels_irs_params;
    cusolverDnIRSParamsCreate( &gels_irs_params );
    cusolverDnIRSInfos_t gels_irs_infos;
    cusolverDnIRSInfosCreate( &gels_irs_infos );

    // Set the main and the low precision of the solver DSgels 
    // D is for double S for single precision thus 
    // main_precision is CUSOLVER_R_FP64,low_precision is CUSOLVER_R_FP32
    cusolverDnIRSParamsSetSolverPrecisions( gels_irs_params,CUSOLVER_R_64F,CUSOLVER_R_32F );
    // Set the refinement solver.
    cusolverDnIRSParamsSetRefinementSolver( gels_irs_params,CUSOLVER_IRS_REFINE_CLASSICAL );
    // Get work buffer size
    cusolver_status = cusolverDnIRSXgels_bufferSize(cusolverH,gels_irs_params,&lwork);
    assert(CUSOLVER_STATUS_SUCCESS == cusolver_status);
    // Allocate workspace
    cudaStat = cudaMalloc((void**)&d_work,lwork);
    assert(cudaSuccess == cudaStat);
    // Run solver
    cusolver_status = cusolverDnIRSXgels(cusolverH,gels_irs_infos,(void *)d_A,(void *)d_B,(void *)d_X,int(cusolver_status));

还请注意，当m> n是一个超额认购的方程组时，您就不能选择RHS然后找到SO，因此，最好选择SOL，生成RHS = A * SOL，然后使用RHS和与SOL比较。

还要注意，LDX应该> = max（m，n）

我通过以下方式修改了您的代码：

#include <stdio.h>
#include <stdlib.h>
#include <assert.h>

#include <cuda_runtime.h>
#include <cusolverDn.h>


#define USE_BUG
typedef double mt;

#ifndef max
#define max(a,b) ((a) > (b) ? (a) : (b))
#endif

void matvec(int m,int n,int nrhs,const mt* A,int lda,mt *X,int ldx,mt *B,int ldb)
{
    mt sum[nrhs];

    for (int row = 0; row < m; row++) {
        for (int r = 0; r < nrhs; r++) sum[r] = 0.0;
        for (int col = 0; col < n; col++) {
            for (int r = 0; r < nrhs; r++){
                sum[r] += A[row + col * lda] * X[col + r*ldx];
            }
        }
        for (int r = 0; r < nrhs; r++) B[row + r*ldb] = sum[r];
    }
}

mt check_solution(int n,mt *ref,int ldr,int ldx)
{
    mt error=0.0;
    for (int r = 0; r < nrhs; r++){
        for (int i = 0; i < n; i++) {
            error = max(error,abs(ref[i+r*ldr] - X[i+r*ldr]));
        }
    }
    return error;
}


void printMatrix(int m,const char* name)
{
    for (int row = 0; row < m; row++) {
        for (int col = 0; col < n; col++) {
            mt Areg = A[row + col * lda];
            printf("%s(%d,%d) = %f\n",name,row + 1,col + 1,Areg);
        }
    }
}





int main(int argc,char*argv[])
{
#ifndef USE_BUG
        // 3x3 example works fine
    const int m = 3;
    const int n = 3;
    mt A[m*n] = { 1.0,4.0,2.0,5.0,1.0,3.0,6.0,1.0 };
    mt sol[n] = { 6.0,15.0,4.0 };
#else
    // 4x3 example d_info/info_gpu returns -1
    const int m = 4;
    const int n = 3;
    mt A[m*n] = { 1.0,2.0 };
    mt sol[n] =   { 6.0,4.0 };
#endif
    mt X[n];
    mt B[m];

    int lda = m;
    int ldb = max(m,n);
    int ldx = max(m,n);
    int nrhs = 1;
    int niter = 0;
    int info_gpu = 0;
    size_t lwork = 0;

    mt *d_A = NULL;
    mt *d_B = NULL;
    mt *d_X = NULL;
    mt *d_work = NULL;
    int* d_info = NULL;

    // compute B = A*sol
    matvec(m,A,sol,B,ldb);

    cusolverDnHandle_t cusolverH = NULL;
    cudaError_t cudaStat = cudaSuccess;
    cusolverStatus_t cusolver_status = CUSOLVER_STATUS_SUCCESS;

    cusolver_status = cusolverDnCreate(&cusolverH);
    assert(CUSOLVER_STATUS_SUCCESS == cusolver_status);

    // Allocate space in the GPU
    cudaStat = cudaMalloc((void**)&d_A,sizeof(mt) * m * n);
    assert(cudaSuccess == cudaStat);

    cudaStat = cudaMalloc((void**)&d_B,sizeof(mt) * m * nrhs);
    assert(cudaSuccess == cudaStat);

    cudaStat = cudaMalloc((void**)&d_X,sizeof(mt) * n * nrhs);
    assert(cudaSuccess == cudaStat);

    cudaStat = cudaMalloc((void**)&d_info,sizeof(int));
    assert(cudaSuccess == cudaStat);

    // Copy matrices into GPU space
    cudaStat = cudaMemcpy(d_A,sizeof(mt) * m * n,cudaMemcpyHostToDevice);
    assert(cudaSuccess == cudaStat);
    cudaStat = cudaMemcpy(d_B,sizeof(mt) * m * nrhs,cudaMemcpyHostToDevice);
    assert(cudaSuccess == cudaStat);

    #if 1
    // =======================================================
    // create the params and info structure for the expert interface
    cusolverDnIRSParams_t gels_irs_params;
    cusolverDnIRSParamsCreate( &gels_irs_params );
    cusolverDnIRSInfos_t gels_irs_infos;
    cusolverDnIRSInfosCreate( &gels_irs_infos );

    // Set the main and the low precision of the solver DSgels 
    // D is for double S for single precision thus 
    // main_precision is CUSOLVER_R_FP64,int(cusolver_status));
    #else

    // Get work buffer size
    cusolver_status = cusolverDnDDgels_bufferSize(cusolverH,int(cusolver_status));
    #endif

    // Sync threads
    cudaStat = cudaDeviceSynchronize();
    assert(cudaSuccess == cudaStat);

    // Copy GPU info
    cudaStat = cudaMemcpy(&info_gpu,d_info,sizeof(int),cudaMemcpyDeviceToHost);
    assert(cudaSuccess == cudaStat);

    // Get solved data
    cudaStat = cudaMemcpy(X,sizeof(mt) * n * nrhs,cudaMemcpyDeviceToHost);
    assert(cudaSuccess == cudaStat);

    printf("after gels: info_gpu = %d\n",info_gpu);
    printf("after gels: niter    = %d\n",niter);
    printf("after gels: error    = %e\n",check_solution(n,X,ldx));
    printMatrix(3,"X");


    if (d_A) cudaFree(d_A);
    if (d_B) cudaFree(d_B);
    if (d_X) cudaFree(d_X);
    if (d_info) cudaFree(d_info);
    if (d_work) cudaFree(d_work);
    if (cusolverH) cusolverDnDestroy(cusolverH);
    cudaDeviceReset();
    return 0;
}

使用nvcc -o test test.cu -lcusolver编译