Testing CUDA 11 cusolverDnDSgels()

Question

Trying to make sense of cusolverDnDSgels function. If I run it with simple 3x3 example as in the docs it works, but when I run it with my data then d_info returns -1 which as the docs says if d_info = -i then i-th argument is not valid.

Bellow I posted the code with 3 by 3 and 4 by 3 matrices where the former works and second doesn't.

As a reference I used this web site calculator 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, 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, 4.0, 2.0, 2.0, 5.0, 1.0, 3.0, 6.0, 1.0, 5.0, 1.0, 2.0 };
    //double B[4] = { 6.0, 15.0, 4.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, m, n, nrhs, d_A, lda, d_B, ldb, d_X, ldx, d_work, 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, d_X, sizeof(double) * n * nrhs, cudaMemcpyDeviceToHost);
    assert(cudaSuccess == cudaStat);
    
    printf("after DDgels: info_gpu = %d\n", info_gpu);
    printMatrix(n, nrhs, X, ldx, "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;
}

Answer 1

Unfortunately, there is an inconsistency in cuSolver setting creating this issue. There is a way to avoid such issue by calling the expert API "cusolverDnIRSXgels" "cusolverDnIRSXgels_bufferSize" that give the user more control.

Thus in your code replace

    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, m, n, nrhs, d_A, lda, d_B, ldb, d_X, ldx, d_work, lwork, &niter, d_info);
    printf("gels status: %d\n", int(cusolver_status));

by

    // 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, m, n, nrhs, &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_params, gels_irs_infos, m, n, nrhs, (void *)d_A, lda, (void *)d_B, ldb, (void *)d_X, ldx, d_work, lwork, &niter, d_info);
    printf("gels status: %d\n", int(cusolver_status));

Also note that when m>n it is an oversubscribed system of equation so you cannot choose RHS and then find SO, it is better for your test to choose a SOL, generate RHS=A*SOL then solve using the RHS and compare to SOL.

Also note that LDX should be >= max(m,n)

I modified your code by the following:

#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, int nrhs, mt *ref, int ldr, mt *X, 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, int n, const mt* A, int lda, 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, 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, 4.0, 2.0, 2.0, 5.0, 1.0, 3.0, 6.0, 1.0, 5.0, 1.0, 2.0 };
    mt sol[n] =   { 6.0, 15.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, n, nrhs, A, lda, sol, ldx, 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, A, sizeof(mt) * m * n, cudaMemcpyHostToDevice);
    assert(cudaSuccess == cudaStat);
    cudaStat = cudaMemcpy(d_B, 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, 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, m, n, nrhs, &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_params, gels_irs_infos, m, n, nrhs, (void *)d_A, lda, (void *)d_B, ldb, (void *)d_X, ldx, d_work, lwork, &niter, d_info);
    printf("gels status: %d\n", int(cusolver_status));
    #else

    // Get work buffer size
    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, m, n, nrhs, d_A, lda, d_B, ldb, d_X, ldx, d_work, lwork, &niter, d_info);
    printf("gels status: %d\n", 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, d_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, nrhs, sol, ldx, X, ldx));
    printMatrix(3, nrhs, X, ldx, "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;
}

compile using nvcc -o test test.cu -lcusolver