HPC - Lecture 12 - Coalescence and Shared Memory GPU

---

Coalescence

Transactions to and from global memory are governed by coalescing rules. The goal of coalescence is to optimize read and write operations in global memory by improving memory access performance.

Definition of Coalescing

Coalescing refers to a set of conditions that allow multiple memory accesses to be merged into a single transaction.

This memory access pattern makes it possible to achieve very high performance when using global memory

Warp and Half-Warp Organization

Kernel memory accesses are said to be coalesced when threads with consecutive identifiers access contiguous memory locations using the same instruction

This follows from the way threads are organized into warps.

• Warp: Groups of 32 threads all executing the same instruction
• Hardware can combine memory accesses when threads with consecutive identifiers access contiguous memory locations using the same instruction.
• Threads in a warp execute the same instruction at the same time. In this way, the hardware can combine their memory accesses into a single transaction.
• Half-Warp: Memory accesses are specifically grouped into 16-thread units for transaction merging.

Guidelines for Coalescing

To achieve high performance, memory accesses should follow these rules:

• Contiguity: Accesses are most efficient when they are contiguous.
• Alignment: The starting address of a memory region must be a multiple of the region’s size (granularity).
• One-to-One Correspondence: The -th thread of a half-warp should access the -th element of a block.
• Thread Participation: Accesses are still coalesced even if some threads do not participate in the operation (inactive threads).

To achieve contiguity:

int indexRow=threadIdx.y + blockIdx.y*blockDim.y;
int indexCol=threadIdx.x + blockIdx.x*blockDim.x;
array_name[indexRow*col+indexCol]= value

An example of succesfull coalescence:

An example of failed coalescence:

Failure Conditions

Coalescing fails when:

• Accesses are non-sequential.
• Access ranges overlap.
• Accesses start at an address that is not a multiple of the required granularity (e.g., starting at 132 instead of 128 for 64-byte granularity).

Exercise: Outer Product of Two Vectors

• Computing for .
• Input: array A of size M and array B of size N, integers M,N
• Output: matrix C of size MxN

Kernel Configuration:

• Use a two-dimensional structure for blocks and threads.
• indexRow = threadIdx.y + blockIdx.y * blockDim.y (Iterates over rows of and elements of ).
• indexCol = threadIdx.x + blockIdx.x * blockDim.x (Iterates over columns of and elements of ).

Memory Allocation with Pitch: For 2D matrices, cudaMallocPitch is used to ensure hardware-friendly alignment:

cudaMallocPitch((void**)&C_device, &pitch, N * sizeof(int), M);

In the kernel, access the element using:

c[indexRow * pitch + indexCol] = a[indexRow] * b[indexCol];

#include ...
int main(int argn, char * argv[])
{
	//variables declaration...
	// test on main arguments ...
	// variables initialization ...
	// host memory allocation(2D-1D)...
	A_host=(int*)malloc(M*sizeof(int));
	B_host=(int*)malloc(N*sizeof(int));
	C_host=(int*)malloc(M*N*sizeof(int));
	copy=(int*)malloc(M*N*sizeof(int));
	// device memory allocation
	cudaMalloc((void**)&A_device, M*sizeof(int));
	cudaMalloc((void**)&B_device, N*sizeof(int));
	cudaMallocPitch((void**)&C_device,&pitch,N*sizeof(int), M);
	
	//data initialization on the host
	//data transfer from host to device
	cudaMemcpy(A_device, A_host, M*sizeof(int),
	cudaMemcpyHostToDevice);
	cudaMemcpy(B_device, B_host, N*sizeof(int),
	cudaMemcpyHostToDevice);
	//kernel launch
	Outer_ProductGPU<<<nBlocks, nThreadsPerBlock>>>
	(A_device,B_device,C_device,pitch/sizeof(int),M,N);
	//results transfer from device to host
	cudaMemcpy2D(copy, N*sizeof(int), C_device,pitch,
	N*sizeof(int), M, cudaMemcpyDeviceToHost);
	//call to the serial version………
	//print results………
	//correctness tests
	//memory de-allocation
	//exit
}
 

Implement the serial algorithm and its parallel version (kernel)

void OuterProductCPU (int*a, int*b, int*c,
int m, int n)
{
	int i,j;
	for(i=0;i<m;i++)
		for(j=0;j<n;j++)
			c[i*n+j]=a[i]*b[j];
}
__global__ void Outer_ProductGPU(int*a, int*b, int*c, int pitch, int m, int n)
{
	int indexRow, indexCol;
	indexRow=threadIdx.y+blockIdx.y*blockDim.y;
	indexCol=threadIdx.x+blockIdx.x*blockDim.x;
	if(indexRow<m && indexCol<n)
		c[indexRow*pitch+indexCol] = a[indexRow]*b[indexCol];
}

---

Shared Memory (SM)

Shared memory is on-chip memory with much lower latency than global memory. It is shared among all threads within the same block.

Stencil of a Vector

A stencil operation computes an output element as a sum of a neighborhood of radius around the corresponding input element:

Pseudocode:

void stencilCPU(int *in, int *out, int n, int k)
{
 int i,j;
 for(i=0;i<n;i++)
	 for(j=-k;j<=k;j++)
		out[i]+=(i+j)<0||(i+j)>=n?0:in[i+j];
}

The i-th element of array B is them sum of the elements in the neighborhood of raiud around the -th element of array .

For example consider a vector with and radius .

and so on.

Approach 1: Without Shared Memory

Let’s compute the stencil of an array of size using blocks each composed of threads.

Each thread works on:

• the -th element identified by threadIdx.x + blockIdx.x*blockDim.x
• and its neighborhood of radius around the element assigned to it.
• Pay attention ! to the out-of-bonds errors for the first and last elements of the input array.

__global__ void stencilWoutSMGPU(int *in, int *out,
int n, int radius)
{
	 int index=threadIdx.x+blockIdx.x*blockDim.x;
	 int i;
	 int value=0;
	 for(i=-radius;i<=radius;i++)
		value+=(index+i)<0||(index+i)>=n?0: in[index+i];
	 out[index]=value;
}
 
// kernel launch
stencilWoutSMGPU <<<nBlocks, nThreadsPerBlock>>> (in_device, out_device_WoutSM, N,radius);

Memory Accesses: Each thread performs reads and 1 write.

• For , this results in 8 global memory accesses per thread.

Approach 2: With Shared Memory

Consider a number of threads of .

Each block shares elements in shared memory (where is the number of threads per block). Each thread operates on the neighborhood around its assigned element, first computing the sum, then storing the result into the output array. During this work, it accesses only shared memory.

• Loading Phase: Each thread loads its central element into SM. The first threads also load the “halo” elements (previous and following ).
• Computational Phase: Threads access only shared memory to compute the sum.

Global Memory Accesses: This approach requires only 4 global memory accesses (3 reads to load SM, 1 write for output) regardless of the value of .

Consider this other examle where the first threads of the block load the previous elements and the following elements into shared memory, with respect to those loaded in the previos step. In the case of out-of-bonds access, the value inserted is zero.

At this step:

__global__ void stencilWithSMGPU(int *in, int*out,
int n, int radius)
{
	extern __shared__ int shMem[]; 
}
 

• extern indicates that the size of the shared memory variable is specified as a parameter when the kernel is launched.

int ShMemSize=(nThreadsPerBlock.x+2*radius)*sizeof(int);
 
stencilWithSMGPU <<<nBlocks, nThreadsPerBlock, ShMemSize>>> (in_device, out_device_conSM, N, radius);

Shared memory allocation - Kernel parameter

If the size of shared memory is already known when the kernel is launched, it is possible to allocate it directly in the kernel launch configuration by specifying a third parameter in the triple angle brackets, indicating the size in bytes of the shared memory you want to use.

//kernel calling
nameKernel<<<nBlocks, nThreadsPerBlock, sizeShMemByte>>> (parameters)

Alternatively, shared memory can be allocated directly inside the kernel at compile-time, but only statically i.e:

__shared__ int var[500];

Kernel Pattern for Shared Memory (Stencil Example)

__global__ void stencilWithSMGPU(int *in, int *out, int n, int radius) {
    extern __shared__ int shMem[];
    int global_idx = threadIdx.x + blockIdx.x * blockDim.x;
    int local_idx = threadIdx.x + radius;
 
    // 1. Copy central element
    shMem[local_idx] = in[global_idx];
 
    // 2. Copy halo elements (previous and following k)
    if (threadIdx.x < radius) {
        shMem[local_idx - radius] = (global_idx - radius < 0) ? 0 : in[global_idx - radius];
        shMem[local_idx + blockDim.x] = (global_idx + blockDim.x >= n) ? 0 : in[global_idx + blockDim.x];
    }
 
    // 3. Synchronize
    __syncthreads();
 
    // 4. Computation
    int value = 0;
    for (int i = -radius; i <= radius; i++) {
        value += shMem[local_idx + i];
    }
 
    // 5. Write to Global Memory
    out[global_idx] = value;
}

Technical Implementation in CUDA

Shared Memory Allocation

1. Dynamic Allocation: Specified at runtime during kernel launch.
   • Kernel: extern __shared__ int shMem[];
   • Launch: kernel<<<nBlocks, nThreads, sizeInBytes>>>(...);
2. Static Allocation: Specified at compile-time.
   • Kernel: __shared__ int var[500];

Thread Synchronization

To prevent race conditions in shared memory, synchronization barriers are required.

CUDA provides two functions for thread synchronization:

• __syncthreads(): Used inside the kernel to synchronize all threads within the same block. Essential after loading data into SM before processing it.
• cudaDeviceSynchronize(): Used on the Host side to wait for the kernel to complete before accessing results.

//finally the computational step…
int value = 0;
 for(i = -radius; i<=radius; i++ )
 value += shMem[local_idx + i];
// … e then the writing procedure in GM
 out[global_idx] = value;

How many accesses to global memory does each thread perform? 1 write to store the output value and 3 reads for uploading data in the shared memory.

Without SM: reads 2*k +1 central element + k before and after k. Writes 1.

• For k=3, 2*3 +1 + 1 = 6+2 = 8

With SM: reads 3 to upload data in the SM, writes 1. For any value of k it is always 4

Efficiency

This strategy is efficient only if the array size N and the number of blocks M satisfy mod(M,N)==0

Dot Product of Two Vectors

• Computing the scalar .
• Input: array A of size N, array B of size M, integers N and M
• Output: scalar k

The dot product of two vectors is the sum of the elements of the vector obtained by pointwise multiplication of the two vectors.

void DotProductArrayCPU(int*a, int *b,int N,int*ris)
{
	 int i;
	 for(i=0;i<N;i++)
		 *ris+=(a[i]*b[i]);
}

Without Shared Memory

Let us denote =number of blocks, =number of threads per block.

We use threads to process an array of size .

Steps:

1. Each thread computes a partial sum of a portion of the vectors.
2. The partial sums computed by each thread are written to a supporting array in Global Memory.
3. The Host (CPU) performs a final serial loop to sum all partial results.

• Bottleneck: The CPU must process elements (Total Blocks Threads per Block).

// To simplify the main, we use a function that handles memory allocation and kernel parameter definitions, then launches the kernel.
 
void Kconfiguration_DotProductArrayWoutSMGPU (int* in1, int* in2, int nBlocks,int nThreadsPerBlock, int N, int *res)
{
	 int *in1_device, *in2_device;
	 int *PartialSums_host, *PartialSums_device;
	 int size=N*sizeof(int);
	 int k= N/(nThreadsPerBlock*nBlocks);
	 
	 cudaMalloc((void**)&in1_device, size);
	 cudaMalloc((void**)&in2_device, size);
	 cudaMalloc((void**)&PartialSums_device,
	 nBlocks*nThreadsPerBlock*sizeof(int));
	 
	 PartialSums_host=(int*)malloc(nBlocks*nThreadsPerBlock*sizeof(int));
	 
 cudaMemcpy(in1_device, in1, size, cudaMemcpyHostToDevice);
 cudaMemcpy(in2_device, in2, size, cudaMemcpyHostToDevice);
 cudaMemset(PartialSums_device, 0, nThreadsPerBlock*nBlocks*sizeof(int));
 DotProductArrayWoutSMGPU<<<nBlocks, nThreadsPerBlock>>>(in1_device, in2_device,PartialSums_device, k,N);
 cudaDeviceSynchronize();
 
 cudaMemcpy(PartialSums_host, PartialSums_device, nBlocks*nThreadsPerBlock*sizeof(int), cudaMemcpyDeviceToHost);
 
 for(int i=0;i<nThreadsPerBlock*nBlocks;i++)
	 *ris+=PartialSums_host[i];
 
 free(PartialSums_host);
 cudaFree(in1_device);
 cudaFree(in2_device);
 cudaFree(PartialSums_device);
}
 

The main function only contains input reading, checks, and output printing.

__global__ void DotProductArrayWoutSMGPU (int *in1, int *in2, int *PartialSums, int k, int N)
{
	int i;
	int PartialSumsIndex = threadIdx.x+blockIdx.x*blockDim.x;
	int inIndex=(threadIdx.x+blockIdx.x*blockDim.x)*k;
	if(inIndex<N)
		for(i=0;i<k;i++)
			if(inIndex+i<N)
				PartialSums[PartialSumsIndex]+= in1[inIndex+i]*in2[inIndex+i];
 else
 break;
}

With Shared Memory

It is possible to modify the previous code in order to exploit the low latency of the Shared Memory and to reduce the number of elements to sum into the host (i.e. by the CPU).

The strategy based on dividing the work into portions and computing partial sums remains the same, but here the M threads of the same block can write their partial sums to an array stored in shared memory.

1. Threads in a block compute partial sums and write them into an array in Shared Memory.
2. Thread 0 of each block sums the values in SM and writes a single block-level result to global memory.
3. The Host only needs to process (number of blocks) partial sums (instead of )

Kernel configuration

void Kconfiguration_DotProductArrayWithSMGPU (int* in1, int* in2, int nBlocks, int nThreadsPerBlock, int stride, int N, int *res)
{
 int *in1_device, *in2_device;
 int *PartialSums_host, *PartialSums_device;
 int size=N*sizeof(int);
 
 cudaMalloc((void**)&in1_device, size);
 cudaMalloc((void**)&in2_device, size);
 cudaMalloc((void**)&PartialSums_device,nBlocks*sizeof(int));
 
 PartialSums_host=(int*)malloc(nBlocks*sizeof(int));
 
 cudaMemcpy(in1_device, in1, size, cudaMemcpyHostToDevice);
 cudaMemcpy(in2_device, in2, size, cudaMemcpyHostToDevice);
 int ShMemSize=nThreadsPerBlock*sizeof(int);
 DotProductArrayWithSMGPU<<<nBlocks, nThreadsPerBlock, ShMemSize>>>
(in1_device, in2_device, PartialSums_device, stride,N);
 cudaThreadSynchronize();
 
 cudaMemcpy(PartialSums_host,PartialSums_device, nBlocks*sizeof(int), cudaMemcpyDeviceToHost);
 for(int i=0;i<nBlocks;i++)
	 *res+=PartialSums_host[i];
	
 free(PartialSums_host);
 cudaFree(in1_device);
 cudaFree(in2_device);
 cudaFree(PartialSums_device);
}

Dot product of two vettors with SM

__global__ void DotProductArrayWithSMGPU (int *in1, int *in2, int *PartialSums, int stride, int N)
{
	 extern __shared__ int shMem[];
	 int i;
	 int PartialSumsIndex=blockIdx.x;
	 int shMemIndex=threadIdx.x;
	 int inIndex = (threadIdx.x+blockIdx.x*blockDim.x)*stride;
	 shMem[shMemIndex]
	 
	 if(inIndex<N){
		for(i=0;i<stride;i++)
			 if(inIndex+i<N)
				 shMem[shMemIndex]+=in1[inIndex+i]*in2[inIndex+i];
			 else
				 break;
		__syncthreads();
		 int value=0;
		 if(threadIdx.x==0){
			 for(i=0;i<blockDim.x;i++)
				 value+=shMem[i];
			 PartialSums[PartialSumsIndex]=value;
		 }
	}
}

Note: Check the computation of If mod(N,L*M)!= 0, appropriate adjustments must be made.

Parallel kernel and the serial version

void RowScaling(int*matrix, int*vector,
int*ris,int rows, int columns)
{
	int i,j;
	for(i=0;i<rows;i++)
	for(j=0;j<columns;j++)
	ris[i*columns+j]= matrix[i*columns+j]*vector[j];
}
 
__global__ void RowScalingGPU(int*mat,int*vec,
int*ris,int M,int N, int pitch)
{
	int indexX=threadIdx.x + blockDim.x*blockIdx.x;
	int indexY=threadIdx.y + blockDim.y*blockIdx.y;
	if(indexX <M && indexY<N)
		ris[indexX*pitch+indexY] = mat[indexX*pitch+indexY] * vec[indexY];
}
 

Exercise

Implement the followings:

• Outer product of two vectors
• Stencil of a vector using SM
• Dot product of two vectors using SM
• Element-wise matrix-vector multiplication using SM

Implement all the CUDA codes with timing and error checks.

---