CUDA's nvvp reports non-ideal memory access pattern, but bandwidth is almost peaking

Question

EDIT: new minimal working example to illustrate the question and better explanation of nvvp's outcome (following suggestions given in the comments).

So, I have crafted a "minimal" working example, which follows:

#include <cuComplex.h>
#include <iostream>

int const n = 512 * 100;

typedef float real;

template < class T >
struct my_complex {
   T x;
   T y;
};

__global__ void set( my_complex< real > * a )
{
   my_complex< real > & d = a[ blockIdx.x * 1024 + threadIdx.x ];
   d = { 1.0f, 0.0f };
}

__global__ void duplicate_whole( my_complex< real > * a )
{
   my_complex< real > & d = a[ blockIdx.x * 1024 + threadIdx.x ];
   d = { 2.0f * d.x, 2.0f * d.y };
}

__global__ void duplicate_half( real * a )
{
   real & d = a[ blockIdx.x * 1024 + threadIdx.x ];
   d *= 2.0f;
}

int main()
{
   my_complex< real > * a;
   cudaMalloc( ( void * * ) & a, sizeof( my_complex< real > ) * n * 1024 );

   set<<< n, 1024 >>>( a );
   cudaDeviceSynchronize();
   duplicate_whole<<< n, 1024 >>>( a );
   cudaDeviceSynchronize();
   duplicate_half<<< 2 * n, 1024 >>>( reinterpret_cast< real * >( a ) );
   cudaDeviceSynchronize();

   my_complex< real > * a_h = new my_complex< real >[ n * 1024 ];
   cudaMemcpy( a_h, a, sizeof( my_complex< real > ) * n * 1024, cudaMemcpyDeviceToHost );

   std::cout << "( " << a_h[ 0 ].x << ", " << a_h[ 0 ].y << " )" << '\t' << "( " << a_h[ n * 1024 - 1 ].x << ", " << a_h[ n * 1024 - 1 ].y << " )"  << std::endl;

   return 0;
}

When I compile and run the above code, kernels duplicate_whole and duplicate_half take just about the same time to run.

However, when I analyze the kernels using nvvp I get different reports for each of the kernels in the following sense. For kernel duplicate_whole, nvvp warns me that at line 23 (d = { 2.0f * d.x, 2.0f * d.y };) the kernel is performing

Global Load L2 Transaction/Access = 8, Ideal Transaction/Access = 4

I agree that I am loading 8 byte words. What I do not understand is why 4 bytes is the ideal word size. In special, there is no performance difference between the kernels.

I suppose that there must be circumstances where this global store access pattern could cause performance degradation. What are these?

And why is that I do not get a performance hit?

I hope that this edit has clarified some unclear points.

+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

I'll start wit some kernel code to exemplify my question, which will follow below

template < class data_t >
__global__ void chirp_factors_multiply( std::complex< data_t > const * chirp_factors,
                                        std::complex< data_t > * data,
                                        int M,
                                        int row_length,
                                        int b,
                                        int i_0
                                        )
{
#ifndef CUGALE_MUL_SHUFFLE
    // Output array length:
    int plane_area = row_length * M;
    // Process element:
    int i = blockIdx.x * row_length + threadIdx.x + i_0;
    my_complex< data_t > const chirp_factor = ref_complex( chirp_factors[ i ] );
    my_complex< data_t > datum;
    my_complex< data_t > datum_new;

    for ( int i_b = 0; i_b < b; ++ i_b )
    {
        my_complex< data_t > & ref_datum = ref_complex( data[ i_b * plane_area + i ] );
        datum = ref_datum;
        datum_new.x = datum.x * chirp_factor.x - datum.y * chirp_factor.y;
        datum_new.y = datum.x * chirp_factor.y + datum.y * chirp_factor.x;
        ref_datum = datum_new;
    }
#else
    // Output array length:
    int plane_area = row_length * M;
    // Element to process:
    int i = blockIdx.x * row_length + ( threadIdx.x + i_0 ) / 2;
    my_complex< data_t > const chirp_factor = ref_complex( chirp_factors[ i ] );

    // Real and imaginary part of datum (not respectively for odd threads):
    data_t datum_a;
    data_t datum_b;

    // Even TIDs will read data in regular order, odd TIDs will read data in inverted order:
    int parity = ( threadIdx.x % 2 );
    int shuffle_dir = 1 - 2 * parity;
    int inwarp_tid = threadIdx.x % warpSize;

    for ( int i_b = 0; i_b < b; ++ i_b )
    {
        int data_idx = i_b * plane_area + i;
        datum_a = reinterpret_cast< data_t * >( data + data_idx )[ parity ];
        datum_b = __shfl_sync( 0xFFFFFFFF, datum_a, inwarp_tid + shuffle_dir, warpSize );

        // Even TIDs compute real part, odd TIDs compute imaginary part:
        reinterpret_cast< data_t * >( data + data_idx )[ parity ] = datum_a * chirp_factor.x - shuffle_dir * datum_b * chirp_factor.y;
    }
#endif // #ifndef CUGALE_MUL_SHUFFLE
}

Let us consider the case where data_t is float, which is memory bandwidth limited. As it can be seen above, there are two versions of the kernel, one which reads/writes 8 bytes (a whole complex number) per thread and another which reads/writes 4 bytes per thread and then shuffles the results so the complex product is computed correctly.

The reason why I have written the version using shuffle is because nvvp insisted that reading 8 bytes per thread was not the best idea because this memory access pattern would be inefficient. This is the case even though in both systems tested (GTX 1050 and GTX Titan Xp) memory bandwidth was very close to theoretical maximum.

Surely enough I knew that no improvement was likely to happen, and this was indeed the case: both kernels take pretty much the same time to run. So, my question is the following:

Why is that nvvp reports that reading 8 bytes would be less efficient than reading 4 bytes per thread? In which circumstances would that be the case?

As a side note, single precision is more important to me, but double is useful in some cases too. Interestingly enough, in the case where data_t is double, there is no execution time difference too between the two kernel versions, even though in this case the kernel is compute bound and the shuffle version performs some more flops than the original version.

Note: the kernels are applied to a row_length * M * b dataset (b images with row_length columns and M lines) and the chirp_factor array is row_length * M. Both kernels run perfecly fine (I can edit the question to show you the calls to both versions if you have doubts about it).

Answer 1

The issue here has to do with how the compiler is processing your code. nvvp is merely dutifully reporting what is happening when you run your code.

If you use the cuobjdump -sass tool on your executable, you will discover that the duplicate_whole routine is doing two 4-byte loads and two 4-byte stores. This is not optimal, partly becuase there is a stride in each load and store (each load and store touches alternate elements in memory).

The reason for this is that the compiler does not know the alignment of your my_complex struct. Your struct would be legal for use in situations that would prevent the compiler from generating a (legal) 8-byte load. As discussed here we can fix this by informing the compiler that we only intend to use the struct in alignment scenarios where a CUDA 8-byte load is legal (i.e. it is "naturally aligned"). The modification to your struct looks like this:

template < class T >
struct  __align__(8) my_complex {
   T x;
   T y;
};

With that change to your code, the compiler generates 8-byte loads for the duplicate_whole kernel, and you should see a different report from the profiler. You should use this sort of decoration only when you understand what it means and are willing to enter into a contract with the compiler that you will ensure this is the case. If you do something unusual, like unusual pointer casting, you can violate your end of the bargain and generate a machine fault.

The reason you don't see much performance difference almost certainly has to do with CUDA load/store behavior and the GPU caches

When you do a strided load, the GPU loads an entire cacheline anyway, even though (in this case) you only need half the elements (the real elements) for that particular load operation. However you need the other half of the elements (the imaginary elements) anyway; they will be loaded on the next instruction, and this instruction most likely hits in the cache, due to the previous load.

On a strided store in this case, writing strided elements in one instruction and the alternate elements in the next instruction will end up using one of the caches as a "coalescing buffer". This isn't coalescing in the typical sense used in CUDA terminology; that sort of coalescing only applies to a single instruction. However the cache "coalescing buffer" behavior allows it to "accumulate" multiple writes to an already-resident line, before that line gets written out or evicted. This is approximately equivalent to "write-back" cache behavior.