Simulations on the GPU sometimes return all zero images

[JanWP]

Simulations performed on the GPU sometimes return images without any signal, where the same simulation performed on the CPU returns reasonable result.

Software environment:
MRiLab v1.3, Jan 10th 2017
Matlab R2016a
Windows 7 Professional SP1
CUDA 7.5
IPP 2017
Microsoft Visual Studio 2012

Hardware environment:
i7 3930K, 32GB ram
GTX Titan Black (GK110, 15 SMXs, 2304 shaders, 192 TMUs, 48 ROPs)
GTX Titan (GK110, 14 SMXs, 2880 shaders, 240 TMUs, 48 ROPs)

To reproduce:
Start MRiLab in Matlab. In my case the GTX Titan Black card is used by default.

1. Load > Load Phantom Example > Brain (High Resolution 215x180x180)
2. Click "<Localizer", then "Update"
3. "Sequence", GradientEcho > PSD_GRE3D, "Accept", "Update", "Scan"
   -> The output shown in the preview contains only noise.
4. "Hardware" tab > NoiseLevel -> set to 0, "Update", "Scan"
   -> The output shown in the preview is an image of all zeros

The same steps performed with the GTX Titan yield the same results.
The same steps performed with the CPU solver yield a reasonable image.
The same steps performed with the standard resolution brain phantom calculated on either one of the GPUs also yield a reasonable image that is visually indistinguishable from the one calculated on the CPU.

Sometimes, even high resolution images can be calculated on the GPU, but I haven't yet figured out the exact circumstances. Very similar simulations sometimes work, and sometimes they don't.

I'm willing to spend time helping to debug this, but I would need some suggestion of where to look or what tests to perform. I consider myself an expert in Matlab and have some limited experience writing C++/CUDA code to create mex files.

Jan

[leoliuf]

Jan, Thank for willing to debug. I would say that the returning zero will likely happen as the CUDA kernel wasn't launching properly. Did you recompile the mex code for your card? Or you have the out of box mex file working? You have better GPU card than mine.

[leoliuf]

If you are debugging using VS, you might consider using nVidia Nsight tool. It works for Mex code after you attach matlab process.

[JanWP]

Thanks for your quick feedback. I cannot tell from the execution time that there is a difference between situations where the simulation works, and where it doesn't work. I did compile everything on this computer. So the cuda kernel is definitely called and doing something in every case, as far as I can tell.

Unfortunately, I have no experience at all with using nVidia Nsight, and while that seems like a potentially useful thing to learn, I prefer not investing time in learning using new tools right now. Besides, my evaluation version of VS expired...

I guess I'll start by looking for call of the kernel in the Matlab code and see what input it receives and returns.

[JanWP]

Tying to debug the hard way. I find that the success or failure of the simulation depends on the x-y size of the object domain, and seems independent of sequence parameters. Basically, the simulation fails, as soon as dimBlockImg.y exceeds 896. For information, on my machine:

deviceProp.multiProcessorCount = 15;
deviceProp.maxThreadsPerBlock = 1024;
deviceProp.regsPerBlock = 65536;

Could it be that there is a limitation other than threads per block or registers per block? Shared memory size maybe (48 KB per multiprocessor in my case)? My understanding is that it does not seem to be an issue of register count, since if the kernel used floor(65536/896)=73 registers instead of the 63 assumed in the calculations of the GPU kernel, dimBlockImg.y=897 should work as well but doesn't.

[JanWP]

Maybe there is an issue with register count, after all. How did you choose the value of 63 registers coded in DoScanAtGPU.cu? My understanding is that this value depends on the optimizations and code reordering taking place in the CUDA compiler and thus on the target architecture.

In DoScanAtGPU.log in the build tree I find the following:

    ptxas info    : 0 bytes gmem
    ptxas info    : Compiling entry function '_Z20BlochKernelNormalGPUfPdPfS0_S0_S0_S0_S0_S0_S0_S0_S0_S0_S0_S0_S0_S0_S0_ffS0_S0_fiiiiiiiiii' for 'sm_20'
    ptxas info    : Function properties for _Z20BlochKernelNormalGPUfPdPfS0_S0_S0_S0_S0_S0_S0_S0_S0_S0_S0_S0_S0_S0_S0_ffS0_S0_fiiiiiiiiii
             0 bytes stack frame, 0 bytes spill stores, 0 bytes spill loads
    ptxas info    : Used 61 registers, 244 bytes cmem[0], 16 bytes cmem[16]
    ...
    ptxas info    : 0 bytes gmem
    ptxas info    : Compiling entry function '_Z20BlochKernelNormalGPUfPdPfS0_S0_S0_S0_S0_S0_S0_S0_S0_S0_S0_S0_S0_S0_S0_ffS0_S0_fiiiiiiiiii' for 'sm_30'
    ptxas info    : Function properties for _Z20BlochKernelNormalGPUfPdPfS0_S0_S0_S0_S0_S0_S0_S0_S0_S0_S0_S0_S0_S0_S0_ffS0_S0_fiiiiiiiiii
        0 bytes stack frame, 0 bytes spill stores, 0 bytes spill loads
    ptxas info    : Used 62 registers, 532 bytes cmem[0], 32 bytes cmem[2]
    ...
    ptxas info    : 0 bytes gmem
    ptxas info    : Compiling entry function '_Z20BlochKernelNormalGPUfPdPfS0_S0_S0_S0_S0_S0_S0_S0_S0_S0_S0_S0_S0_S0_S0_ffS0_S0_fiiiiiiiiii' for 'sm_35'
    ptxas info    : Function properties for _Z20BlochKernelNormalGPUfPdPfS0_S0_S0_S0_S0_S0_S0_S0_S0_S0_S0_S0_S0_S0_S0_ffS0_S0_fiiiiiiiiii
        0 bytes stack frame, 0 bytes spill stores, 0 bytes spill loads
    ptxas info    : Used 72 registers, 532 bytes cmem[0], 32 bytes cmem[2]
    ...
    ptxas info    : 0 bytes gmem
    ptxas info    : Compiling entry function '_Z20BlochKernelNormalGPUfPdPfS0_S0_S0_S0_S0_S0_S0_S0_S0_S0_S0_S0_S0_S0_S0_ffS0_S0_fiiiiiiiiii' for 'sm_50'
    ptxas info    : Function properties for _Z20BlochKernelNormalGPUfPdPfS0_S0_S0_S0_S0_S0_S0_S0_S0_S0_S0_S0_S0_S0_S0_ffS0_S0_fiiiiiiiiii
        0 bytes stack frame, 0 bytes spill stores, 0 bytes spill loads
    ptxas info    : Used 61 registers, 532 bytes cmem[0], 44 bytes cmem[2]

My cards have a compute capability of 3.5, which means that actually 72 registers will be used. So I guess I have a choice of either changing the hardcoded value in DoScanAtGPU.cu, or to pass a --maxrregcount limit at compile time. Where would be the best place to put such a compile option?

[JanWP]

Solved

I changed Lib\src\engine\GPUEngine\CMakeLists.txt to include --maxrregcount=63:

IF (IPP_FOUND) # use ipp
    # compile kernels for shader models 20 30 35 50, support double-precision floating-point operation
    set(CUDA_NVCC_FLAGS -Xcompiler -fPIC -use_fast_math -gencode=arch=compute_20,code="sm_20,compute_20"  
                          -gencode=arch=compute_30,code="sm_30,compute_30" 
                          -gencode=arch=compute_35,code="sm_35,compute_35" 
                          -gencode=arch=compute_50,code="sm_50,compute_50" 
                          --ptxas-options=-v --maxrregcount=63 -DMATLAB_MEX_FILE -DIPP)
else (IPP_FOUND) # use framewave
    # compile kernels for shader models 20 21 30 35 50, support double-precision floating-point operation
    set(CUDA_NVCC_FLAGS -Xcompiler -fPIC -use_fast_math -gencode=arch=compute_20,code="sm_20,compute_20"
                          -gencode=arch=compute_30,code="sm_30,compute_30" 
                          -gencode=arch=compute_35,code="sm_35,compute_35" 
                          -gencode=arch=compute_50,code="sm_50,compute_50" 
                          --ptxas-options=-v --maxrregcount=63 -DMATLAB_MEX_FILE -DFW)
endif (IPP_FOUND)

Now the code produces correct results up to the maximum of dimBlockImg.y=1024.

BTW, my compiler complains that -fPIC is an unknown option.