Day 14

GPU Kernels Learning Journey

Description

### File: `cmpFHD.cu` **Summary:** Implemented the FHD (Fully-Hybrid Domain) algorithm for non-Cartesian magnetic resonance imaging (MRI) reconstruction in CUDA. The code focuses on optimizing the parallelism structure to handle iterative reconstruction efficiently, aiming to balance computational load while reducing memory footprint. **Learned:** - Gained insights into non-Cartesian MRI imaging techniques and their relevance in modern medical imaging applications. - Developed an understanding of iterative reconstruction methods and how parallelization can significantly improve performance in reconstructing images from non-Cartesian data. - Implemented optimizations to address common challenges in MRI reconstruction, such as memory bandwidth limitations and computational heavy-lifting. ### File: `cmpFHD_real_image.cu` **Summary:** Built upon the previous implementation of the FHD algorithm to include real image reading and processing capabilities. This version takes an actual image, applies the FHD reconstruction algorithm, and outputs the reconstructed image, demonstrating practical applicability of the CUDA code. **Learned:** - Expanded the previous understanding of memory management and kernel optimization by integrating real-world data processing into the workflow. - Familiarized myself with image I/O operations in CUDA, allowing for the handling of real data as input for reconstruction algorithms. ### Reading: - Completed **Chapter 14** of the PMPP book. - Delved into the case study of non-Cartesian magnetic resonance imaging, which provided: - Background on the principles and necessities driving advancements in MRI technology. - A comprehensive look at iterative reconstruction techniques that enhance image quality using statistical estimation methods. - Detailed steps on optimizing the kernel parallelism structure to maximize performance and minimize memory constraints in handling MRI data. - Insights into experimental performance tuning, particularly the advantages of leveraging hardware trigonometry functions to achieve rapid computations. --- ### Day 15 #### File: `flash_attention_backprop.cu` **Summary:** Implemented the backpropagation for Flash Attention in CUDA, continuing from the forward pass developed earlier. The backpropagation step computes the gradients required for training the attention mechanism. However, a small issue arose where some of the gradients are outputting as zero at certain points, which will be addressed and fixed in the coming days. **Learned:** - Explored the process of backpropagation in the context of Flash Attention, including the calculation of gradients for the attention weights and input matrices. - Worked on integrating gradient calculation with memory optimization techniques to maintain efficiency, consistent with the original forward pass. - Identified potential issues related to numerical stability when dealing with gradient flow in CUDA, specifically in the attention layer. --- #### File: `cnn.cu` **Summary:** Developed a Convolutional Neural Network (CNN) implementation in CUDA, including both forward and backward passes with pooling layers. Used the unrolling trick for improved performance in the backward pass, optimizing the matrix operations involved. **Learned:** - Implemented the core components of a CNN in CUDA, including convolutions, activations, pooling layers, and backpropagation. - Utilized the unrolling trick to optimize it, improving the performance of matrix multiplications and gradient calculations. - Gained deeper understanding of the computational requirements for CNN training on GPUs and the importance of efficient memory access patterns and parallelism in deep learning. --- ### Reading: - **Chapter 15:** *Application Case Study—Molecular Visualization and Analysis* - Delved into the background and practical aspects of molecular visualization in parallel computing. - Learned about the importance of thread granularity adjustments and memory coalescing in visualizing large-scale molecular structures using CUDA. - **Chapter 16:** *Application Case Study—Machine Learning* - Focused on Convolutional Neural Networks (ConvNets) and their implementation in CUDA. - Covered key concepts such as basic layers, backpropagation, and the reduction of convolutional layers to matrix multiplication for optimization. - Explored the cuDNN library and its use in accelerating deep learning operations. - **Chapter 17:** *Parallel Programming and Computational Thinking* - Studied the core principles of parallel computing, including problem decomposition, algorithm selection, and computational thinking. - Focused on strategies for optimizing memory locality and shared memory usage in parallel applications. --- ### Day 16 #### Code: `NaiveBayes.cu`, `NaiveBayesKernel.cuh`, `NaiveBayesTrain.cuh`, `NaiveBayesTrain.cpp`, and `main.cpp` **Summary:** Implemented a CUDA-accelerated Naive Bayes classifier, focusing on the training and inference stages. Leveraging shared memory to maximize computational efficiency, the implementation is structured to divide work among threads for parallelized data processing of feature probabilities. **Components Developed:** 1. **`NaiveBayes.cu`**: - This file contains the CUDA kernel responsible for calculating feature likelihoods and class probabilities in parallel. Shared memory was used where possible to minimize global memory access penalties. - Optimized kernel launches to balance between grid and block dimensions for datasets with high dimensionality. 2. **`NaiveBayesKernel.cuh`**: - Header file declaring the kernel functions, ensuring modularity and separation of concerns in code structure. 3. **`NaiveBayesTrain.cuh`**: - Declared the host-side training function, encapsulating the logic to copy data to the GPU, launch CUDA kernels, and retrieve results. 4. **`NaiveBayesTrain.cpp`**: - Implemented the host-side training process, providing pre-processing for input data and managing memory transfers between CPU and GPU. 5. **`main.cpp`**: - Entry point of the program, performing tasks like loading data, splitting datasets for training and testing, and evaluating model performance after training. --- #### Blog Update - Updated My blog with an important information about using NVCC in Colab. - Link: [Learning CUDA with a Weak GPU (or No GPU at All)](https://hamdi.bearblog.dev/learning-cuda-with-a-weak-gpu-or-no-gpu-at-all-yes-you-can/) --- ### Day 17 #### Code: `vec_cublas.cu` **Summary:** Today, I implemented vector addition using the cuBLAS library in CUDA. By leveraging the optimized linear algebra routines provided by cuBLAS, this implementation achieves highly efficient computation of the vector addition operation `C = A + B` for two input vectors `A` and `B`. The addition was performed using the `cublasSaxpy` function, which computes scaled vector addition. **Key Concepts Implemented:** - Used the `cublasSaxpy` function to perform the vector addition in the format `C = alpha * A + B` where `alpha` is a scaling factor. In this case, `alpha` was set to `1.0` to achieve the simple addition of `A` and `B`. - Managed the cuBLAS library handle for the operation. --- **Learned:** - **cuBLAS Basics:** - Gained an introduction to the cuBLAS library and its capabilities for high-performance linear algebra operations. - Learned how to use cuBLAS functions like `cublasSaxpy` for vector addition and understood its parameters. - **cuBLAS Handle Management:** - Understood how to create and destroy a cuBLAS handle using `cublasCreate` and `cublasDestroy`. This is critical for managing state across cuBLAS calls. - **Functionality of `cublasSaxpy`:** - Reviewed the underlying algorithm and implementation of the AXPY operation, which computes `y = a*x + y` for real vectors `x` and `y` and scalar `a`. ---

Code Files

cmpFHD.cu

#include <cuda_runtime.h>
#include <math.h>
#include <stdio.h>
#include <random>
#include <iostream>

#define FHD_THREADS_PER_BLOCK 256
#define PI 3.14159265358979323846
#define CHUNK_SIZE 256

__constant__ float kx_c[CHUNK_SIZE], ky_c[CHUNK_SIZE], kz_c[CHUNK_SIZE];

__global__ void cmpFHd(float* rPhi, float* iPhi, float* phiMag,
                       float* x, float* y, float* z,
                       float* rMu, float* iMu, int M) {
    int n = blockIdx.x * FHD_THREADS_PER_BLOCK + threadIdx.x;

    float xn_r = x[n]; 
    float yn_r = y[n]; 
    float zn_r = z[n];

    float rFhDn_r = rPhi[n]; 
    float iFhDn_r = iPhi[n];

    for (int m = 0; m < M; m++) {
        float expFhD = 2 * PI * (kx_c[m] * xn_r + ky_c[m] * yn_r + kz_c[m] * zn_r);
        
        float cArg = __cosf(expFhD);
        float sArg = __sinf(expFhD);

        rFhDn_r += rMu[m] * cArg - iMu[m] * sArg;
        iFhDn_r += iMu[m] * cArg + rMu[m] * sArg;
    }

    rPhi[n] = rFhDn_r;
    iPhi[n] = iFhDn_r;
    phiMag[n] = sqrtf(rFhDn_r * rFhDn_r + iFhDn_r * iFhDn_r);
}

int main() {
    int N = 1024; // Define problem size
    int M = 1024; // Number of samples

    std::cout << "Starting program..." << std::endl;

    float *x, *y, *z, *rMu, *iMu, *rPhi, *iPhi, *phiMag;
    
    // Allocate memory and check for errors
    cudaError_t cudaStatus;
    // cudaMallocManaged takes care of memory allocation on both CPU and GPU
    cudaStatus = cudaMallocManaged(&x, N * sizeof(float));
    if (cudaStatus != cudaSuccess) {
        std::cout << "cudaMallocManaged failed for x: " << cudaGetErrorString(cudaStatus) << std::endl;
        return 1;
    }

    cudaStatus = cudaMallocManaged(&y, N * sizeof(float));
    cudaStatus = cudaMallocManaged(&z, N * sizeof(float));
    cudaStatus = cudaMallocManaged(&rMu, M * sizeof(float));
    cudaStatus = cudaMallocManaged(&iMu, M * sizeof(float));
    cudaStatus = cudaMallocManaged(&rPhi, N * sizeof(float));
    cudaStatus = cudaMallocManaged(&iPhi, N * sizeof(float));
    cudaStatus = cudaMallocManaged(&phiMag, N * sizeof(float));

    std::cout << "Memory allocated successfully" << std::endl;

    // Initialize random number generator
    std::random_device rd;
    std::mt19937 gen(rd());
    std::uniform_real_distribution<float> dis(-1.0f, 1.0f);

    // Initialize input arrays with random values
    for (int i = 0; i < N; i++) {
        x[i] = dis(gen);
        y[i] = dis(gen);
        z[i] = dis(gen);
        rPhi[i] = 0.0f;
        iPhi[i] = 0.0f;
        phiMag[i] = 0.0f;
    }

    // Initialize rMu and iMu
    for (int i = 0; i < M; i++) {
        rMu[i] = dis(gen);
        iMu[i] = dis(gen);
    }

    std::cout << "Data initialized successfully" << std::endl;

    // Print some initial values
    std::cout << "\nInitial values:" << std::endl;
    for (int i = 0; i < 5; i++) {
        std::cout << "x[" << i << "] = " << x[i] << ", y[" << i << "] = " << y[i] << ", z[" << i << "] = " << z[i] << std::endl;
    }

    // Process data in chunks
    for (int i = 0; i < M / CHUNK_SIZE; i++) {
        std::cout << "\nProcessing chunk " << i + 1 << " of " << M / CHUNK_SIZE << std::endl;

        cudaStatus = cudaMemcpyToSymbol(kx_c, &x[i * CHUNK_SIZE], CHUNK_SIZE * sizeof(float));
        if (cudaStatus != cudaSuccess) {
            std::cout << "cudaMemcpyToSymbol failed: " << cudaGetErrorString(cudaStatus) << std::endl;
            return 1;
        }

        cudaStatus = cudaMemcpyToSymbol(ky_c, &y[i * CHUNK_SIZE], CHUNK_SIZE * sizeof(float));
        cudaStatus = cudaMemcpyToSymbol(kz_c, &z[i * CHUNK_SIZE], CHUNK_SIZE * sizeof(float));

        // Launch kernel
        cmpFHd<<<N / FHD_THREADS_PER_BLOCK, FHD_THREADS_PER_BLOCK>>>(
            rPhi, iPhi, phiMag, x, y, z, rMu, iMu, CHUNK_SIZE);
        
        cudaStatus = cudaGetLastError();
        if (cudaStatus != cudaSuccess) {
            std::cout << "Kernel launch failed: " << cudaGetErrorString(cudaStatus) << std::endl;
            return 1;
        }

        cudaStatus = cudaDeviceSynchronize();
        if (cudaStatus != cudaSuccess) {
            std::cout << "cudaDeviceSynchronize failed: " << cudaGetErrorString(cudaStatus) << std::endl;
            return 1;
        }
    }

    std::cout << "\nComputation completed. Results:" << std::endl;
    for (int i = 0; i < 5; i++) {
        std::cout << "rPhi[" << i << "] = " << rPhi[i] 
                 << ", iPhi[" << i << "] = " << iPhi[i] 
                 << ", phiMag[" << i << "] = " << phiMag[i] << std::endl;
    }

    // Free memory
    cudaFree(x);
    cudaFree(y);
    cudaFree(z);
    cudaFree(rMu);
    cudaFree(iMu);
    cudaFree(rPhi);
    cudaFree(iPhi);
    cudaFree(phiMag);

    std::cout << "\nProgram completed successfully" << std::endl;
    return 0;
}

cmpFHd_real_image.cu

/*
nvcc cmpFHd_real_image.cu -o cmpFHd_real_image `pkg-config --cflags --libs opencv4`

and don't forget to put the image file in the same directory.
*/



#include <cuda_runtime.h>
#include <opencv2/opencv.hpp>
#include <iostream>

#define FHD_THREADS_PER_BLOCK 256
#define PI 3.14159265358979323846
#define CHUNK_SIZE 256

using namespace cv;
using namespace std;

__constant__ float kx_c[CHUNK_SIZE], ky_c[CHUNK_SIZE], kz_c[CHUNK_SIZE];

__global__ void cmpFHd(float* rPhi, float* iPhi, float* phiMag,
                       float* x, float* y, float* z,
                       float* rMu, float* iMu, int M) {
    int n = blockIdx.x * FHD_THREADS_PER_BLOCK + threadIdx.x;
    
    float xn_r = x[n]; 
    float yn_r = y[n]; 
    float zn_r = z[n];

    float rFhDn_r = rPhi[n]; 
    float iFhDn_r = iPhi[n];

    for (int m = 0; m < M; m++) {
        float expFhD = 2 * PI * (kx_c[m] * xn_r + ky_c[m] * yn_r + kz_c[m] * zn_r);
        
        float cArg = __cosf(expFhD);
        float sArg = __sinf(expFhD);

        rFhDn_r += rMu[m] * cArg - iMu[m] * sArg;
        iFhDn_r += iMu[m] * cArg + rMu[m] * sArg;
    }

    rPhi[n] = rFhDn_r;
    iPhi[n] = iFhDn_r;
}

int main() {
    // Load an image using OpenCV
    Mat image = imread("lena_gray.png", IMREAD_GRAYSCALE);
    if (image.empty()) {
        cerr << "Error: Could not open the image!" << endl;
        return -1;
    }

    // Normalize image to range [0,1]
    image.convertTo(image, CV_32F, 1.0 / 255);

    int N = image.rows * image.cols; // Number of pixels
    int M = 256; // Number of frequency components
  
    float *x, *y, *z, *rMu, *iMu, *rPhi, *iPhi, *phiMag;
    
    // Allocate CUDA memory
    cudaMallocManaged(&x, N * sizeof(float));
    cudaMallocManaged(&y, N * sizeof(float));
    cudaMallocManaged(&z, N * sizeof(float));
    cudaMallocManaged(&rMu, M * sizeof(float));
    cudaMallocManaged(&iMu, M * sizeof(float));
    cudaMallocManaged(&rPhi, N * sizeof(float));
    cudaMallocManaged(&iPhi, N * sizeof(float));
    cudaMallocManaged(&phiMag, N * sizeof(float));

    // Initialize x, y coordinates from image pixels
    for (int i = 0; i < image.rows; i++) {
        for (int j = 0; j < image.cols; j++) {
            int idx = i * image.cols + j;
            x[idx] = (float)j / image.cols;  // Normalize to [0,1]
            y[idx] = (float)i / image.rows;  // Normalize to [0,1]
            z[idx] = image.at<float>(i, j);  // Use intensity as "z"
            rPhi[idx] = z[idx];              // Initial real part
            iPhi[idx] = 0.0f;                // Initial imaginary part
        }
    }

    // Randomly initialize rMu and iMu
    for (int i = 0; i < M; i++) {
        rMu[i] = static_cast<float>(rand()) / RAND_MAX;
        iMu[i] = static_cast<float>(rand()) / RAND_MAX;
    }

    // Copy chunks of kx, ky, kz to constant memory
    for (int i = 0; i < M / CHUNK_SIZE; i++) {
        cudaMemcpyToSymbol(kx_c, &x[i * CHUNK_SIZE], CHUNK_SIZE * sizeof(float));
        cudaMemcpyToSymbol(ky_c, &y[i * CHUNK_SIZE], CHUNK_SIZE * sizeof(float));
        cudaMemcpyToSymbol(kz_c, &z[i * CHUNK_SIZE], CHUNK_SIZE * sizeof(float));

        // Launch CUDA kernel
        cmpFHd<<<N / FHD_THREADS_PER_BLOCK, FHD_THREADS_PER_BLOCK>>>(rPhi, iPhi, phiMag, x, y, z, rMu, iMu, CHUNK_SIZE);
        cudaDeviceSynchronize();
    }

    // Convert results back to image format
    Mat outputImage(image.rows, image.cols, CV_32F);
    for (int i = 0; i < image.rows; i++) {
        for (int j = 0; j < image.cols; j++) {
            int idx = i * image.cols + j;
            outputImage.at<float>(i, j) = sqrt(rPhi[idx] * rPhi[idx] + iPhi[idx] * iPhi[idx]);
        }
    }

    // Normalize and save output image
    normalize(outputImage, outputImage, 0, 255, NORM_MINMAX);
    outputImage.convertTo(outputImage, CV_8U);
    imwrite("output.jpg", outputImage);

    // Free memory
    cudaFree(x);
    cudaFree(y);
    cudaFree(z);
    cudaFree(rMu);
    cudaFree(iMu);
    cudaFree(rPhi);
    cudaFree(iPhi);
    cudaFree(phiMag);

    cout << "Processed image saved as output.jpg" << endl;
    return 0;
}