CUBLAS: Incorrect inversion for matrix with zero pivot

by

There seems to be a bug in the current CUBLAS library implementation of cublas<t>getrfBatched for matrices of dimension (n) such that 3<=n<=16, when there is a “zero pivot” as you say.

A possible workaround is to “identity-extend” your A matrix to be inverted, when n<17, to a size of 17×17 (using matlab nomenclature):

 LU = getrf( [A 0 ; 0 I]);

continuing, you can then use cublas<t>getriBatched in an “ordinary” fashion:

 invA = getri( LU(1:3,1:3) )

(You can also leave everything at n=17, call getri that way, and then extract the result as the first 3×3 rows and columns of invA.)

Here is a fully worked example, borrowing from the code you supplied, showing the inversion of your supplied 3×3 zero_pivot matrix, using the zero_pivot_war matrix as an “identity-extended” workaround:

$ cat t340.cu
#include <cstdio>
#include <cstdlib>
#include <cuda_runtime.h>
#include <cublas_v2.h>

#define cudacall(call)                                                                                                          \
    do                                                                                                                          \
    {                                                                                                                           \
        cudaError_t err = (call);                                                                                               \
        if(cudaSuccess != err)                                                                                                  \
        {                                                                                                                       \
            fprintf(stderr,"CUDA Error:\nFile = %s\nLine = %d\nReason = %s\n", __FILE__, __LINE__, cudaGetErrorString(err));    \
            cudaDeviceReset();                                                                                                  \
            exit(EXIT_FAILURE);                                                                                                 \
        }                                                                                                                       \
    }                                                                                                                           \
    while (0)

#define cublascall(call)                                                                                        \
    do                                                                                                          \
    {                                                                                                           \
        cublasStatus_t status = (call);                                                                         \
        if(CUBLAS_STATUS_SUCCESS != status)                                                                     \
        {                                                                                                       \
            fprintf(stderr,"CUBLAS Error:\nFile = %s\nLine = %d\nCode = %d\n", __FILE__, __LINE__, status);     \
            cudaDeviceReset();                                                                                  \
            exit(EXIT_FAILURE);                                                                                 \
        }                                                                                                       \
                                                                                                                \
    }                                                                                                           \
    while(0)


void invert_device(float* src_d, float* dst_d, int n)
{
    cublasHandle_t handle;
    cublascall(cublasCreate_v2(&handle));

    int batchSize = 1;

    int *P, *INFO;

    cudacall(cudaMalloc<int>(&P,17 * batchSize * sizeof(int)));
    cudacall(cudaMalloc<int>(&INFO,batchSize * sizeof(int)));

    int lda = 17;

    float *A[] = { src_d };
    float** A_d;
    cudacall(cudaMalloc<float*>(&A_d,sizeof(A)));
    cudacall(cudaMemcpy(A_d,A,sizeof(A),cudaMemcpyHostToDevice));

    cublascall(cublasSgetrfBatched(handle,17,A_d,lda,P,INFO,batchSize));

    int INFOh = 0;
    cudacall(cudaMemcpy(&INFOh,INFO,sizeof(int),cudaMemcpyDeviceToHost));

    if(INFOh == 17)
    {
        fprintf(stderr, "Factorization Failed: Matrix is singular\n");
        cudaDeviceReset();
        exit(EXIT_FAILURE);
    }

    float* C[] = { dst_d };
    float** C_d;
    cudacall(cudaMalloc<float*>(&C_d,sizeof(C)));
    cudacall(cudaMemcpy(C_d,C,sizeof(C),cudaMemcpyHostToDevice));

    cublascall(cublasSgetriBatched(handle,n,A_d,lda,P,C_d,n,INFO,batchSize));

    cudacall(cudaMemcpy(&INFOh,INFO,sizeof(int),cudaMemcpyDeviceToHost));

    if(INFOh != 0)
    {
        fprintf(stderr, "Inversion Failed: Matrix is singular\n");
        cudaDeviceReset();
        exit(EXIT_FAILURE);
    }

    cudaFree(P), cudaFree(INFO), cublasDestroy_v2(handle);
}

void invert(float* src, float* dst, int n)
{
    float* src_d, *dst_d;

    cudacall(cudaMalloc<float>(&src_d,17 * 17 * sizeof(float)));
    cudacall(cudaMemcpy(src_d,src,17 * 17 * sizeof(float),cudaMemcpyHostToDevice));
    cudacall(cudaMalloc<float>(&dst_d,n * n * sizeof(float)));

    invert_device(src_d,dst_d,n);

    cudacall(cudaMemcpy(dst,dst_d,n * n * sizeof(float),cudaMemcpyDeviceToHost));

    cudaFree(src_d), cudaFree(dst_d);
}

void test_invert()
{
    const int n = 3;

    //Random matrix with full pivots
/*    float full_pivots[n*n] = { 0.5, 3, 4,
                                1, 3, 10,
                                4 , 9, 16 };

    //Almost same as above matrix with first pivot zero
    float zero_pivot[n*n] = { 0, 3, 4,
                              1, 3, 10,
                              4 , 9, 16 };

    float zero_pivot_col_major[n*n] = { 0, 1, 4,
                                        3, 3, 9,
                                        4 , 10, 16 };
*/
    float zero_pivot_war[17*17] = {
                                        0,3,4,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
                                        1,3,10,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
                                        4,9,16,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
                                        0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,
                                        0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,
                                        0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,
                                        0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,
                                        0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,
                                        0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,
                                        0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,
                                        0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
                                        0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,
                                        0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,
                                        0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,
                                        0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,
                                        0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,
                                        0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1 };
/*
    float another_zero_pivot[n*n] = { 0, 3, 4,
                                      1, 5, 6,
                                      9, 8, 2 };

    float another_full_pivot[n * n] = { 22, 3, 4,
                                        1, 5, 6,
                                        9, 8, 2 };

    float singular[n*n] = {1,2,3,
                           4,5,6,
                           7,8,9};
*/
    float result[n*n];

    //Select matrix by setting "a"
    float* a = zero_pivot_war;

    fprintf(stdout, "Input:\n\n");
    for(int i=0; i<n; i++)
    {
        for(int j=0; j<n; j++)
            fprintf(stdout,"%f\t",a[i*17+j]);
        fprintf(stdout,"\n");
    }

    fprintf(stdout,"\n\n");

    invert(a,result,n);

    fprintf(stdout, "Inverse:\n\n");
    for(int i=0; i<n; i++)
    {
        for(int j=0; j<n; j++)
            fprintf(stdout,"%f\t",result[i*n+j]);
        fprintf(stdout,"\n");
    }

}

int main()
{
    test_invert();

//    int n;  scanf("%d",&n);
    return 0;
}
$ nvcc -arch=sm_20 -o t340 t340.cu -lcublas
$ cuda-memcheck ./t340
========= CUDA-MEMCHECK
Input:

0.000000        3.000000        4.000000
1.000000        3.000000        10.000000
4.000000        9.000000        16.000000


Inverse:

-0.700000       -0.200000       0.300000
0.400000        -0.266667       0.066667
-0.050000       0.200000        -0.050000
========= ERROR SUMMARY: 0 errors
$

The above result appears to me to be correct based on a simple test elsewhere.

I don’t have any further technical details about the nature of the possible bug in CUBLAS. From what I can tell, it is present in both CUDA 5.5 and CUDA 6.0 RC. Detailed bug discussions for NVIDIA-supplied assets (e.g. CUBLAS library) should be taken up on the NVIDIA developer forums or directly at the bug filing portal on developer.nvidia.com (you must be a registered developer to file a bug).

More Related Contents:

  1. What is the canonical way to check for errors using the CUDA runtime API?
  2. CUDA compute capability requirements
  3. Cuda atomics change flag
  4. GPU Emulator for CUDA programming without the hardware [closed]
  5. How is CUDA memory managed?
  6. Modifying registry to increase GPU timeout, windows 7
  7. How can I implement a custom atomic function involving several variables?
  8. Using GPU from a docker container?
  9. How to measure the inner kernel time in NVIDIA CUDA?
  10. Passing Host Function as a function pointer in __global__ OR __device__ function in CUDA
  11. How do CUDA blocks/warps/threads map onto CUDA cores?
  12. How to verify CuDNN installation?
  13. What is the purpose of using multiple “arch” flags in Nvidia’s NVCC compiler?
  14. How to get the CUDA version?
  15. How do I use Nvidia Multi-process Service (MPS) to run multiple non-MPI CUDA applications?
  16. Copy an object to device?
  17. What kind of variables consume registers in CUDA?
  18. CUDA: How to use -arch and -code and SM vs COMPUTE
  19. In CUDA, what is memory coalescing, and how is it achieved?
  20. How to use 2D Arrays in CUDA?
  21. How can I check the progress of matrix multiplication?
  22. cudaMemset() – does it set bytes or integers?
  23. Using maximum shared memory in Cuda
  24. Cuda Mutex, why deadlock?
  25. Inconsistency of IDs between ‘nvidia-smi -L’ and cuDeviceGetName()
  26. CUDA allocation alignment is 256 bytes – seriously?
  27. printf() in my CUDA kernel doesn’t result produce any output
  28. Multiply Rectangular Matrices in CUDA
  29. Which Compute Capability is supported by which CUDA versions?
  30. “invalid configuration argument ” error for the call of CUDA kernel?

Leave a Comment