Cuda In-Situ memory race issue for algorithms such as convolution of morphologicam dilation

Question

I wrote a dilation kernel in CUDA and it works well when my input and my output images are different buffers, but I am facing what I understand to be a memory race issue when I call my kernel in an in-situ case, i.e. the input and the output buffers point to the same memory location.

I tried :

a. using cooperative groups,

b. using a mutex and an atomic addition but as suggested in this paper and in several sources on the web,

c. using a lock-free inter-block synchronization, the synchronization proposed in this same paper.

All my attempts failed because :

a. did not work because my input buffer is a const pointer and I have a compilation error when I have to cast it into a void* parameter (which makes sense), so I could not go further.

b. did not work because I faced a wierd behaviour : I have 16x16 blocks, each with 32x32 threads. Synchronizing the blocks should increase the mutex to 256 but the program blocks after 48 atomic additions.

c. did not work because it seams to be no inter-block synchronization, although the code I used directly from the paper seems good to me. I could improve a little the race effect by adding some __syncthreads()

This is the dilation function ;

template <typename T>
__global__ void GenericDilate2dImg_knl(const ImageSizeInfo imgSizeInfo,
                                       volatile int* syncArrayIn, volatile int* syncArrayOut, 
                                       const unsigned long localSizeX, const unsigned long localSizeY,
                                       const int borderPolicyType, const T outOfImageValue,
                                       const struct StructuringElementInfo seInfo,
                                       const T* pInBuf, T* pOutBuf)
{
    // Extract sizeX, sizeY, etc. from imgSizeInfo
    SPLIT_SIZES_FROM_STRUCT(imgSizeInfo)

    // Declare the shared buffer pSharedBuf
    extern __shared__ char pSharedMem[];
    T* pSharedBuf = reinterpret_cast<T*>(pSharedMem);

    const unsigned long x = blockDim.x * blockIdx.x + threadIdx.x;
    const unsigned long y = blockDim.y * blockIdx.y + threadIdx.y;
    const unsigned long planIdx = blockDim.z * blockIdx.z + threadIdx.z;

    const unsigned long nbPlans = sizeZ * sizeC * sizeT;

    const unsigned long idx = x + y * sizeX + planIdx * sizeX*sizeY;

    // Copy the input image data into shared memory
    if (x < blockDim.x * gridDim.x && y < blockDim.y * gridDim.y && planIdx < blockDim.z * gridDim.z) {
        copyDataToSharedMemory2d(pInBuf, sizeX, sizeY, planIdx,
                                 localSizeX, localSizeY, 
                                 seInfo._paddingX, seInfo._paddingY,
                                 borderPolicyType, outOfImageValue,
                                 pSharedBuf);
    }

    // Wait to ensure that the copy is terminated
    if (pInBuf == pOutBuf) {
        // Grid synchronization for in-situ case
        //__gpu_sync(gridDim.x * gridDim.y);        // Use a mutex
        __gpu_sync2(1, syncArrayIn, syncArrayOut);  // Use a lock-free barrier
    }
    else
        // The input and ouput buffers point to different data 
        // -> we simply need to synchronize the threads inside the block
        __syncthreads();

    // Compute the convolution for pixels inside the image
    if (x < sizeX && y < sizeY && planIdx < nbPlans) {
        T vMax = 0;
        for (unsigned int curCoefIdx = 0; curCoefIdx < seInfo._nbOffsets; ++curCoefIdx) {
            const unsigned int sx = threadIdx.x + seInfo._paddingX + seInfo._pOffsetsX[curCoefIdx];
            const unsigned int sy = threadIdx.y + seInfo._paddingY + seInfo._pOffsetsY[curCoefIdx];
            const unsigned long sidx = sx + sy * localSizeX;
            const T curVal = pSharedBuf[sidx];
            vMax = (vMax > curVal ? vMax : curVal);
        }

        // Round the result
        pOutBuf[idx] = vMax;
    }
}

My function to copy from global to shared memory is :

template <typename T>
__device__ void copyDataToSharedMemory2d(const T* pInBuf,
                                         const unsigned long sizeX, const unsigned long sizeY, const unsigned long planIdx,
                                         const unsigned long localSizeX, const unsigned long localSizeY,
                                         const int paddingX, const int paddingY,
                                         const int borderPolicyType, const T outOfImageValue,
                                         T* pSharedBuf)
{
    const int x = blockDim.x * blockIdx.x + threadIdx.x;
    const int y = blockDim.y * blockIdx.y + threadIdx.y;
    const int localX = threadIdx.x;
    const int localY = threadIdx.y;

    // Fill the shared buffer tile by tile
    // A tile is related to the group size
    const unsigned int groupSizeX = blockDim.x;
    const unsigned int groupSizeY = blockDim.y;

    // For each tile
    for (int offsetY = 0; offsetY < localSizeY; offsetY += groupSizeY) {
        int curLocalY = localY + offsetY;
        int curGlobalY = y + offsetY - paddingY;
        for (int offsetX = 0; offsetX < localSizeX; offsetX += groupSizeX) {
            int curLocalX = localX + offsetX;
            int curGlobalX = x + offsetX - paddingX;

            // If the current coordinate is inside the shared sub-image
            if (curLocalX < localSizeX && curLocalY < localSizeY) {
                const int idx = curLocalX + curLocalY * localSizeX;
                pSharedBuf[idx] = getPixel2d(pInBuf, sizeX, sizeY, curGlobalX, curGlobalY, planIdx, borderPolicyType, outOfImageValue);
            }
        }
    }
}

Where getPixel2d allows me to manage the data out of the image:

template <typename T>
__device__
T getPixel2d(const T* pInBuf,
             const unsigned long sizeX, const unsigned long sizeY,
             const int x, const int y, const int z,
             const int borderPolicyType, const T outOfImageValue)
{
    int x_inside = x;
    if (x < 0 || x >= sizeX) {
        switch (borderPolicyType) {
        case 0://outside the image, there is a constant value
            return outOfImageValue;
        case 1://outside the image, we propagate the data at the image borders
            if (x < 0)
                x_inside = 0;
            else // x >= sizeX
                x_inside = sizeX - 1;
            break;
        case 2://Miror effect
            if (x < 0)
                x_inside = -(x + 1);
            else // x >= sizeX
                x_inside = sizeX - ((x - sizeX) + 1);
            break;
        }
    }

    // y-coordinate inside the image
    int y_inside = y;
    if (y < 0 || y >= sizeY) {
        switch (borderPolicyType) {
        case 0://outside the image, there is a constant value
            return outOfImageValue;
        case 1://outside the image, we propagate the data at the image borders
            if (y < 0)
                y_inside = 0;
            else // y >= sizeY
                y_inside = sizeY - 1;
            break;
        case 2://Miror effect
            if (y < 0)
                y_inside = -(y + 1);
            else // y >= sizeY
                y_inside = sizeY - ((y - sizeY) + 1);
            break;
        default: break;
        }
    }

    return pInBuf[x_inside + y_inside * sizeX + z * sizeX * sizeY];
}

and now, here are my inter-block synchronization functions :

// Using a mutex
__device__ volatile int g_mutex;
__device__ void __gpu_sync(int goalVal) {
    //thread ID in a block
    int tid_in_block = threadIdx.x * blockDim.y + threadIdx.y;
    // only thread 0 is used for synchronization
    if (tid_in_block == 0) {
        atomicAdd((int*)&g_mutex, 1);
        printf("[%d] %d Vs %d\n", blockIdx.x * gridDim.y + blockIdx.y, g_mutex, goalVal);
        //only when all blocks add 1 to g_mutex
        //will g_mutex equal to goalVal
        while (g_mutex </*!=*/ goalVal) {
            ;//Do nothing here
        }
    }
    __syncthreads();
}

// Lock-free barrier
__device__ void __gpu_sync2(int goalVal, volatile int* Arrayin, volatile int* Arrayout) {
    // thread ID in a block
    int tid_in_blk = threadIdx.x * blockDim.y + threadIdx.y;
    int nBlockNum = gridDim.x * gridDim.y;
    int bid = blockIdx.x * gridDim.y + blockIdx.y;
    // only thread 0 is used for synchronization
    if (tid_in_blk == 0) {
        Arrayin[bid] = goalVal;
    }
    if (bid == 1) {
        if (tid_in_blk < nBlockNum) {
            while (Arrayin[tid_in_blk] != goalVal) {
                ;//Do nothing here
            }
        }
        __syncthreads();
        if (tid_in_blk < nBlockNum) {
            Arrayout[tid_in_blk] = goalVal;
        }
    }
    if (tid_in_blk == 0) {
        while (Arrayout[bid] != goalVal) {
            ;//Do nothing here
        }
    }
    __syncthreads();
}

The image I get for in-situ calculation is :

I used a 11x15 structuring emelent and the size of the shared buffer is (nbThreadsPerBlock+2*paddindX) * (nbThreadsPerBlock+2*paddindY). The wrong result (showed by the arrows) appears at the top of some blocks, but always at the same location and with the same values. I'd expect a more random result for memory race effect...

Is there a better approach to manage in-situ calculation or any reason that would prevent the grid synchronization to work?

EDIT The size of the image I used is 510x509 and I run my code on a NVidia Quadro RTX 5000.

Answer 1

I would normally suggest minimal reproducible example for a question like this, as well as an indication of the GPU you are running on, but we can probably proceed without that. In short, what you are trying to do will not work reliably, as you've already discovered.

You have chosen a thread strategy of assigning one thread in your grid per output point:

    pOutBuf[idx] = vMax;

which is sensible and fine. I imagine based on this:

I have 16x16 blocks, each with 32x32 threads.

that your input images are 512x512 (16x32 threads in each direction, one thread per output point).

And as you've already stated, you have 256 blocks (each of 1024 threads) in your grid. Furthermore, for the in-situ case, we can simplify your kernel to the following pseudo-code:

__global__ void GenericDilate2dImg_knl(...){
  read_in_image();
  grid_wide_sync();
  write_out_image();
}

For such a methodology to work, then, the read_in_image() step must be able to read the entire image, before any writing occurs. However your methodology will not work in the general case, and evidently not on your specific GPU, either. In order to read in the entire image as per above, we must have every threadblock in the grid simultaneously resident on the SMs in your GPU. All 256 blocks need to be deposited, and running on an SM. But the GPU provides no inherent guarantees of such a thing. If your GPU has, for example 24 SMs in it, each of which can hold a maximum of 2048 threads, then your GPU would have a "running" or "instantaneous" capacity of 24*2048 threads, or 48 of your threadblocks. There would not be enough room for all 256 threadblocks to be running. Not only does your algorithm depend on that, but all 3 of your grid sync methods depend on that notion as well.

The fact that your 2nd grid sync method stops after 48 "atomic additions" suggested the example numbers above to me. It's a plausible proximal explanation for why that method may have failed that way: your GPU only allowed 48 of your threadblocks to be resident, and the other 208 threadblocks were waiting in the wings, not yet deposited on any SM, and therefore not allowing any of their threads to run. Those threads in those 208 threadblocks need to run to pick up the relevant input data, as well as to satisfy the requirements of the grid-wide sync. But they are not running, because they are waiting for room to open up on a SM. And room never opens up on a SM, because the full SMs have threadblocks that are waiting at the grid sync point. So you have deadlock.

This problem is not easily solvable in the general case. Any grid sync mechanism, including cooperative groups, has an inherent requirement that all threadblocks be actually simultaneously schedulable on your particular GPU. Therefore in the general case, where we don't know the data set size or the GPU we will be running on, the problem is quite difficult.

One possible approach is to divide your input data set into regions, and have your kernel process a region at a time. This may require multiple grid syncs, one to handle the in/out division in each region, and one to handle the progression of the kernel as it steps through regions. You would also have to handle the region edges carefully.

Another possible approach if you know the specifics of the data set size and the GPU you are running on, is just to make sure you are running on a GPU "large enough" to handle the data set size. For example, an A100 GPU could probably have as many 216 blocks simultaneously resident, so for that case you could handle a somewhat smaller image size, perhaps 14x32=448 height and 448 width dimensions.

Given that these approaches for in-place or in-situ work for this particular example require considerable complexity, I personally would be strongly motivated to use the methodology where output is different than input. That approach will likely run noticeably quicker as well. A grid wide sync is not a "free" construct from a performance perspective.