How do I coalesce global reads in chunks with size the same as block size?

user_1446510 0
Viewed 304

I have a long array of float2's in global memory from which each block of my kernel reads 512 continuous values in each chunk (there is no discernible pattern to the specific chunk any given block accesses however). So I was thinking that I should be able to coalesce those 512 float2 loads from global memory. In an attempt to do that, I load the pointer to the first element out of the 512 into shared memory and then have each thread increment the pointer by its threadID to read the float2 from that address. However, this doesn't seem to be coalescing the reads when I check with Nsight Compute. I've marked the lines below where it says I have uncoalesced global accesses.

Along with sequential access, the best practices guide mentions alignment requirement but I am using a vector type, which the guide says are automatically aligned (in addition to them corresponding to 8 bytes per value). So I'm not sure why my reads aren't coalesced. The ratio of real to ideal global L2 sectors is around 1.12 for the read at the top and 1.09 for the bottom, which I'm assuming aren't that bad, but I still want to optimize those reads if possible.

One slight complication in my code that I didn't mention above is that each of those 512 reads are done nsegs times in a loop but I'm not sure if that would affect the read patterns. Another strange thing is that, in my real full-code, Nsight says I have uncoalesced read in the second read (where I increase the pointer by threadID and read that address) but does not mark the read from global to shared memory as uncoalesced. Unlike that, in the example code below, Nsight marks both of them as uncoalesced reads. Also, my full-code uses the __shfl_down_sync warpReduce method (called by blockReduceSum) but I'm doing atomic sums below to make the example code shorter.

#include <stdio.h>

#define gpuErrchk(ans) { gpuAssert((ans), __FILE__, __LINE__); }

inline void gpuAssert(cudaError_t code, const char *file, int line, bool abort = true) {
   if (code != cudaSuccess) {
      fprintf(stderr, "GPUassert: %s %s %d\n", cudaGetErrorString(code), file, line);
      if (abort) exit(code);
   }
}

#define nbins 512
#define nsegs 340
#define ntemplates 100

__device__ float2 * d;

__global__ void kernel(float * d_all_sums, int * d_template_indices) {

    float final = 0.0;

    // blockIdx.x gives template index, threadIdx.x normally gets bin index but use it to mean segment index for loading global memory to shared.

    __shared__ float2 * power_at_0_pointers[nsegs];

    if (threadIdx.x < nsegs)
        power_at_0_pointers[threadIdx.x] = &d[ __ldg(&d_template_indices[blockIdx.x * nsegs + threadIdx.x]) ];  // Uncoalesced here. Real to ideal ratio ~ 1.12.

    __syncthreads();

    __shared__ float power_sum;

    for (int i = 0; i < nsegs; i++) {
        __shared__ float powers[nbins];

        float2 * pow_first_bin = power_at_0_pointers[i];
        float2 input_power_c = *(pow_first_bin + threadIdx.x);       // Uncoalesced here. Real to ideal ratio ~ 1.09.
        
        float power = input_power_c.x * input_power_c.x + input_power_c.y * input_power_c.y;

        power = (2 * power - 1.0) / 2.0;
        powers[threadIdx.x] = power;

        atomicAdd(&power_sum, powers[threadIdx.x]);

        final += power / power_sum;
    }

    if (threadIdx.x == 0)
        d_all_sums[blockIdx.x] = final;
}

int random(int min, int max){
   return min + rand() / (RAND_MAX / (max - min + 1) + 1);
}

int main(){

    float2 * h;
    float * h_all_sums;
    float2 * d_ptr;
    float * d_all_sums;
    int * h_template_indices, * d_template_indices;

    h_template_indices = (int *) malloc(ntemplates * nsegs * sizeof(int));
    h = (float2 *) malloc(nbins * nsegs * ntemplates * sizeof(float2));
    h_all_sums = (float *) malloc(ntemplates * sizeof(float));
    memset(h_all_sums, 0.0, ntemplates);

    for (int k = 0; k < ntemplates; k++) {
        for (int i = 0; i < nsegs; i++) {
            h_template_indices[k * nsegs + i] = random(0, ntemplates * nsegs - nbins);
            for (int j = 0; j < nbins; j++)
                h[k * nbins * nsegs + i * nbins + j] = make_float2(100 * (float) rand() / (float)(RAND_MAX), 100 * (float) rand() / (float)(RAND_MAX));
        }
    }

    gpuErrchk( cudaMalloc((void**) &d_ptr, nbins * nsegs * ntemplates * sizeof(float2)) );
    gpuErrchk( cudaMemcpy(d_ptr, h, nbins * nsegs * ntemplates * sizeof(float2), cudaMemcpyHostToDevice) );
    gpuErrchk( cudaMemcpyToSymbol(d, &d_ptr, sizeof(float2*)) );

    gpuErrchk( cudaMalloc( (void**) &d_all_sums, ntemplates * sizeof(float) ) );

    gpuErrchk( cudaMalloc((void**) &d_template_indices, nsegs * ntemplates * sizeof(int)) );
    gpuErrchk( cudaMemcpy(d_template_indices, h_template_indices, nsegs * ntemplates * sizeof(int), cudaMemcpyHostToDevice) );

    kernel<<<ntemplates, nbins>>>(d_all_sums, d_template_indices);

    gpuErrchk( cudaMemcpy(h_all_sums, d_all_sums, ntemplates * sizeof(float), cudaMemcpyDeviceToHost) );

    FILE *f = fopen("test_output.txt", "w");
    if (f != NULL) {
        for (int k = 0; k < ntemplates; k++)
            fprintf(f, "k = %d; power = %f.\n", k, h_all_sums[k]);
    }

    fclose( f );

    gpuErrchk( cudaFree(d_ptr) );
    gpuErrchk( cudaFree(d_all_sums) );
    gpuErrchk( cudaFree(d_template_indices) );

    free( h );
    free( h_all_sums );
    free( h_template_indices );

    gpuErrchk( cudaPeekAtLastError() );
    gpuErrchk( cudaDeviceSynchronize() );

    printf("All done.\n");
}

EDIT: I've included the result of Nsight Compute about coalescence below. I've also included code that matches my full code more closely, in terms of how block reduction is done (shuffle instead of atomics). With that change, now Nsight says I have uncoalesced access only at one stop (another one in a library but that's obviously out of my control). The difference in summing/reduction method, for some reason, seemed to have made the read from global to shared memory to coalesced?

Nsight Compute report of global memory coalescence.

Source and SASS for uncoalesced access

#include <stdio.h>

#define gpuErrchk(ans) { gpuAssert((ans), __FILE__, __LINE__); }

inline void gpuAssert(cudaError_t code, const char *file, int line, bool abort = true) {
   if (code != cudaSuccess) {
      fprintf(stderr, "GPUassert: %s %s %d\n", cudaGetErrorString(code), file, line);
      if (abort) exit(code);
   }
}

__inline__ __device__ float warpReduceSum(float val) {
  for (int offset = warpSize/2; offset > 0; offset /= 2)
    val += __shfl_down_sync(0xffffffff, val, offset);
  return val;
}

__inline__ __device__ float blockReduceSum(float val) {
  static __shared__ float shared[32];
  int lane = threadIdx.x % warpSize;
  int wid = threadIdx.x / warpSize;
  val = warpReduceSum(val);
  if (lane == 0) shared[wid] = val;
  __syncthreads();
  val = (threadIdx.x < blockDim.x / warpSize) ? shared[lane] : float(0.0);
  if (wid == 0) val = warpReduceSum(val);
  return val;
}

#define nbins 512
#define nsegs 340
#define ntemplates 100

__device__ float2 * d;

__global__ void kernel(float * d_all_sums, int * d_template_indices) {

    float final = 0.0;

    // blockIdx.x gives template index, threadIdx.x normally gets bin index but use it to mean segment index for loading global memory to shared.

    __shared__ float2 * power_at_0_pointers[nsegs];

    if (threadIdx.x < nsegs)
        power_at_0_pointers[threadIdx.x] = &d[ __ldg(&d_template_indices[blockIdx.x * nsegs + threadIdx.x]) ];

    __syncthreads();

    for (int i = 0; i < nsegs; i++) {
        __shared__ float powers[nbins];

        float2 * pow_first_bin = power_at_0_pointers[i];
        float2 input_power_c = *(pow_first_bin + threadIdx.x);       // Uncoalesced here.
        
        float power = input_power_c.x * input_power_c.x + input_power_c.y * input_power_c.y;

        power = (2 * power - 1.0) / 2.0;
        powers[threadIdx.x] = 1.0;

        float power_sum = blockReduceSum(powers[threadIdx.x]);

        final += power / power_sum;
    }

    if (threadIdx.x == 0)
        d_all_sums[blockIdx.x] = final;
}

int random(int min, int max){
   return min + rand() / (RAND_MAX / (max - min + 1) + 1);
}

int main(){

    float2 * h;
    float * h_all_sums;
    float2 * d_ptr;
    float * d_all_sums;
    int * h_template_indices, * d_template_indices;

    h_template_indices = (int *) malloc(ntemplates * nsegs * sizeof(int));
    h = (float2 *) malloc(nbins * nsegs * ntemplates * sizeof(float2));
    h_all_sums = (float *) malloc(ntemplates * sizeof(float));
    memset(h_all_sums, 0.0, ntemplates);

    for (int k = 0; k < ntemplates; k++) {
        for (int i = 0; i < nsegs; i++) {
            h_template_indices[k * nsegs + i] = random(0, ntemplates * nsegs - nbins);
            for (int j = 0; j < nbins; j++)
                h[k * nbins * nsegs + i * nbins + j] = make_float2(100 * (float) rand() / (float)(RAND_MAX), 100 * (float) rand() / (float)(RAND_MAX));
        }
    }

    gpuErrchk( cudaMalloc((void**) &d_ptr, nbins * nsegs * ntemplates * sizeof(float2)) );
    gpuErrchk( cudaMemcpy(d_ptr, h, nbins * nsegs * ntemplates * sizeof(float2), cudaMemcpyHostToDevice) );
    gpuErrchk( cudaMemcpyToSymbol(d, &d_ptr, sizeof(float2*)) );

    gpuErrchk( cudaMalloc( (void**) &d_all_sums, ntemplates * sizeof(float) ) );

    gpuErrchk( cudaMalloc((void**) &d_template_indices, nsegs * ntemplates * sizeof(int)) );
    gpuErrchk( cudaMemcpy(d_template_indices, h_template_indices, nsegs * ntemplates * sizeof(int), cudaMemcpyHostToDevice) );

    kernel<<<ntemplates, nbins>>>(d_all_sums, d_template_indices);

    gpuErrchk( cudaMemcpy(h_all_sums, d_all_sums, ntemplates * sizeof(float), cudaMemcpyDeviceToHost) );

    FILE *f = fopen("test_output.txt", "w");
    if (f != NULL) {
        for (int k = 0; k < ntemplates; k++)
            fprintf(f, "k = %d; power = %f.\n", k, h_all_sums[k]);
    }

    fclose( f );

    gpuErrchk( cudaFree(d_ptr) );
    gpuErrchk( cudaFree(d_all_sums) );
    gpuErrchk( cudaFree(d_template_indices) );

    free( h );
    free( h_all_sums );
    free( h_template_indices );

    gpuErrchk( cudaPeekAtLastError() );
    gpuErrchk( cudaDeviceSynchronize() );

    printf("All done.\n");
}

EDIT2: Added memory chart and tables as requested.

Nsight memory chart

Nsight memory table

cuda
1 Answers

So I'm not sure why my reads aren't coalesced.

I believe you are conflating two different things. You haven't actually shown the nsight compute output that you are using to make your judgements, but I believe you are mixing up coalescing with efficiency. This statement:

The ratio of real to ideal global L2 sectors is around 1.12 for the read at the top and 1.09 for the bottom,

is really about efficiency.

Let's review some definitions:

  • coalescing - This refers to the ability to group requests together based on locality of the addresses those requests refer to. When the addresses of requests from each thread in a warp can be grouped together, the GPU memory controller will coalesce those requests to create some smaller number of transactions, to service those requests

  • efficiency - for the usage here, we are referring to the bytes necessary to service those requests, divided by the bytes that were actually retrieved to service those requests

One way that efficiency can be reduced is if coalescing is reduced. However there are other ways that efficiency can be reduced, and your code is exercising one of those, in terms of alignment.

When the memory controller needs data (lets say from DRAM, but caches also have granular subdivisions) it is not possible to request an arbitrary number of bytes. The DRAM memory (or a cache) serves up data in groups called segments (caches: lines or sectors). The memory controller will request one or more segments from DRAM memory to serve the needs of particular read requests.

What happens when we have unaligned reads? Here's a pictorial example:

segments:      |  1  |  2  |  3  |  4  |  5  | ...
aligned read:  ^^^^^^^^^^^^^^^^^^^^^^^^

In the case of the aligned read, we see that it is both coalesced and aligned. That means that all the data that is needed is adjacent to each other (so it coalesces nicely) and also that the start and end points of the data region requested align with segment boundaries. Therefore the memory controller will retrieve 4 segments to service the above aligned read request, and that is exactly what would be required anyway (in terms of bytes), so the efficiency is 100% or a 1.0 ratio.

segments:      |  1  |  2  |  3  |  4  |  5  | ...
unaligned read:  ^^^^^^^^^^^^^^^^^^^^^^^^

In the case of the unaligned read, the amount of bytes requested is the same as the aligned read, but now since we don't align with segment boundaries, the memory controller must retrieve segments 1-5 (a total of 5 segments rather than 4) to service the needs of the read request. This is still "highly" coalesced - there is only one address group, and all bytes requested are adjacent to each other. But the efficiency is reduced because now the memory controller must retrieve 5 segments, instead of 4, to serve up the same number of bytes requested by the warp request. The ratio of actual to necessary here would be 5/4 or 1.25.

Note that the above discussion does not take into account what may already be in the caches, and neither does the profiler when making these kinds of measurement. The above discussion is equally applicable to cache behavior/measurements, however, simply replacing "lines" or "sectors" for "segments".

Your line of code here:

for (int i = 0; i < nsegs; i++) {
    __shared__ float powers[nbins];

    float2 * pow_first_bin = power_at_0_pointers[i];
    float2 input_power_c = *(pow_first_bin + threadIdx.x);       // Uncoalesced here

is experiencing this kind of "efficiency loss" when measured by the profiler. From a coalescing perspective, I would say that that code coalesces "nicely" because the + threadIdx.x in the index construction pretty much guarantees that: adjacent threads (in x) will pull from adjacent locations in memory. That is the recipe for "ideal" coalescing. However, the for-loop offsets the pointers by an index value i which means the read "pattern" (i.e. the coalesced group of addresses) will "walk" accross the pattern of L2 sectors (or memory segments) as the loop progresses. Except for those locations where the i index offset is actually aligned with a sector/segment boundary, the read request at that particular i offset will require additional sector(s) from the L2 cache.

So is this a problem?

I would generally say no, for a few reasons:

  1. This is a commonly occurring pattern, and difficult to refactor.
  2. There is still what I would refer to as "nearly ideal" coalescing.
  3. Concerns about efficiency mostly ignore actual cache behavior.

Let's delve into item 3 above more deeply.

The general purpose of a cache is to reduce the penalty associated with going to the next level in the memory hierarchy to retrieve data. For the L2 cache, the next level in the memory hierarchy is DRAM memory (generally speaking). For the "walking" pattern that you have, that pattern will populate the L2 cache with the needed data pretty quickly, so that subsequent iterations in i generally don't have to retrieve much (if anything) from DRAM memory, in order to service that loop as it proceeds. Therefore the L2 cache is doing its job. It is mostly insulating that loop from the cost of DRAM fetches, as it proceeds. The fact that you sometimes have to retrieve extra sectors from the L2 as the loop proceeds is a far less pressing concern than if you would have to go to DRAM for each loop iteration, for all the data.

Along with sequential access, the best practices guide mentions alignment requirement but I am using a vector type, which the guide says are automatically aligned.

vector types are not automatically aligned, when the alignment we are referring to is sector/segment alignment. That should hopefully be obvious from the previous discussion.

In order for vector types to be useful, they must be naturally aligned. This means that the address of the element itself must line up to a memory boundary that is a whole-number multiple of the byte size of the vector type. A vector type that consists of 8 bytes must be located at address 0, or 8, or 16, or 24, or 32, or 40, etc. This is not the same kind of alignment being discussed here and elsewhere.

My opinion is you're going to have difficulty "optimizing" this code if you simply focus on the line the profiler is reporting for efficiency loss. And asking why there is efficiency loss there (such as what I have tried to describe here) in my opinion will not be a very fruitful avenue to optimizing this particular code. Instead I would try to refactor the code as a whole, in the context of what you are trying to do, none of which is clear to me.

The two areas I would focus on without knowing the higher level purpose of this code algorithmically, would be:

  1. Seek to get rid of atomic usage
  2. Wherever there is data reuse, see if shared memory usage can provide a benefit.

Getting rid of atomic usage often requires refactoring the algorithm significantly. Right now the code focuses on processing data from an input perspective. One possible approach to get rid of atomics is to focus on asking "what is the final result in each location" and then write code that produces that, perhaps one thread per output or one warp or block per output point. By restructuring code this way, you can usually get rid of atomic usage. Another alternative is to look at canonical parallel reduction methods to get rid of atomics.

Shared memory refactoring in the presence of data reuse across threads is one of the most common CUDA tutorial topics. You'll find many treatments of that.

For an orderly introduction to some of the concepts I have discussed here, you may wish to review (at least the first 4 sessions of) this CUDA tutorial series.

user_16959600
Related