Error using pagefun with mldivide

Hello all,

I am trying to solve a system of linear equations A*x=B on the GPU. Example code:

l_z = 6401;
B = zeros(80, 1, l_z, 'gpuArray');
B(80,1) = 1;

A = randn(80,80,l_z, 'gpuArray') + 1i* randn(80,80,l_z, 'gpuArray');

x = pagefun(@mldivide, A, B);

Which causes Matlab to throw the following error:

Error using gpuArray/pagefun
Call to zgetrsBatched in CUBLAS failed with error status: CUBLAS_STATUS_INTERNAL_ERROR.

Error in problem_GPU (line 7)
x = pagefun(@mldivide, A, B);

Since this error is somewhat cryptic, I tried several (random) things and was able to narrow it down a bit. The error occcurs only if l_z is greater than 682. That threshold stays the same even when I change the dimensions of A and B from 80 to 100, or when I make B a square matrix. Also, it makes no difference if A is complex or not.

I have absolutely no idea what the problem could be here and would be grateful for some help.

I am running Matlab R2021a with the latest updates and also updated my graphics driver.

>> gpuDevice()

ans = 

  CUDADevice with properties:

                      Name: 'Quadro P2200'
                     Index: 1
         ComputeCapability: '6.1'
            SupportsDouble: 1
             DriverVersion: 11.6000
            ToolkitVersion: 11
        MaxThreadsPerBlock: 1024
          MaxShmemPerBlock: 49152
        MaxThreadBlockSize: [1024 1024 64]
               MaxGridSize: [2.1475e+09 65535 65535]
                 SIMDWidth: 32
               TotalMemory: 5.3685e+09
           AvailableMemory: 4.3240e+09
       MultiprocessorCount: 10
              ClockRateKHz: 1493000
               ComputeMode: 'Default'
      GPUOverlapsTransfers: 1
    KernelExecutionTimeout: 1
          CanMapHostMemory: 1
           DeviceSupported: 1
           DeviceAvailable: 1
            DeviceSelected: 1

[Damian]

Of course there is the possibility of a workaround:

temp = pagefun(@inv, A);
x = pagefun(@mtimes, temp, B);

However, this doubles the computing time and, as I assume, is discouraged for a reason in the help section for inv():

x = A\b is computed differently than x = inv(A)*b and is recommended for solving systems of linear equations.

Answer 1

[Joss Knight]

This is a known bug in cublas introduced in CUDA 11 and fixed in the next MATLAB release. The workaround, as you correctly discovered, is to break up your computation into multiple calls with no more than 682 pages. Alternatively if reducing the size of your matrices is an option, they need to be reduced to size 64x64 or below.

Since this is being detected by customers we will publish a bug report to make the bug and its workaround more discoverable. Apologies for the inconvenience.

[Damian]

This clears it up. Thanks for the explanation!