Day 7

GPU Kernels Learning Journey

Description

### File: `one_d_convolution.cu` **Summary:** Implemented a simple 1D convolution algorithm using CUDA. This involved sliding a kernel (or filter) over an input array and computing the weighted sum of elements. Each thread was assigned to compute the convolution at a specific position in the output array. **Learned:** - Basics of 1D convolution in parallel, including mapping threads to positions in the output array. - How to handle boundary conditions (halo cells) when the kernel partially overlaps the input array bounds. - Importance of memory layout and contiguous access for kernel weights and input arrays to maximize performance. --- ### File: `one_d_convolution_with_tiling.cu` **Summary:** Implemented an optimized version of the 1D convolution algorithm using tiling and shared memory. Divided the input array into tiles and loaded data into shared memory, minimizing global memory accesses for better performance. Used halo cells to handle edge cases where kernel overlap extended into neighboring tiles. **Learned:** - Tiling in CUDA: Dividing input data into manageable chunks and leveraging shared memory to reduce global memory latency. - Use of **halo cells** to ensure correctness at tile boundaries during convolution. - How to balance computation and memory usage in tiled algorithms to improve performance. - Proper synchronization of threads within a block (using `__syncthreads()`) to ensure data consistency in shared memory. --- ### File: `2d_convolution_with_tiling.cu` **Summary:** Implemented a 2D convolution algorithm with tiling optimization using CUDA. Divided the input matrix into tiles and leveraged shared memory to minimize global memory accesses, ensuring efficient computation of the convolution kernel across the matrix. Handled boundary conditions using halo cells to process edges and corners correctly. **Learned:** - Extended tiling techniques from 1D to 2D data structures for efficient parallel computation. - Optimized global memory access by using shared memory for each tile. - Synchronization of threads for consistent shared memory usage within a block (`__syncthreads()` for proper execution order). - Efficient handling of edge cases and boundary cells in 2D convolution. --- ### Reading: - Read **Chapter 7** of the PMPP book. - Learned about parallel patterns for convolution, including basic algorithms, memory optimizations with constant and shared memory, and tiling techniques with halo cells for 1D and 2D convolution.

Code Files

one_d_convolution.cu

#include <iostream>
#include <cuda_runtime.h>
#define Mask_width 5  
__constant__ float M[Mask_width];
__global__ void oned_convolution_kernel(const float* A,float* C,int n ){
//without tiling 
int threadId=threadIdx.x;
int i=blockDim.x*blockIdx.x+threadId;

if (i<n){
float result=0.0f;
for (int k=-1*Mask_width/2;k<Mask_width/2+1;k++) {
  printf("%.i",k);
  if (i+k>=0 && i+k<n) {
    
  result+=A[i+k]*M[k+Mask_width/2];
  
}}
C[i]=result;

}
}





// Host function to check for CUDA errors
void checkCudaError(const char* message) {
    cudaError_t error = cudaGetLastError();
    if (error != cudaSuccess) {
        std::cerr << message << " - CUDA Error: " << cudaGetErrorString(error) << std::endl;
        exit(EXIT_FAILURE);
    }
}


int main(){
 
  int n=10;
  float A[n],C[n];
  float d_M[Mask_width];
  
   for (int i=0; i<Mask_width;i++){
    d_M[i]=i;

  }
  for (int i=0; i<n;i++){
    A[i]=i;

  }

  float *d_a,*d_c;
  cudaMalloc(&d_a,n*sizeof(float));
  cudaMalloc(&d_c,n*sizeof(float));
  cudaMemcpy(d_a,A,n*sizeof(float),cudaMemcpyHostToDevice);
  checkCudaError("Failed to copy input data to device");
  cudaMemcpyToSymbol(M,d_M,Mask_width*sizeof(float));
  checkCudaError("Failed to copy mask data to device");
  dim3 dimBlock(32);
  dim3 dimGrid((n + dimBlock.x - 1) / dimBlock.x);
  oned_convolution_kernel<<<dimGrid, dimBlock>>>(d_a,d_c,n);
  checkCudaError("Failed to execute the kernel");
  cudaDeviceSynchronize();
  cudaMemcpy(C,d_c,n*sizeof(float),cudaMemcpyDeviceToHost);
checkCudaError("Failed to copy output data to host");
  cudaFree(d_a);
  cudaFree(d_c);
  

  //printing the results
  printf("A:\n");
  for (int i=0; i<n;i++){
    printf("%.2f ", A[i]);

  }
  printf("\n");
   printf("\nd_m:\n");
    for (int i = 0; i < Mask_width; i++) {

            printf("%.2f ", d_M[i]);

    }
  printf("\n");
  printf("\nC:\n");
    for (int i = 0; i < n; i++) {

            printf("%.2f ", C[i]);

    }
  printf("\n");
}

2d_convolution_with_tiling.cu

#include <stdio.h>
#include <iostream>

// I'm assuming that the mask and the matrix to be square for simplicity
#define Mask_width 5
#define shared_size (32 + Mask_width - 1)
__constant__ float M[Mask_width][Mask_width];

__global__ void twod_convolution_kernel(const float* A, float* C, int n) {
    int threadx = threadIdx.x;
    int thready = threadIdx.y;
    int i = blockDim.x * blockIdx.x + threadx;
    int j = blockDim.y * blockIdx.y + thready;
    
    __shared__ float S_A[shared_size][shared_size];

    // Load main data
    if ((i < n) && (j < n)) {
        S_A[threadx + Mask_width/2][thready + Mask_width/2] = A[i*n+j];
    }

    // Load left halo
    if (threadx < Mask_width/2) {
        int left_idx = blockIdx.x * blockDim.x - (Mask_width/2) + threadx;
        if (left_idx >= 0 && j < n) {
            S_A[threadx][thready + Mask_width/2] = A[left_idx*n+j];
        }
        else {
            S_A[threadx][thready + Mask_width/2] = 0.0f;
        }
    }

    // Load right halo
    if (threadx < Mask_width/2) {
        int right_idx = blockIdx.x * blockDim.x + blockDim.x + threadx;
        if (right_idx < n && j < n) {
            S_A[threadx + blockDim.x + Mask_width/2][thready + Mask_width/2] = A[right_idx*n+j];
        }
        else {
            S_A[threadx + blockDim.x + Mask_width/2][thready + Mask_width/2] = 0.0f;
        }
    }

    // Load top halo
    if (thready < Mask_width/2) {
        int top_idy = j - (Mask_width/2) + thready;
        if (top_idy >= 0 && i < n) {
            S_A[threadx + Mask_width/2][thready] = A[i*n+top_idy];
        }
        else {
            S_A[threadx + Mask_width/2][thready] = 0.0f;
        }
    }

    // Load bottom halo
    if (thready < Mask_width/2) {
        int bottom_idy = j + blockDim.y + thready;
        if (bottom_idy < n && i < n) {
            S_A[threadx + Mask_width/2][thready + blockDim.y + Mask_width/2] = A[i*n+bottom_idy];
        }
        else {
            S_A[threadx + Mask_width/2][thready + blockDim.y + Mask_width/2] = 0.0f;
        }
    }

    __syncthreads();

    if ((i < n) && (j < n)) {
        float result = 0.0f;
        for (int k = 0; k < Mask_width; k++) {
            for (int x = 0; x < Mask_width; x++) {
                result += S_A[threadx + k][thready + x] * M[k][x];
            }
        }
        C[i*n+j] = result;
    }
}

void checkCudaError(const char* message) {
    cudaError_t error = cudaGetLastError();
    if (error != cudaSuccess) {
        fprintf(stderr, "%s - CUDA Error: %s\n", message, cudaGetErrorString(error));
        exit(EXIT_FAILURE);
    }
}

int main() {
    int n = 10;
    float *h_A = (float*)malloc(n * n * sizeof(float));
    float *h_C = (float*)malloc(n * n * sizeof(float));
    float d_M[Mask_width][Mask_width];

    for (int i = 0; i < Mask_width; i++) {
        for (int j = 0; j < Mask_width; j++) {
            d_M[i][j] = 5;
        }
    }

    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            h_A[i*n + j] = 3;
        }
    }

    float *d_a, *d_c;
    cudaMalloc(&d_a, n*n*sizeof(float));
    cudaMalloc(&d_c, n*n*sizeof(float));
    cudaMemcpy(d_a, h_A, n*n*sizeof(float), cudaMemcpyHostToDevice);
    checkCudaError("Failed to copy input data to device");
    cudaMemcpyToSymbol(M, d_M, Mask_width*Mask_width*sizeof(float));
    checkCudaError("Failed to copy mask data to device");

    dim3 dimBlock(32, 32);
    dim3 dimGrid((n + dimBlock.x - 1) / dimBlock.x, (n + dimBlock.y - 1) / dimBlock.y);
    twod_convolution_kernel<<<dimGrid, dimBlock>>>(d_a, d_c, n);
    checkCudaError("Failed to execute the kernel");
    
    cudaDeviceSynchronize();
    cudaMemcpy(h_C, d_c, n*n*sizeof(float), cudaMemcpyDeviceToHost);
    checkCudaError("Failed to copy output data to host");

    // Print results
    printf("Results:\n");
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            printf("%.2f ", h_C[i*n + j]);
        }
        printf("\n");
    }

    // Clean up
    cudaFree(d_a);
    cudaFree(d_c);
    free(h_A);
    free(h_C);

    return 0;
}

one_d_convolution_with_tiling.cu

#include <iostream>
#include <cuda_runtime.h>
#define Mask_width 5  
__constant__ float M[Mask_width];

__global__ void oned_convolution_tiling_kernel(const float* A, float* C, int n) {
    int threadId = threadIdx.x;
    int i = blockDim.x * blockIdx.x + threadId;
    
    __shared__ float S_A[32 + Mask_width - 1];
    
    // Load main data
    if (i < n) {
        S_A[threadId + Mask_width/2] = A[i];
    }
    
    // Load left halo
    if (threadId < Mask_width/2) {
        int left_idx = blockIdx.x * blockDim.x - (Mask_width/2) + threadId;
        if (left_idx >= 0) {
            S_A[threadId] = A[left_idx];
        }
        else {
            S_A[threadId] = 0.0f;
        }
    }
    
    // Load right halo
    if (threadId < Mask_width/2) {
        int right_idx = blockIdx.x * blockDim.x + blockDim.x + threadId;
        if (right_idx < n) {
            S_A[threadId + blockDim.x + Mask_width/2] = A[right_idx];
        }
        else {
            S_A[threadId + blockDim.x + Mask_width/2] = 0.0f;
        }
    }
    
    __syncthreads();
    
    if (i < n) {
        float result = 0.0f;
        for (int k = 0; k < Mask_width; k++) {
            int idx = threadId + k;
            if ((i + k - Mask_width/2) >= 0 && (i + k - Mask_width/2) < n) {
                result += S_A[idx] * M[k];
            }
        }
        C[i] = result;
    }
}




// Host function to check for CUDA errors
void checkCudaError(const char* message) {
    cudaError_t error = cudaGetLastError();
    if (error != cudaSuccess) {
        std::cerr << message << " - CUDA Error: " << cudaGetErrorString(error) << std::endl;
        exit(EXIT_FAILURE);
    }
}


int main(){
 
  int n=10;
  float A[n],C[n];
  float d_M[Mask_width];
  
   for (int i=0; i<Mask_width;i++){
    d_M[i]=i;

  }
  for (int i=0; i<n;i++){
    A[i]=i;

  }

  float *d_a,*d_c;
  cudaMalloc(&d_a,n*sizeof(float));
  cudaMalloc(&d_c,n*sizeof(float));
  cudaMemcpy(d_a,A,n*sizeof(float),cudaMemcpyHostToDevice);
  checkCudaError("Failed to copy input data to device");
  cudaMemcpyToSymbol(M,d_M,Mask_width*sizeof(float));
  checkCudaError("Failed to copy mask data to device");
  dim3 dimBlock(32);
  dim3 dimGrid((n + dimBlock.x - 1) / dimBlock.x);
  oned_convolution_tiling_kernel<<<dimGrid, dimBlock>>>(d_a,d_c,n);
  checkCudaError("Failed to execute the kernel");
  cudaDeviceSynchronize();
  cudaMemcpy(C,d_c,n*sizeof(float),cudaMemcpyDeviceToHost);
checkCudaError("Failed to copy output data to host");
  cudaFree(d_a);
  cudaFree(d_c);
  

  //printing the results
  printf("A:\n");
  for (int i=0; i<n;i++){
    printf("%.2f ", A[i]);

  }
  printf("\n");
   printf("\nd_m:\n");
    for (int i = 0; i < Mask_width; i++) {

            printf("%.2f ", d_M[i]);

    }
  printf("\n");
  printf("\nC:\n");
    for (int i = 0; i < n; i++) {

            printf("%.2f ", C[i]);

    }
  printf("\n");
}