📜Kommentierte CUDA-Beispiele

1#include <cstdio>
2#include <cstdlib>
3#include <cuda_runtime.h>
4 
5// Jede CUDA-API-Funktion liefert einen Fehlercode – immer prüfen!
6#define CUDA_CHECK(call)                                                  \
7    do {                                                                  \
8        cudaError_t err_ = (call);                                        \
9        if (err_ != cudaSuccess) {                                        \
10            fprintf(stderr, "CUDA-Fehler: %s (%s:%d)\n",                  \
11                    cudaGetErrorString(err_), __FILE__, __LINE__);        \
12            exit(EXIT_FAILURE);                                           \
13        }                                                                 \
14    } while (0)
15 
16__global__ void vecAdd(const float* a, const float* b, float* c, int n) {
17    int i = blockIdx.x * blockDim.x + threadIdx.x;   // globaler Index
18    if (i < n) c[i] = a[i] + b[i];                   // Wächter gegen Überlauf
19}
20 
21int main() {
22    const int n = 1 << 20;                           // 1 048 576 Elemente
23    const size_t bytes = n * sizeof(float);
24 
25    float* h_a = (float*)malloc(bytes);
26    float* h_b = (float*)malloc(bytes);
27    float* h_c = (float*)malloc(bytes);
28    for (int i = 0; i < n; ++i) { h_a[i] = 1.0f; h_b[i] = 2.0f; }
29 
30    float *d_a, *d_b, *d_c;
31    CUDA_CHECK(cudaMalloc(&d_a, bytes));
32    CUDA_CHECK(cudaMalloc(&d_b, bytes));
33    CUDA_CHECK(cudaMalloc(&d_c, bytes));
34    CUDA_CHECK(cudaMemcpy(d_a, h_a, bytes, cudaMemcpyHostToDevice));
35    CUDA_CHECK(cudaMemcpy(d_b, h_b, bytes, cudaMemcpyHostToDevice));
36 
37    const int block = 256;                           // Vielfaches von 32
38    const int grid = (n + block - 1) / block;        // = 4096 Blöcke
39    vecAdd<<<grid, block>>>(d_a, d_b, d_c, n);
40    CUDA_CHECK(cudaGetLastError());                  // Fehler beim Start?
41    CUDA_CHECK(cudaMemcpy(h_c, d_c, bytes, cudaMemcpyDeviceToHost));  // wartet auf den Kernel
42 
43    for (int i = 0; i < n; ++i)
44        if (h_c[i] != 3.0f) { printf("Fehler bei %d\n", i); return 1; }
45    printf("OK: alle %d Werte = 3\n", n);
46 
47    cudaFree(d_a); cudaFree(d_b); cudaFree(d_c);
48    free(h_a); free(h_b); free(h_c);
49    return 0;
50}

1__inline__ __device__ float warpReduceSum(float v) {
2    // 5 Schritte: 16, 8, 4, 2, 1 – danach hält Lane 0 die Summe des Warps
3    for (int offset = warpSize / 2; offset > 0; offset >>= 1)
4        v += __shfl_down_sync(0xffffffff, v, offset);
5    return v;
6}
7 
8__global__ void reduceSum(const float* __restrict__ in, float* out, int n) {
9    __shared__ float warpSums[32];                 // max. 1024/32 Warps je Block
10    float v = 0.0f;
11    // Grid-Stride: ein festes Grid bearbeitet beliebig große Eingaben
12    for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < n; i += blockDim.x * gridDim.x)
13        v += in[i];
14 
15    int lane = threadIdx.x % warpSize;
16    int warp = threadIdx.x / warpSize;
17    v = warpReduceSum(v);
18    if (lane == 0) warpSums[warp] = v;
19    __syncthreads();
20 
21    if (warp == 0) {
22        v = (lane < blockDim.x / warpSize) ? warpSums[lane] : 0.0f;
23        v = warpReduceSum(v);
24        if (lane == 0) atomicAdd(out, v);          // ein Atomic je Block
25    }
26}
27 
28// Start (out vorher auf 0 setzen!):
29// cudaMemset(d_out, 0, sizeof(float));
30// reduceSum<<<numSMs * 4, 256>>>(d_in, d_out, n);

1#define TILE 16
2 
3__global__ void matmulTiled(const float* __restrict__ A, const float* __restrict__ B,
4                            float* __restrict__ C, int N) {
5    __shared__ float As[TILE][TILE];
6    __shared__ float Bs[TILE][TILE];
7    int tx = threadIdx.x, ty = threadIdx.y;
8    int row = blockIdx.y * TILE + ty;
9    int col = blockIdx.x * TILE + tx;
10    float sum = 0.0f;
11 
12    for (int p = 0; p < (N + TILE - 1) / TILE; ++p) {
13        int aCol = p * TILE + tx;
14        int bRow = p * TILE + ty;
15        // benachbarte Threads (tx) lesen benachbarte Adressen → coalesced
16        As[ty][tx] = (row < N && aCol < N) ? A[row * N + aCol] : 0.0f;
17        Bs[ty][tx] = (bRow < N && col < N) ? B[bRow * N + col] : 0.0f;
18        __syncthreads();
19 
20        #pragma unroll
21        for (int k = 0; k < TILE; ++k)
22            sum += As[ty][k] * Bs[k][tx];   // As: Broadcast, Bs: aufeinanderfolgende Bänke
23        __syncthreads();
24    }
25    if (row < N && col < N) C[row * N + col] = sum;
26}
27 
28// dim3 block(TILE, TILE);                                   // 256 Threads
29// dim3 grid((N + TILE - 1) / TILE, (N + TILE - 1) / TILE);
30// matmulTiled<<<grid, block>>>(d_A, d_B, d_C, N);

1cudaEvent_t start, stop;
2cudaEventCreate(&start);
3cudaEventCreate(&stop);
4 
5cudaEventRecord(start);                    // Marke in den Stream
6vecAdd<<<grid, block>>>(d_a, d_b, d_c, n);
7cudaEventRecord(stop);
8cudaEventSynchronize(stop);                // warten, bis die Marke erreicht ist
9 
10float ms = 0.0f;
11cudaEventElapsedTime(&ms, start, stop);
12double bytes = 3.0 * n * sizeof(float);    // 2 × lesen + 1 × schreiben
13printf("%.3f ms, effektive Bandbreite %.1f GB/s\n", ms, bytes / ms / 1e6);
14 
15cudaEventDestroy(start);
16cudaEventDestroy(stop);

1cudaDeviceProp p;
2cudaGetDeviceProperties(&p, 0);
3printf("%s: CC %d.%d, %d SMs, %zu MB, Warp %d\n",
4       p.name, p.major, p.minor, p.multiProcessorCount, p.totalGlobalMem >> 20, p.warpSize);
5printf("max. %d Threads/Block, %d Threads/SM, %zu KB Shared/Block\n",
6       p.maxThreadsPerBlock, p.maxThreadsPerMultiProcessor, p.sharedMemPerBlock >> 10);
7 
8int minGrid = 0, blockSize = 0;
9cudaOccupancyMaxPotentialBlockSize(&minGrid, &blockSize, vecAdd, 0, 0);
10 
11int blocksPerSM = 0;
12cudaOccupancyMaxActiveBlocksPerMultiprocessor(&blocksPerSM, vecAdd, blockSize, 0);
13float occ = (float)(blocksPerSM * blockSize) / p.maxThreadsPerMultiProcessor;
14printf("Blockgröße %d → %d Blöcke/SM, Occupancy %.0f %%\n", blockSize, blocksPerSM, occ * 100);

1__global__ void saxpy(int n, float a, const float* x, float* y) {
2    for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < n; i += blockDim.x * gridDim.x)
3        y[i] = a * x[i] + y[i];
4}
5 
6int main() {
7    const int n = 1 << 24;
8    float *x, *y;
9    cudaMallocManaged(&x, n * sizeof(float));
10    cudaMallocManaged(&y, n * sizeof(float));
11    for (int i = 0; i < n; ++i) { x[i] = 1.0f; y[i] = 2.0f; }   // CPU schreibt
12 
13    int dev = 0, numSMs = 0;
14    cudaGetDevice(&dev);
15    cudaDeviceGetAttribute(&numSMs, cudaDevAttrMultiProcessorCount, dev);
16    saxpy<<<32 * numSMs, 256>>>(n, 2.0f, x, y);                 // GPU rechnet
17    cudaDeviceSynchronize();                                      // Pflicht vor CPU-Zugriff
18 
19    printf("y[0] = %f (erwartet 4.0)\n", y[0]);
20    cudaFree(x);
21    cudaFree(y);
22    return 0;
23}