Resolve 1D shared memory bank conflict with paddling

Accelerated Computing CUDA CUDA Programming and Performance
,
1

Hi Experts,

I am studying how to resolve bank conflict by paddling.
I want to store 72 cuDoubleComplex values into shared memory, and this is my memory layout.

image
image1845×979 90.4 KB

The columns are bank number from 0 to 31, I know that double complex is 16 bytes.

As a beginner on this topic, the memory arrangement I can currently think of is like this but my question is, I have 72 elements to store into shared memory but this way can only save 32 elements.

These data need to be reused, and 72 elements should fit into shared memory.

Currently, my only solution for bank conflicts is paddling, which, although it wastes additional space, seems to be the most straightforward solution.

Are there any other methods to do this?

2

There is no way to assess shared memory bank conflicts based strictly on storage ordering. An access pattern must also be supplied. For a specific load from or store to shared memory, what address/location in shared memory will each thread in a warp access?

As long as you have adjacent/contiguous locations/addresses across the threads in a warp, you need not worry about bank conflicts, regardless of the size of the data type (at least up to 16 bytes per thread, naturally aligned). That is not the only bank-conflict-free pattern that is possible.

Unit 4 of this online training series covers shared bank conflicts.

3

Hi @Robert_Crovella,

I watched this video before, let me explain my scenario more detailed.

I arranged several blocks, each of them has 64 threads.
I will first store these 72 elements in shared memory to avoid repeated access to global memory.

Each block will perform a multiplication of a 1x8 array with an 8x8 array.

It’s clear that this approach will lead to bank conflicts. I also checked with Nsight Compute.
Can they be avoided, or is shared memory unsuitable for this scenario?

image
image1091×555 15 KB

4

bank conflicts in shared memory != shared memory unsuitable

Are these bank conflicts significant for overall application performance? Does a higher-performing alternative design variant exist?

5

well, since you are loading shared memory, if you transpose when you load, there would be no subsequent bank conflicts - you would be multiplying row-wise across threads by row-wise across threads

alternatively, a shared array with columnar access can generally be accessed without bank conflicts by adding one column.

are you using the matrix more than once? If not, I doubt use of shared is going to help.

A block strikes me as overkill for that operation. I’m not sure it merits more than a warp. I hope you are not running blocks of a single warp.

If you are only using the matrix once, I would probably suggest something like this:

#include <thrust/complex.h>

// R = VxM
const int Vl = 8; // maximum 32
using mt = thrust::complex<double>;
__global__ void k(mt *V, mt *M, mt *R){
  __shared__ mt sV[Vl];
  if (threadIdx.x < Vl) sV[threadIdx.x] = V[threadIdx.x];
  __syncwarp();
  if (threadIdx.x < Vl){
    mt ts = {0};
    for (int i = 0; i < Vl; i++)
      ts += sV[i]*M[i*Vl+threadIdx.x];
    R[threadIdx.x] = ts;}
}

This pretty much assumes you have something else useful for more than 8 threads to do, e.g. wanting to do many of the vector-matrix multiplies.

If the only thing you wanted to do in the entirety of your GPU code is this single Vector-matrix multiply, then that problem is too small to be interesting on a GPU. It’s possible in that case another realization might be “slightly faster”, like loading all the values at once, then doing some warp-shuffle gymnastics, but that type of realization would not make much sense to me at all without more context about what else (if anything) you are trying to accomplish.

6

are you using the matrix more than once

Yes, below is my kernel implementation which is MUSIC DOA estimation.
It will scan the azimuth field of view (-60~60 degree) and elevation field of view (60~95 degree).
At each pair of (azimuth, elevation), it will calculate (a^H)(UU^H)(a).

a (8x1): steering vector of certain pair of (azimuth, elevation)
H: Hermitian operator
U: noise subspace of signal
Ucov (8x8): covariance matrix of noise subspace

I merged UU^H into Ucov, which means the covariance matrix of noise subspace to simplify the multiplication.

So in each pair of (azimuth, elevation), after I calculated 8 steering vectors, I saved them into shared_steering_v.

The covariance matrix of noise subspace (8*8) also saved into shared_mem shared_noise_cov.

And both shared_steering_v and shared_noise_cov will be used more than once.

cfg.n_ant = 8;

__global__ void spectrum_scanning(double *spectrum_out, cuDoubleComplex *noiseSpace, DOAConfig cfg, int azimuth_begin, int azimuth_end, int elevation_begin, int elevation_end) {
    
    extern __shared__ cuDoubleComplex shared_mem[];

    cuDoubleComplex *shared_steering_v = shared_mem;
    cuDoubleComplex *shared_noise_cov  = shared_mem + cfg.n_ant;

    double d      = cfg.d;
    double lambda = cfg.lambda;

    int tid = threadIdx.x;
    int bid = blockIdx.x * gridDim.y + blockIdx.y;

    int row = tid / cfg.n_ant;
    int col = tid % cfg.n_ant;
    int col_majored_idx = row + col * cfg.n_ant;

    int curr_azimuth   = azimuth_begin + blockIdx.x;
    int curr_elevation = elevation_begin + blockIdx.y;

    if (tid < cfg.n_ant) {

        double azimuth   = curr_azimuth * M_PI / 180.0;
        double elevation = curr_elevation * M_PI / 180.0;
        
        double sin_phi = sin(azimuth);
        double cos_phi = cos(azimuth);
        double sin_theta = sin(elevation);
        double cos_theta = cos(elevation);

        double phaseShift = 2 * M_PI * d * 
                            (cfg.antenna[tid][0] * (cos_phi * sin_theta) + cfg.antenna[tid][1] * (sin_phi * sin_theta) + cfg.antenna[tid][2] * cos_theta) /
                            lambda;

        shared_steering_v[tid] = make_cuDoubleComplex(cos(phaseShift), sin(phaseShift));

    }
    if (tid < cfg.n_ant * cfg.n_ant) {
        shared_noise_cov[tid] = noiseSpace[tid];
    }
    __syncthreads();

    if (tid < cfg.n_ant * cfg.n_ant) {
        shared_noise_cov[col_majored_idx] = cuCmul(shared_noise_cov[col_majored_idx], cuConj(shared_steering_v[col]));
        
    }
    __syncthreads();

    // reduce axis-Y
    int blockSize = cfg.n_ant * cfg.n_ant;
    for (int stride = blockSize / 2; stride >= cfg.n_ant; stride >>= 1) {
        if (col_majored_idx < stride) {
            shared_noise_cov[col_majored_idx] = cuCadd(shared_noise_cov[col_majored_idx], shared_noise_cov[col_majored_idx + stride]);
        }
        __syncthreads();
    }

    if (tid < cfg.n_ant) {
        shared_steering_v[tid] = cuCmul(shared_noise_cov[tid], shared_steering_v[tid]);
    }
    __syncthreads();
    
    int steer_size = cfg.n_ant;
    for (int stride = steer_size / 2; stride > 0; stride >>= 1) {
        if (tid < stride) {
            shared_steering_v[tid] = cuCadd(shared_steering_v[tid], shared_steering_v[tid + stride]);
        }
        __syncthreads();
    }

    if (tid == 0) {
        spectrum_out[bid] = 10 * log10(1.0 / cuCabs(shared_steering_v[0]));
    }
}

The way I launch the kernel:

Each time I scan 11 degrees of azimuth FOV, and 36 degree of elevation FOV, so the dim of blockNum_ is (11, 36).
The thread number of each block is 64 for the size of noise subspace covariance matrix is (8*8).

    dim3 blockNum_(11, 36);
    int threadNum_ = 64;
    int sharedMem  = m_cfg.n_ant * (m_cfg.n_ant + 1) * sizeof(cuDoubleComplex);

    // #pragma unroll
    for (int azimuth_row = 0; azimuth_row < (AZIMUTH_END - AZIMUTH_BEGIN); azimuth_row += blockNum_.x) {
        int offset = azimuth_row * blockNum_.y;
        spectrum_scanning<<<blockNum_, threadNum_, sharedMem, m_streams[azimuth_row % n_stream]>>>(
            m_spectrum.data().get() + offset,
            m_noise_space_cov.data().get(),
            m_cfg,
            AZIMUTH_BEGIN + azimuth_row,
            AZIMUTH_BEGIN + azimuth_row + blockNum_.x -1,
            ELEVATION_BEGIN, ELEVATION_END
        );
    }
    
    for (int i = 0; i < n_stream; ++i) {
        cudaStreamSynchronize(m_streams[i]);
    }

Thanks for your kind help.

7

Recomendation (performance): Use sincos().

Recommendation (performance and accuracy): Remove the multiplication with M_PI from of the computation of phaseShift, then use sincospi() to compute shared_steering_v.

8

Would you mind sharing a complete example which can be compiled and run?
Which GPU are you using?
Is n_ant always set to 8, or do you also have other configurations?

9

I thought of a possibility:

configuring 64 threads per thread block will inevitably cause bank conflicts (according to Nsight Compute calculations, there are currently over 40,000).

To achieve the best runtime, will this become a trade-off between the number of threads (parallelization level) and the number of bank conflicts?

For example, I might write several versions: 64 threads per block, 32 threads per block, 16 threads per block, and 8 threads per block, and bank conflicts might decrease accordingly. Perhaps a sweet spot for runtime can be found.

@striker159
I am using RTX4060.

I can provide a compileable version, but since the test data part is somewhat troublesome, I’d like to focus primarily on the algorithmic ideas and discussion for this topic.

10

The input data for the compilable version could be random numbers to test the kernel.