Vector Addition and CUDA Memory Management

GPUs excel at data-parallel computations, where identical arithmetic operations are applied across vast datasets. To understand GPU memory spaces and kernel execution, we examine vector addition:

$\mathbf{z} = \mathbf{x} + \mathbf{y} \implies z[i] = x[i] + y[i] \quad (0 \le i < n)$

Because GPUs contain significantly more 32-bit floating-point units (FP32) than 64-bit units (FP64), high-performance GPU codes standardly use single-precision float arrays.

---

6.7 Global Thread Indexing

In a 1D grid of 1D blocks, each thread calculates its unique global index using the block dimension, block index, and thread index:

$\text{rank} = \text{blockDim.x} \times \text{blockIdx.x} + \text{threadIdx.x}$

Table 6.4: Global Thread Ranks in a Grid (4 Blocks, 5 Threads per Block)

Block (blockIdx.x)	Thread 0	Thread 1	Thread 2	Thread 3	Thread 4
Block 0	$0$	$1$	$2$	$3$	$4$
Block 1	$5$	$6$	$7$	$8$	$9$
Block 2	$10$	$11$	$12$	$13$	$14$
Block 3	$15$	$16$	$17$	$18$	$19$

---

Handling Arbitrary Problem Sizes: The Boundary Check

In practical applications, the vector size $n$ is rarely an exact multiple of the block size (e.g., $n = 997$). To process all elements, the grid must launch $\lceil n / \text{th\_per\_blk} \rceil$ blocks:

$\text{blk\_ct} = \dfrac{n + \text{th\_per\_blk} - 1}{\text{th\_per\_blk}}$

This launches more total threads than array elements. Every kernel must guard memory accesses with a boundary check:

int my_elt = blockDim.x * blockIdx.x + threadIdx.x;
if (my_elt < n) {
    z[my_elt] = x[my_elt] + y[my_elt];
}

Threads with my_elt >= n simply idle and safely skip memory operations.

---

6.8 Unified Memory: cudaMallocManaged

Introduced in CUDA 6.0, Unified Memory creates a single, managed address space shared between the CPU and GPU. Pointers allocated with cudaMallocManaged can be dereferenced by both host code and device kernels:

__host__ cudaError_t cudaMallocManaged(
    void**   devPtr  /* out: pointer to allocated memory */,
    size_t   size    /* in: size in bytes                */,
    unsigned flags   /* in: access flags (optional)      */
);

The underlying CUDA driver automatically migrates memory pages on-demand across the PCIe/NVLink bus to the processor currently accessing them.

Below is the complete C/CUDA source for Program 6.3:

/* Program 6.3: Vector addition using Unified Memory */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <cuda.h>

/* CUDA Kernel */
__global__ void Vec_add(const float x[], const float y[], float z[], const int n) {
    int my_elt = blockDim.x * blockIdx.x + threadIdx.x;
    if (my_elt < n) {
        z[my_elt] = x[my_elt] + y[my_elt];
    }
}

/* Helper prototypes */
void Allocate_vectors(float** x_p, float** y_p, float** z_p, float** cz_p, int n);
void Free_vectors(float* x, float* y, float* z, float* cz);

int main(int argc, char* argv[]) {
    int n = 1000000;
    int th_per_blk = 256;
    int blk_ct = (n + th_per_blk - 1) / th_per_blk;
    float *x, *y, *z, *cz;

    /* Allocate memory using Unified Memory */
    Allocate_vectors(&x, &y, &z, &cz, n);

    /* Initialize vectors on CPU */
    for (int i = 0; i < n; i++) {
        x[i] = 1.0f;
        y[i] = 2.0f;
    }

    /* Launch kernel on GPU */
    Vec_add<<<blk_ct, th_per_blk>>>(x, y, z, n);
    cudaDeviceSynchronize();

    /* Verify result on CPU */
    printf("z[0] = %f, z[n-1] = %f\n", z[0], z[n - 1]);

    Free_vectors(x, y, z, cz);
    return 0;
}

void Allocate_vectors(float** x_p, float** y_p, float** z_p, float** cz_p, int n) {
    /* Allocated in Unified Memory: accessible by host and device */
    cudaMallocManaged(x_p, n * sizeof(float));
    cudaMallocManaged(y_p, n * sizeof(float));
    cudaMallocManaged(z_p, n * sizeof(float));

    /* Allocated in host RAM: CPU only */
    *cz_p = (float*)malloc(n * sizeof(float));
}

void Free_vectors(float* x, float* y, float* z, float* cz) {
    cudaFree(x);
    cudaFree(y);
    cudaFree(z);
    free(cz);
}

Unified Memory Constraints

• Unified Memory requires a 64-bit operating system and a device with Compute Capability $\ge 3.0$.
• On devices with Compute Capability $< 6.0$, host and device cannot access managed memory concurrently: the host must wait for cudaDeviceSynchronize() before touching managed memory while a kernel runs.

---

6.9 Explicit Memory Management: cudaMalloc and cudaMemcpy

On legacy devices, or in performance-critical code where developers require precise control over when data transfers occur, explicit memory management is used.

Host memory and device memory are maintained using separate pointer sets:

• Host pointers (hx, hy, hz): Points to system RAM allocated with malloc().
• Device pointers (dx, dy, dz): Points to GPU VRAM allocated with cudaMalloc().

/* Allocates linear memory on the GPU device */
__host__ __device__ cudaError_t cudaMalloc(void** devPtr, size_t size);

/* Copies data across the PCIe interconnect */
__host__ cudaError_t cudaMemcpy(
    void*               dest,
    const void*         src,
    size_t              count,
    enum cudaMemcpyKind kind
);

Transfer Directions (cudaMemcpyKind):

• cudaMemcpyHostToDevice: Transfers data from CPU RAM to GPU VRAM.
• cudaMemcpyDeviceToHost: Transfers data from GPU VRAM back to CPU RAM.

flowchart LR
  subgraph Explicit["Figure 6.3: Explicit CUDA Memory Transfer Flow"]
      direction LR
      subgraph HostRAM["Host RAM"]
          HX["hx, hy (Inputs)"]
          HZ["hz (Result)"]
      end
      PCIE["PCIe Bus"]
      subgraph DeviceVRAM["Device VRAM"]
          DX["dx, dy (Device Inputs)"]
          DZ["dz (Device Output)"]
          KERN["Vec_add<<<>>>"]
      end
      
      HX -->|cudaMemcpyHostToDevice| PCIE --> DX
      DX --> KERN --> DZ
      DZ --> PCIE -->|cudaMemcpyDeviceToHost| HZ
  end

Complete Explicit Memory Implementation: Program 6.9

/* Program 6.9: Vector addition with explicit memory copies */
#include <stdio.h>
#include <stdlib.h>
#include <cuda.h>

__global__ void Vec_add(const float x[], const float y[], float z[], const int n) {
    int my_elt = blockDim.x * blockIdx.x + threadIdx.x;
    if (my_elt < n) {
        z[my_elt] = x[my_elt] + y[my_elt];
    }
}

int main(int argc, char* argv[]) {
    int n = 1000000;
    int th_per_blk = 256;
    int blk_ct = (n + th_per_blk - 1) / th_per_blk;

    /* 1. Allocate Host memory */
    float *hx = (float*)malloc(n * sizeof(float));
    float *hy = (float*)malloc(n * sizeof(float));
    float *hz = (float*)malloc(n * sizeof(float));

    for (int i = 0; i < n; i++) {
        hx[i] = 1.0f;
        hy[i] = 2.0f;
    }

    /* 2. Allocate Device memory */
    float *dx, *dy, *dz;
    cudaMalloc(&dx, n * sizeof(float));
    cudaMalloc(&dy, n * sizeof(float));
    cudaMalloc(&dz, n * sizeof(float));

    /* 3. Transfer input data: Host -> Device */
    cudaMemcpy(dx, hx, n * sizeof(float), cudaMemcpyHostToDevice);
    cudaMemcpy(dy, hy, n * sizeof(float), cudaMemcpyHostToDevice);

    /* 4. Launch kernel */
    Vec_add<<<blk_ct, th_per_blk>>>(dx, dy, dz, n);

    /* 5. Transfer result: Device -> Host (cudaMemcpy is synchronous) */
    cudaMemcpy(hz, dz, n * sizeof(float), cudaMemcpyDeviceToHost);

    printf("Result hz[0] = %f\n", hz[0]);

    /* 6. Free memory */
    cudaFree(dx); cudaFree(dy); cudaFree(dz);
    free(hx); free(hy); free(hz);

    return 0;
}

Synchronous Nature of cudaMemcpy

Calls to cudaMemcpy() are synchronous with respect to the host: the function waits for any executing kernels on the device to finish before performing the transfer, and blocks until the memory transfer is complete. Therefore, calling cudaDeviceSynchronize() prior to cudaMemcpy is redundant.

---

6.10 Returning Results from CUDA Kernels

CUDA kernels have a void return type and cannot return scalar values via conventional return statements. Furthermore, passing the address of a host stack variable fails:

/* RUNTIME CRASH: Memory address &sum is invalid on the GPU! */
int sum = -5;
Add<<<1, 1>>>(2, 3, &sum); /* Fails: &sum resides in host CPU stack memory */

Three Patterns to Return Scalar Results:

1. Unified Memory Pointer:

   int* sum_p;
   cudaMallocManaged(&sum_p, sizeof(int));
   Add<<<1, 1>>>(2, 3, sum_p);
   cudaDeviceSynchronize();
   printf("Sum = %d\n", *sum_p);
   cudaFree(sum_p);

2. Explicit Device Memory Allocation & Copy:

   int *h_sum = (int*)malloc(sizeof(int));
   int *d_sum;
   cudaMalloc(&d_sum, sizeof(int));
   Add<<<1, 1>>>(2, 3, d_sum);
   cudaMemcpy(h_sum, d_sum, sizeof(int), cudaMemcpyDeviceToHost);

3. Global Managed Variables (__managed__):

   __managed__ int sum;
   
   __global__ void Add(int x, int y) {
       sum = x + y;
   }
   
   int main(void) {
       Add<<<1, 1>>>(2, 3);
       cudaDeviceSynchronize();
       printf("Sum = %d\n", sum);
       return 0;
   }