CUDA Hands-on 완벽 가이드: GPU 컴퓨팅의 모든 것

GPU (GPC - Graphics Processing Cluster)
├── SM (Streaming Multiprocessor) × N개
│   ├── CUDA Core (FP32/INT32) × 128개 (Hopper 기준)
│   ├── Tensor Core × 4개 (4세대, Hopper 기준)
│   ├── RT Core (Ray Tracing) × 1개
│   ├── Warp Scheduler × 4개
│   ├── Register File (256 KB)
│   ├── L1 Cache / Shared Memory (공유, 최대 228 KB)
│   └── SFU (Special Function Unit)
├── L2 Cache (공유)
├── Memory Controller
└── HBM (High Bandwidth Memory)

속도:  Register > Shared/L1 > L2 > Global (HBM)
크기:  Global (HBM) > L2 > Shared/L1 > Register

Grid (커널 실행 단위)
├── Block (0,0)          Block (1,0)          Block (2,0)
│   ├── Warp 0           ├── Warp 0           ├── Warp 0
│   │   ├── Thread 0     │   ├── Thread 0     │   ├── Thread 0
│   │   ├── Thread 1     │   ├── Thread 1     │   ├── Thread 1
│   │   ├── ...          │   ├── ...          │   ├── ...
│   │   └── Thread 31    │   └── Thread 31    │   └── Thread 31
│   ├── Warp 1           ├── Warp 1           ├── Warp 1
│   └── ...              └── ...              └── ...
├── Block (0,1)          Block (1,1)          ...
└── ...

# 1. NVIDIA 드라이버 확인
nvidia-smi

# 2. CUDA Keyring 추가 (Ubuntu 22.04 예시)
wget https://developer.download.nvidia.com/compute/cuda/repos/ubuntu2204/x86_64/cuda-keyring_1.1-1_all.deb
sudo dpkg -i cuda-keyring_1.1-1_all.deb
sudo apt-get update

# 3. CUDA Toolkit 설치
sudo apt-get install cuda-toolkit-13-1

# 4. 환경 변수 설정
export PATH=/usr/local/cuda-13.1/bin:$PATH
export LD_LIBRARY_PATH=/usr/local/cuda-13.1/lib64:$LD_LIBRARY_PATH

# 5. 설치 확인
nvcc --version

# 1. NVIDIA 공식 사이트에서 CUDA Toolkit 다운로드
# https://developer.nvidia.com/cuda-downloads

# 2. 설치 후 환경 변수 확인
echo %CUDA_PATH%
# 일반적으로 C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v13.1

# 3. 확인
nvcc --version
nvidia-smi

# 기본 정보 출력
nvidia-smi

# 출력 예시:
# +-----------------------------------------------------------------------------------------+
# | NVIDIA-SMI 560.35.03    Driver Version: 560.35.03    CUDA Version: 13.1                 |
# |-----------------------------------------+------------------------+----------------------|
# | GPU  Name                 Persistence-M | Bus-Id          Disp.A | Volatile Uncorr. ECC |
# | Fan  Temp   Perf          Pwr:Usage/Cap |           Memory-Usage | GPU-Util  Compute M. |
# |=========================================+========================+======================|
# |   0  NVIDIA H100 80GB HBM3          On | 00000000:3B:00.0   Off |                    0 |
# | N/A   32C    P0              72W / 700W |    1234MiB / 81559MiB  |      0%      Default |
# +-----------------------------------------+------------------------+----------------------+

# 지속적 모니터링 (1초 간격)
nvidia-smi -l 1

# dmon - GPU 메트릭 모니터링
nvidia-smi dmon -s pucvmet -d 1

# 특정 GPU 정보
nvidia-smi -i 0 -q

# GPU 프로세스 확인
nvidia-smi pmon -i 0

# GPU 클럭, 메모리 정보
nvidia-smi -q -d CLOCK,MEMORY

# CSV 형식으로 출력 (스크립트용)
nvidia-smi --query-gpu=name,temperature.gpu,utilization.gpu,memory.used,memory.total \
  --format=csv,noheader,nounits

# 기본 컴파일
nvcc -o hello hello.cu

# 특정 아키텍처 타겟
nvcc -arch=sm_90 -o kernel kernel.cu   # Hopper H100
nvcc -arch=sm_80 -o kernel kernel.cu   # Ampere A100

# 디버그 모드
nvcc -g -G -o debug_kernel kernel.cu

# 최적화 레벨
nvcc -O3 -o optimized kernel.cu

# PTX 코드 생성 (중간 표현)
nvcc -ptx kernel.cu

# 여러 아키텍처 동시 지원 (Fat Binary)
nvcc -gencode arch=compute_80,code=sm_80 \
     -gencode arch=compute_90,code=sm_90 \
     -o multi_arch kernel.cu

# 라이브러리 링크
nvcc -lcublas -lcurand -o linked_kernel kernel.cu

# 상세 컴파일 정보
nvcc --resource-usage -o kernel kernel.cu

# 1. 드라이버 지원 CUDA 버전 (Driver API)
nvidia-smi
# 우측 상단에 "CUDA Version: 13.1" 표시
# 이 버전은 드라이버가 지원하는 최대 CUDA 런타임 버전

# 2. CUDA Toolkit 버전 (Runtime API)
nvcc --version
# "release 13.1" 표시
# 실제 설치된 CUDA Toolkit 버전

# 3. 런타임에서 확인
python3 -c "import torch; print(torch.version.cuda)"

# 1. NVIDIA Container Toolkit 설치 (Ubuntu)
curl -fsSL https://nvidia.github.io/libnvidia-container/gpgkey | \
  sudo gpg --dearmor -o /usr/share/keyrings/nvidia-container-toolkit-keyring.gpg

curl -s -L https://nvidia.github.io/libnvidia-container/stable/deb/nvidia-container-toolkit.list | \
  sed 's#deb https://#deb [signed-by=/usr/share/keyrings/nvidia-container-toolkit-keyring.gpg] https://#g' | \
  sudo tee /etc/apt/sources.list.d/nvidia-container-toolkit.list

sudo apt-get update
sudo apt-get install -y nvidia-container-toolkit

# 2. Docker 런타임 설정
sudo nvidia-ctk runtime configure --runtime=docker
sudo systemctl restart docker

# 3. GPU 컨테이너 실행
docker run --rm --gpus all nvidia/cuda:13.1.0-base-ubuntu22.04 nvidia-smi

# 특정 GPU만 사용
docker run --rm --gpus '"device=0,1"' nvidia/cuda:13.1.0-base-ubuntu22.04 nvidia-smi

# docker compose에서 GPU 사용
# docker-compose.yml:

# docker-compose.yml
services:
  gpu-app:
    image: nvidia/cuda:13.1.0-devel-ubuntu22.04
    deploy:
      resources:
        reservations:
          devices:
            - driver: nvidia
              count: all
              capabilities: [gpu]
    command: nvidia-smi

# PyTorch (CUDA 12.4 빌드 예시)
pip install torch torchvision torchaudio --index-url https://download.pytorch.org/whl/cu124

# conda를 사용하는 경우
conda install pytorch torchvision torchaudio pytorch-cuda=12.4 -c pytorch -c nvidia

# TensorFlow (GPU 자동 감지)
pip install tensorflow[and-cuda]

# CUDA 사용 가능 여부 확인 (Python)
python3 -c "
import torch
print(f'PyTorch version: {torch.__version__}')
print(f'CUDA available: {torch.cuda.is_available()}')
print(f'CUDA version: {torch.version.cuda}')
print(f'GPU count: {torch.cuda.device_count()}')
if torch.cuda.is_available():
    print(f'GPU name: {torch.cuda.get_device_name(0)}')
"

┌─────────────────────────────────────────────────────────┐
│  Host (CPU)                                             │
│  ┌─────────────────────────────────────────────────────┐│
│  │ 1. 데이터 준비 (Host Memory)                       ││
│  │ 2. GPU 메모리 할당 (cudaMalloc)                    ││
│  │ 3. 데이터 전송: Host -> Device (cudaMemcpy)        ││
│  │ 4. 커널 실행 (<<<grid, block>>>)                   ││
│  │ 5. 결과 전송: Device -> Host (cudaMemcpy)          ││
│  │ 6. GPU 메모리 해제 (cudaFree)                      ││
│  └─────────────────────────────────────────────────────┘│
│                        ↕ PCIe / NVLink                  │
│  ┌─────────────────────────────────────────────────────┐│
│  │  Device (GPU)                                       ││
│  │  - 커널 함수 병렬 실행                              ││
│  │  - 수천 스레드가 동시 처리                           ││
│  └─────────────────────────────────────────────────────┘│
└─────────────────────────────────────────────────────────┘

// __global__: 커널 함수 - Host에서 호출, Device에서 실행
__global__ void myKernel(float* data, int n) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < n) {
        data[idx] *= 2.0f;
    }
}

// __device__: Device 전용 헬퍼 함수
__device__ float square(float x) {
    return x * x;
}

// __host__ __device__: CPU/GPU 겸용 함수
__host__ __device__ float add(float a, float b) {
    return a + b;
}

// 커널 호출 문법
// kernel<<<gridDim, blockDim, sharedMemSize, stream>>>(args...);
//
// gridDim:       Grid 내 Block 수 (dim3)
// blockDim:      Block 내 Thread 수 (dim3)
// sharedMemSize: 동적 Shared Memory 크기 (바이트, 선택)
// stream:        실행 스트림 (선택)

// 1D 예시: 256개 스레드, 1개 블록
myKernel<<<1, 256>>>(d_data, n);

// 1D 예시: N개 데이터를 256 스레드 블록으로 처리
int threadsPerBlock = 256;
int blocksPerGrid = (n + threadsPerBlock - 1) / threadsPerBlock;
myKernel<<<blocksPerGrid, threadsPerBlock>>>(d_data, n);

// 2D 예시: 이미지 처리
dim3 blockDim(16, 16);       // 16x16 = 256 스레드/블록
dim3 gridDim(
    (width + 15) / 16,
    (height + 15) / 16
);
imageKernel<<<gridDim, blockDim>>>(d_image, width, height);

// 내장 변수 (Built-in Variables)
threadIdx.x, threadIdx.y, threadIdx.z  // Block 내 Thread 인덱스
blockIdx.x,  blockIdx.y,  blockIdx.z   // Grid 내 Block 인덱스
blockDim.x,  blockDim.y,  blockDim.z   // Block의 차원 크기
gridDim.x,   gridDim.y,   gridDim.z    // Grid의 차원 크기

// 글로벌 Thread ID 계산 (1D)
int globalIdx = blockIdx.x * blockDim.x + threadIdx.x;

// 글로벌 Thread ID 계산 (2D)
int col = blockIdx.x * blockDim.x + threadIdx.x;
int row = blockIdx.y * blockDim.y + threadIdx.y;
int globalIdx = row * width + col;

// 전체 스레드 수 (Grid-stride loop용)
int totalThreads = gridDim.x * blockDim.x;

// GPU 메모리 할당
float* d_data;
cudaMalloc((void**)&d_data, n * sizeof(float));

// Host -> Device 복사
cudaMemcpy(d_data, h_data, n * sizeof(float), cudaMemcpyHostToDevice);

// Device -> Host 복사
cudaMemcpy(h_data, d_data, n * sizeof(float), cudaMemcpyDeviceToHost);

// GPU 메모리 해제
cudaFree(d_data);

// GPU 메모리 초기화
cudaMemset(d_data, 0, n * sizeof(float));

// 에러 체크 매크로 (필수!)
#define CUDA_CHECK(call) \
    do { \
        cudaError_t err = call; \
        if (err != cudaSuccess) { \
            fprintf(stderr, "CUDA Error at %s:%d - %s\n", \
                    __FILE__, __LINE__, cudaGetErrorString(err)); \
            exit(EXIT_FAILURE); \
        } \
    } while(0)

// 사용 예
CUDA_CHECK(cudaMalloc((void**)&d_data, n * sizeof(float)));
CUDA_CHECK(cudaMemcpy(d_data, h_data, n * sizeof(float), cudaMemcpyHostToDevice));

// hello_cuda.cu
#include <stdio.h>

__global__ void helloKernel() {
    int tid = blockIdx.x * blockDim.x + threadIdx.x;
    printf("Hello from GPU thread %d (block %d, thread %d)\n",
           tid, blockIdx.x, threadIdx.x);
}

int main() {
    // 2개 블록, 각 4개 스레드 = 총 8개 스레드
    helloKernel<<<2, 4>>>();

    // GPU 작업 완료 대기
    cudaDeviceSynchronize();

    printf("Hello from CPU!\n");
    return 0;
}

# 컴파일 및 실행
nvcc -o hello hello_cuda.cu
./hello

// vector_add.cu
#include <stdio.h>
#include <stdlib.h>
#include <cuda_runtime.h>

#define CUDA_CHECK(call) \
    do { \
        cudaError_t err = call; \
        if (err != cudaSuccess) { \
            fprintf(stderr, "CUDA Error: %s at %s:%d\n", \
                    cudaGetErrorString(err), __FILE__, __LINE__); \
            exit(EXIT_FAILURE); \
        } \
    } while(0)

// GPU 커널: 벡터 덧셈
__global__ void vectorAdd(const float* A, const float* B, float* C, int n) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < n) {
        C[idx] = A[idx] + B[idx];
    }
}

int main() {
    const int N = 1 << 20;  // 1M 원소
    size_t bytes = N * sizeof(float);

    // Host 메모리 할당 및 초기화
    float *h_A = (float*)malloc(bytes);
    float *h_B = (float*)malloc(bytes);
    float *h_C = (float*)malloc(bytes);

    for (int i = 0; i < N; i++) {
        h_A[i] = (float)i;
        h_B[i] = (float)(i * 2);
    }

    // Device 메모리 할당
    float *d_A, *d_B, *d_C;
    CUDA_CHECK(cudaMalloc(&d_A, bytes));
    CUDA_CHECK(cudaMalloc(&d_B, bytes));
    CUDA_CHECK(cudaMalloc(&d_C, bytes));

    // Host -> Device 복사
    CUDA_CHECK(cudaMemcpy(d_A, h_A, bytes, cudaMemcpyHostToDevice));
    CUDA_CHECK(cudaMemcpy(d_B, h_B, bytes, cudaMemcpyHostToDevice));

    // 커널 실행
    int threadsPerBlock = 256;
    int blocksPerGrid = (N + threadsPerBlock - 1) / threadsPerBlock;
    vectorAdd<<<blocksPerGrid, threadsPerBlock>>>(d_A, d_B, d_C, N);

    // 에러 체크
    CUDA_CHECK(cudaGetLastError());
    CUDA_CHECK(cudaDeviceSynchronize());

    // Device -> Host 복사
    CUDA_CHECK(cudaMemcpy(h_C, d_C, bytes, cudaMemcpyDeviceToHost));

    // 결과 검증
    for (int i = 0; i < N; i++) {
        if (fabs(h_C[i] - (h_A[i] + h_B[i])) > 1e-5) {
            fprintf(stderr, "Verification failed at index %d!\n", i);
            exit(EXIT_FAILURE);
        }
    }
    printf("Vector addition of %d elements: SUCCESS\n", N);

    // 메모리 해제
    cudaFree(d_A); cudaFree(d_B); cudaFree(d_C);
    free(h_A); free(h_B); free(h_C);

    return 0;
}

// matmul.cu
#include <stdio.h>
#include <cuda_runtime.h>

#define TILE_SIZE 16

// 기본 행렬 곱셈 커널
__global__ void matMulBasic(const float* A, const float* B, float* C,
                            int M, int N, int K) {
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    int col = blockIdx.x * blockDim.x + threadIdx.x;

    if (row < M && col < N) {
        float sum = 0.0f;
        for (int k = 0; k < K; k++) {
            sum += A[row * K + k] * B[k * N + col];
        }
        C[row * N + col] = sum;
    }
}

// Shared Memory를 활용한 Tiled 행렬 곱셈 (최적화 버전)
__global__ void matMulTiled(const float* A, const float* B, float* C,
                            int M, int N, int K) {
    __shared__ float tileA[TILE_SIZE][TILE_SIZE];
    __shared__ float tileB[TILE_SIZE][TILE_SIZE];

    int row = blockIdx.y * TILE_SIZE + threadIdx.y;
    int col = blockIdx.x * TILE_SIZE + threadIdx.x;
    float sum = 0.0f;

    // 타일 단위로 반복
    for (int t = 0; t < (K + TILE_SIZE - 1) / TILE_SIZE; t++) {
        // Shared Memory로 타일 로드
        if (row < M && (t * TILE_SIZE + threadIdx.x) < K)
            tileA[threadIdx.y][threadIdx.x] = A[row * K + t * TILE_SIZE + threadIdx.x];
        else
            tileA[threadIdx.y][threadIdx.x] = 0.0f;

        if ((t * TILE_SIZE + threadIdx.y) < K && col < N)
            tileB[threadIdx.y][threadIdx.x] = B[(t * TILE_SIZE + threadIdx.y) * N + col];
        else
            tileB[threadIdx.y][threadIdx.x] = 0.0f;

        __syncthreads();  // 타일 로드 완료 대기

        // 타일 내 곱셈-누적
        for (int k = 0; k < TILE_SIZE; k++) {
            sum += tileA[threadIdx.y][k] * tileB[k][threadIdx.x];
        }

        __syncthreads();  // 다음 타일 로드 전 동기화
    }

    if (row < M && col < N) {
        C[row * N + col] = sum;
    }
}

int main() {
    const int M = 1024, N = 1024, K = 1024;

    // ... (메모리 할당, 초기화, 데이터 전송 생략)

    dim3 blockDim(TILE_SIZE, TILE_SIZE);
    dim3 gridDim((N + TILE_SIZE - 1) / TILE_SIZE,
                 (M + TILE_SIZE - 1) / TILE_SIZE);

    matMulTiled<<<gridDim, blockDim>>>(d_A, d_B, d_C, M, N, K);

    // ... (결과 복사, 검증, 메모리 해제 생략)
    return 0;
}

// Coalesced Access (좋은 패턴)
// Warp 내 스레드들이 연속 메모리에 접근
__global__ void coalesced(float* data, int n) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < n) {
        data[idx] = data[idx] * 2.0f;  // 연속 접근
    }
}

// Non-coalesced Access (나쁜 패턴)
// 스레드가 stride를 가지고 비연속적으로 접근
__global__ void strided(float* data, int n, int stride) {
    int idx = (blockIdx.x * blockDim.x + threadIdx.x) * stride;
    if (idx < n) {
        data[idx] = data[idx] * 2.0f;  // stride만큼 건너뛴 접근
    }
}

// AoS(Array of Structures) vs SoA(Structure of Arrays)
// AoS: Non-coalesced (나쁨)
struct ParticleAoS {
    float x, y, z;
    float vx, vy, vz;
};
// 스레드 0: particle[0].x, 스레드 1: particle[1].x -> stride 접근

// SoA: Coalesced (좋음)
struct ParticlesSoA {
    float* x;  float* y;  float* z;
    float* vx; float* vy; float* vz;
};
// 스레드 0: x[0], 스레드 1: x[1] -> 연속 접근

// Shared Memory 사용 예시
__global__ void sharedMemExample(float* input, float* output, int n) {
    // 정적 할당
    __shared__ float sharedData[256];

    int tid = threadIdx.x;
    int gid = blockIdx.x * blockDim.x + threadIdx.x;

    // Global -> Shared 로드
    if (gid < n) {
        sharedData[tid] = input[gid];
    }
    __syncthreads();  // 모든 스레드가 로드 완료할 때까지 대기

    // Shared Memory에서 연산 (이웃 원소 합산 예시)
    if (tid > 0 && tid < blockDim.x - 1 && gid < n) {
        output[gid] = sharedData[tid - 1] + sharedData[tid] + sharedData[tid + 1];
    }
}

// Bank Conflict 회피: 패딩 기법
// 문제: 32x32 행렬의 열(column) 접근 시 모든 스레드가 같은 Bank에 매핑
__shared__ float tile[32][32];       // Bank Conflict 발생!
__shared__ float tile[32][32 + 1];   // 패딩으로 Bank Conflict 회피

// Unified Memory 사용
float* data;
cudaMallocManaged(&data, N * sizeof(float));

// CPU에서 초기화
for (int i = 0; i < N; i++) {
    data[i] = (float)i;
}

// GPU 커널 실행 - 별도 cudaMemcpy 불필요!
myKernel<<<blocks, threads>>>(data, N);
cudaDeviceSynchronize();

// CPU에서 결과 바로 접근
printf("Result: %f\n", data[0]);

// 해제
cudaFree(data);

// 메모리 프리페치 힌트 (성능 최적화)
cudaMemPrefetchAsync(data, N * sizeof(float), deviceId);  // GPU로 프리페치
cudaMemPrefetchAsync(data, N * sizeof(float), cudaCpuDeviceId);  // CPU로 프리페치

// Pinned Memory 할당
float* h_pinned;
cudaMallocHost(&h_pinned, N * sizeof(float));  // 또는 cudaHostAlloc

// 일반 메모리 대비 전송 속도 비교
// Pageable:  ~12 GB/s (PCIe Gen4)
// Pinned:    ~25 GB/s (PCIe Gen4) - 약 2배 빠름

// 비동기 전송에 필수
cudaMemcpyAsync(d_data, h_pinned, bytes, cudaMemcpyHostToDevice, stream);

// 해제
cudaFreeHost(h_pinned);

// Stream-ordered Memory Allocator (CUDA 11.2+)
float* d_data;
cudaMallocAsync(&d_data, bytes, stream);

// 사용 후 비동기 해제
cudaFreeAsync(d_data, stream);

// Memory Pool 설정
cudaMemPool_t mempool;
cudaDeviceGetDefaultMemPool(&mempool, deviceId);

// Pool 크기 제한 설정
uint64_t threshold = 1ULL << 30;  // 1 GB
cudaMemPoolSetAttribute(mempool, cudaMemPoolAttrReleaseThreshold, &threshold);

Default Stream (Stream 0):
[MemcpyH2D] -> [Kernel A] -> [MemcpyD2H]
                                           (순차 실행, 파이프라인 미활용)

Multi-Stream:
Stream 1: [MemcpyH2D_1] -> [Kernel_1] -> [MemcpyD2H_1]
Stream 2:     [MemcpyH2D_2] -> [Kernel_2] -> [MemcpyD2H_2]
Stream 3:         [MemcpyH2D_3] -> [Kernel_3] -> [MemcpyD2H_3]
                                           (오버랩 실행, GPU 활용 극대화)

// Stream 생성
cudaStream_t stream1, stream2;
cudaStreamCreate(&stream1);
cudaStreamCreate(&stream2);

// 비동기 메모리 전송 (Pinned Memory 필수!)
cudaMemcpyAsync(d_A, h_A, bytes, cudaMemcpyHostToDevice, stream1);
cudaMemcpyAsync(d_B, h_B, bytes, cudaMemcpyHostToDevice, stream2);

// 커널 실행 (특정 Stream에서)
kernelA<<<grid, block, 0, stream1>>>(d_A);
kernelB<<<grid, block, 0, stream2>>>(d_B);

// 결과 비동기 복사
cudaMemcpyAsync(h_A, d_A, bytes, cudaMemcpyDeviceToHost, stream1);
cudaMemcpyAsync(h_B, d_B, bytes, cudaMemcpyDeviceToHost, stream2);

// Stream 동기화
cudaStreamSynchronize(stream1);
cudaStreamSynchronize(stream2);

// Stream 해제
cudaStreamDestroy(stream1);
cudaStreamDestroy(stream2);

const int NUM_STREAMS = 4;
const int CHUNK_SIZE = N / NUM_STREAMS;

cudaStream_t streams[NUM_STREAMS];
for (int i = 0; i < NUM_STREAMS; i++) {
    cudaStreamCreate(&streams[i]);
}

// Pinned Memory 할당
float *h_in, *h_out;
cudaMallocHost(&h_in, N * sizeof(float));
cudaMallocHost(&h_out, N * sizeof(float));

float *d_in, *d_out;
cudaMalloc(&d_in, N * sizeof(float));
cudaMalloc(&d_out, N * sizeof(float));

// 파이프라인 실행
for (int i = 0; i < NUM_STREAMS; i++) {
    int offset = i * CHUNK_SIZE;
    size_t chunkBytes = CHUNK_SIZE * sizeof(float);

    // 1. Host -> Device (비동기)
    cudaMemcpyAsync(d_in + offset, h_in + offset,
                    chunkBytes, cudaMemcpyHostToDevice, streams[i]);

    // 2. 커널 실행
    int blocks = (CHUNK_SIZE + 255) / 256;
    processKernel<<<blocks, 256, 0, streams[i]>>>(
        d_in + offset, d_out + offset, CHUNK_SIZE);

    // 3. Device -> Host (비동기)
    cudaMemcpyAsync(h_out + offset, d_out + offset,
                    chunkBytes, cudaMemcpyDeviceToHost, streams[i]);
}

// 전체 완료 대기
cudaDeviceSynchronize();

cudaEvent_t start, stop;
cudaEventCreate(&start);
cudaEventCreate(&stop);

// 타이밍 시작
cudaEventRecord(start);

// 커널 실행
myKernel<<<grid, block>>>(d_data, N);

// 타이밍 종료
cudaEventRecord(stop);
cudaEventSynchronize(stop);

// 경과 시간 계산
float milliseconds = 0;
cudaEventElapsedTime(&milliseconds, start, stop);
printf("Kernel execution time: %.3f ms\n", milliseconds);

// 대역폭 계산
float bandwidth = (bytes_read + bytes_written) / (milliseconds * 1e6);
printf("Effective bandwidth: %.2f GB/s\n", bandwidth);

cudaEventDestroy(start);
cudaEventDestroy(stop);

// Occupancy 계산 API
int blockSize;
int minGridSize;

// 최적 블록 크기 자동 계산
cudaOccupancyMaxPotentialBlockSize(&minGridSize, &blockSize, myKernel, 0, 0);
printf("Optimal block size: %d\n", blockSize);
printf("Minimum grid size: %d\n", minGridSize);

// 특정 블록 크기에서의 Occupancy 확인
int maxActiveBlocks;
cudaOccupancyMaxActiveBlocksPerMultiprocessor(
    &maxActiveBlocks, myKernel, blockSize, 0);
printf("Max active blocks per SM: %d\n", maxActiveBlocks);

// 나쁜 예: Warp Divergence 발생
__global__ void divergentKernel(float* data, int n) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < n) {
        // threadIdx.x가 짝수/홀수에 따라 분기
        // -> 같은 Warp 내에서 절반만 실행
        if (threadIdx.x % 2 == 0) {
            data[idx] = data[idx] * 2.0f;
        } else {
            data[idx] = data[idx] + 1.0f;
        }
    }
}

// 좋은 예: Warp 단위로 분기
__global__ void convergentKernel(float* data, int n) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    int warpId = threadIdx.x / 32;  // Warp ID로 분기
    if (idx < n) {
        // 같은 Warp 내 모든 스레드가 같은 경로
        if (warpId % 2 == 0) {
            data[idx] = data[idx] * 2.0f;
        } else {
            data[idx] = data[idx] + 1.0f;
        }
    }
}

// 수동 Loop Unrolling
__global__ void unrolled(float* data, int n) {
    int idx = blockIdx.x * blockDim.x * 4 + threadIdx.x;

    // 4배 Unrolling
    if (idx < n)         data[idx] *= 2.0f;
    if (idx + 256 < n)   data[idx + 256] *= 2.0f;
    if (idx + 512 < n)   data[idx + 512] *= 2.0f;
    if (idx + 768 < n)   data[idx + 768] *= 2.0f;
}

// 컴파일러 지시자를 이용한 Unrolling
__global__ void pragmaUnrolled(float* A, float* B, float* C, int N) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    float sum = 0.0f;

    #pragma unroll 8
    for (int k = 0; k < N; k++) {
        sum += A[idx * N + k] * B[k];
    }
    C[idx] = sum;
}

__global__ void gridStrideKernel(float* data, int n) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    int stride = gridDim.x * blockDim.x;

    // 각 스레드가 stride만큼 건너뛰며 여러 원소 처리
    for (int i = idx; i < n; i += stride) {
        data[i] = data[i] * 2.0f;
    }
}

// 장점:
// 1. 데이터 크기에 관계없이 동작
// 2. 블록/그리드 크기를 자유롭게 조절 가능
// 3. 스레드당 작업량 증가로 커널 시작 오버헤드 감소

# Nsight Compute (커널 레벨 프로파일링)
ncu --set full -o profile_report ./my_cuda_app

# 특정 커널만 프로파일링
ncu --kernel-name myKernel --launch-skip 0 --launch-count 1 ./my_cuda_app

# Nsight Systems (시스템 레벨 프로파일링)
nsys profile --trace=cuda,nvtx -o timeline_report ./my_cuda_app

# 주요 확인 메트릭
# - SM Throughput: SM 활용률
# - Memory Throughput: 메모리 대역폭 활용률
# - Achieved Occupancy: 실제 Occupancy
# - Warp Stall Reasons: Warp 정지 원인
# - L1/L2 Hit Rate: 캐시 적중률

import pycuda.autoinit
import pycuda.driver as cuda
from pycuda.compiler import SourceModule
import numpy as np

# CUDA 커널 정의 (C 코드를 문자열로 작성)
mod = SourceModule("""
__global__ void multiply(float *dest, float *a, float *b, int n) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < n) {
        dest[idx] = a[idx] * b[idx];
    }
}
""")

multiply = mod.get_function("multiply")

# 데이터 준비
n = 1024 * 1024
a = np.random.randn(n).astype(np.float32)
b = np.random.randn(n).astype(np.float32)
dest = np.zeros_like(a)

# 커널 실행
block_size = 256
grid_size = (n + block_size - 1) // block_size

multiply(
    cuda.Out(dest), cuda.In(a), cuda.In(b), np.int32(n),
    block=(block_size, 1, 1), grid=(grid_size, 1)
)

# 결과 검증
assert np.allclose(dest, a * b)
print("PyCUDA multiply: SUCCESS")

from numba import cuda
import numpy as np
import math

# CUDA 커널 정의
@cuda.jit
def vector_add(a, b, c):
    idx = cuda.grid(1)  # 1D 글로벌 인덱스
    if idx < a.shape[0]:
        c[idx] = a[idx] + b[idx]

# 데이터 준비
n = 1_000_000
a = np.random.randn(n).astype(np.float32)
b = np.random.randn(n).astype(np.float32)
c = np.zeros(n, dtype=np.float32)

# Device로 전송
d_a = cuda.to_device(a)
d_b = cuda.to_device(b)
d_c = cuda.to_device(c)

# 커널 실행
threads_per_block = 256
blocks_per_grid = math.ceil(n / threads_per_block)
vector_add[blocks_per_grid, threads_per_block](d_a, d_b, d_c)

# 결과 가져오기
result = d_c.copy_to_host()
assert np.allclose(result, a + b)
print("Numba CUDA vector_add: SUCCESS")


# 2D 커널 예시: 행렬 곱셈
@cuda.jit
def matmul_kernel(A, B, C):
    row, col = cuda.grid(2)  # 2D 글로벌 인덱스
    if row < C.shape[0] and col < C.shape[1]:
        tmp = 0.0
        for k in range(A.shape[1]):
            tmp += A[row, k] * B[k, col]
        C[row, col] = tmp

# Shared Memory 사용
@cuda.jit
def matmul_shared(A, B, C):
    TILE = 16
    sA = cuda.shared.array(shape=(TILE, TILE), dtype=np.float32)
    sB = cuda.shared.array(shape=(TILE, TILE), dtype=np.float32)

    tx = cuda.threadIdx.x
    ty = cuda.threadIdx.y
    row = cuda.blockIdx.y * TILE + ty
    col = cuda.blockIdx.x * TILE + tx

    tmp = 0.0
    for t in range((A.shape[1] + TILE - 1) // TILE):
        if row < A.shape[0] and (t * TILE + tx) < A.shape[1]:
            sA[ty, tx] = A[row, t * TILE + tx]
        else:
            sA[ty, tx] = 0.0

        if (t * TILE + ty) < B.shape[0] and col < B.shape[1]:
            sB[ty, tx] = B[t * TILE + ty, col]
        else:
            sB[ty, tx] = 0.0

        cuda.syncthreads()
        for k in range(TILE):
            tmp += sA[ty, k] * sB[k, tx]
        cuda.syncthreads()

    if row < C.shape[0] and col < C.shape[1]:
        C[row, col] = tmp

import cupy as cp
import numpy as np
import time

# NumPy와 동일한 API
a_gpu = cp.random.randn(10000, 10000, dtype=cp.float32)
b_gpu = cp.random.randn(10000, 10000, dtype=cp.float32)

# GPU 행렬 곱셈
start = time.time()
c_gpu = cp.dot(a_gpu, b_gpu)
cp.cuda.Stream.null.synchronize()  # GPU 작업 완료 대기
gpu_time = time.time() - start

# CPU 비교
a_cpu = cp.asnumpy(a_gpu)
b_cpu = cp.asnumpy(b_gpu)
start = time.time()
c_cpu = np.dot(a_cpu, b_cpu)
cpu_time = time.time() - start

print(f"GPU: {gpu_time:.4f}s, CPU: {cpu_time:.4f}s")
print(f"Speedup: {cpu_time / gpu_time:.1f}x")

# CuPy에서 커스텀 CUDA 커널 사용
custom_kernel = cp.RawKernel(r'''
extern "C" __global__
void relu_kernel(const float* input, float* output, int n) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < n) {
        output[idx] = fmaxf(input[idx], 0.0f);
    }
}
''', 'relu_kernel')

n = 1_000_000
x = cp.random.randn(n, dtype=cp.float32)
y = cp.zeros(n, dtype=cp.float32)

block_size = 256
grid_size = (n + block_size - 1) // block_size
custom_kernel((grid_size,), (block_size,), (x, y, n))

# 검증
assert cp.allclose(y, cp.maximum(x, 0))

import torch
import torch.nn as nn

# ===== 기본 GPU 사용 =====

# GPU 사용 가능 확인
print(f"CUDA available: {torch.cuda.is_available()}")
print(f"CUDA version: {torch.version.cuda}")
print(f"GPU count: {torch.cuda.device_count()}")
print(f"Current device: {torch.cuda.current_device()}")
print(f"GPU name: {torch.cuda.get_device_name(0)}")

# 텐서를 GPU로 이동
device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')

# 방법 1: .to(device)
x = torch.randn(1000, 1000)
x_gpu = x.to(device)

# 방법 2: .cuda()
x_gpu = x.cuda()

# 방법 3: 직접 GPU에서 생성
x_gpu = torch.randn(1000, 1000, device='cuda')

# 특정 GPU 지정
x_gpu1 = x.to('cuda:1')  # 두 번째 GPU

# ===== 메모리 관리 =====

# 메모리 사용량 확인
print(torch.cuda.memory_summary())
print(f"Allocated: {torch.cuda.memory_allocated() / 1e9:.2f} GB")
print(f"Cached: {torch.cuda.memory_reserved() / 1e9:.2f} GB")

# 캐시 비우기
torch.cuda.empty_cache()

# 메모리 스냅샷 (디버깅용)
torch.cuda.memory._record_memory_history()
# ... 코드 실행 ...
snapshot = torch.cuda.memory._snapshot()
torch.cuda.memory._dump_snapshot("memory_snapshot.pickle")


# ===== Mixed Precision Training =====
from torch.amp import autocast, GradScaler

model = nn.Linear(1024, 1024).cuda()
optimizer = torch.optim.Adam(model.parameters(), lr=1e-3)
scaler = GradScaler('cuda')

for data, target in dataloader:
    data, target = data.cuda(), target.cuda()
    optimizer.zero_grad()

    # autocas로 FP16 자동 전환
    with autocast('cuda'):
        output = model(data)
        loss = nn.functional.mse_loss(output, target)

    # GradScaler로 FP16 그래디언트 스케일링
    scaler.scale(loss).backward()
    scaler.step(optimizer)
    scaler.update()

# Mixed Precision의 이점:
# - 메모리 사용량 약 50% 감소
# - 학습 속도 1.5~3배 향상 (Tensor Core 활용)
# - 정확도 손실 최소

import tensorflow as tf

# GPU 확인
print("GPUs available:", tf.config.list_physical_devices('GPU'))

# GPU 메모리 점진적 할당 (권장)
gpus = tf.config.list_physical_devices('GPU')
for gpu in gpus:
    tf.config.experimental.set_memory_growth(gpu, True)

# 또는 메모리 제한 설정
tf.config.set_logical_device_configuration(
    gpus[0],
    [tf.config.LogicalDeviceConfiguration(memory_limit=8192)]  # 8GB
)

# 특정 GPU 사용
with tf.device('/GPU:0'):
    a = tf.random.normal([1000, 1000])
    b = tf.random.normal([1000, 1000])
    c = tf.matmul(a, b)

# Mixed Precision
tf.keras.mixed_precision.set_global_policy('mixed_float16')

# cuDF: pandas 대체 GPU DataFrame
import cudf

# pandas와 동일한 API
df = cudf.read_csv('large_data.csv')
result = df.groupby('category').agg({'value': 'mean', 'count': 'sum'})
filtered = df[df['value'] > 100]

# pandas와 상호 변환
import pandas as pd
pdf = df.to_pandas()   # cuDF -> pandas
gdf = cudf.from_pandas(pdf)  # pandas -> cuDF


# cuML: scikit-learn 대체 GPU ML
import cuml
from cuml.ensemble import RandomForestClassifier
from cuml.cluster import KMeans

# scikit-learn과 동일한 API
rf = RandomForestClassifier(n_estimators=100, max_depth=16)
rf.fit(X_train, y_train)
predictions = rf.predict(X_test)

# Zero-code-change 가속 (cuml.accel)
# scikit-learn 코드를 수정 없이 GPU로 실행
import cuml.accel
cuml.accel.install()

from sklearn.ensemble import RandomForestClassifier  # 자동으로 GPU 가속!

from numba import cuda, njit
import numpy as np
import time
import math

# CPU 버전 (Numba JIT)
@njit
def vector_ops_cpu(a, b, c):
    for i in range(a.shape[0]):
        c[i] = math.sqrt(a[i] ** 2 + b[i] ** 2) * math.sin(a[i]) + math.log(abs(b[i]) + 1)

# GPU 버전 (CUDA)
@cuda.jit
def vector_ops_gpu(a, b, c):
    idx = cuda.grid(1)
    if idx < a.shape[0]:
        c[idx] = math.sqrt(a[idx] ** 2 + b[idx] ** 2) * math.sin(a[idx]) + math.log(abs(b[idx]) + 1)

# 벤치마크
sizes = [100_000, 1_000_000, 10_000_000, 100_000_000]

for n in sizes:
    a = np.random.randn(n).astype(np.float32)
    b = np.random.randn(n).astype(np.float32)
    c_cpu = np.zeros(n, dtype=np.float32)
    c_gpu = np.zeros(n, dtype=np.float32)

    # CPU
    start = time.time()
    vector_ops_cpu(a, b, c_cpu)
    cpu_time = time.time() - start

    # GPU
    d_a = cuda.to_device(a)
    d_b = cuda.to_device(b)
    d_c = cuda.to_device(c_gpu)

    threads = 256
    blocks = math.ceil(n / threads)

    # Warm-up
    vector_ops_gpu[blocks, threads](d_a, d_b, d_c)
    cuda.synchronize()

    start = time.time()
    vector_ops_gpu[blocks, threads](d_a, d_b, d_c)
    cuda.synchronize()
    gpu_time = time.time() - start

    c_gpu = d_c.copy_to_host()

    print(f"N={n:>12,}: CPU={cpu_time:.4f}s, GPU={gpu_time:.4f}s, "
          f"Speedup={cpu_time/gpu_time:.1f}x")

import cupy as cp
import numpy as np
from PIL import Image

def gpu_gaussian_blur(image_array, kernel_size=5, sigma=1.0):
    """GPU 가우시안 블러 구현"""
    # 가우시안 커널 생성
    ax = cp.arange(-kernel_size // 2 + 1., kernel_size // 2 + 1.)
    xx, yy = cp.meshgrid(ax, ax)
    kernel = cp.exp(-(xx**2 + yy**2) / (2. * sigma**2))
    kernel = kernel / kernel.sum()

    img_gpu = cp.asarray(image_array, dtype=cp.float32)
    result = cp.zeros_like(img_gpu)

    pad = kernel_size // 2
    # 채널별 처리
    for c in range(img_gpu.shape[2]):
        padded = cp.pad(img_gpu[:, :, c], pad, mode='reflect')
        # 2D Convolution (CuPy FFT 활용)
        from cupyx.scipy.ndimage import convolve
        result[:, :, c] = convolve(img_gpu[:, :, c], kernel)

    return cp.asnumpy(result.clip(0, 255).astype(cp.uint8))


# CuPy RawKernel로 직접 구현한 Sobel Edge Detection
sobel_kernel = cp.RawKernel(r'''
extern "C" __global__
void sobel_filter(const float* input, float* output,
                  int width, int height) {
    int x = blockIdx.x * blockDim.x + threadIdx.x;
    int y = blockIdx.y * blockDim.y + threadIdx.y;

    if (x > 0 && x < width - 1 && y > 0 && y < height - 1) {
        // Sobel X
        float gx = -1.0f * input[(y-1)*width + (x-1)]
                   + 1.0f * input[(y-1)*width + (x+1)]
                   - 2.0f * input[y*width + (x-1)]
                   + 2.0f * input[y*width + (x+1)]
                   - 1.0f * input[(y+1)*width + (x-1)]
                   + 1.0f * input[(y+1)*width + (x+1)];

        // Sobel Y
        float gy = -1.0f * input[(y-1)*width + (x-1)]
                   - 2.0f * input[(y-1)*width + x]
                   - 1.0f * input[(y-1)*width + (x+1)]
                   + 1.0f * input[(y+1)*width + (x-1)]
                   + 2.0f * input[(y+1)*width + x]
                   + 1.0f * input[(y+1)*width + (x+1)];

        output[y * width + x] = sqrtf(gx * gx + gy * gy);
    }
}
''', 'sobel_filter')

def gpu_sobel_edge_detection(gray_image):
    """GPU Sobel 엣지 검출"""
    h, w = gray_image.shape
    img_gpu = cp.asarray(gray_image, dtype=cp.float32)
    output_gpu = cp.zeros_like(img_gpu)

    block = (16, 16)
    grid = ((w + 15) // 16, (h + 15) // 16)

    sobel_kernel(grid, block, (img_gpu, output_gpu, w, h))

    return cp.asnumpy(output_gpu.clip(0, 255).astype(cp.uint8))

# 사용 예
# img = np.array(Image.open('photo.jpg'))
# blurred = gpu_gaussian_blur(img, kernel_size=11, sigma=3.0)
# gray = np.mean(img, axis=2).astype(np.float32)
# edges = gpu_sobel_edge_detection(gray)

import torch
import torch.nn as nn
import torch.optim as optim
from torch.utils.data import DataLoader, TensorDataset
from torch.amp import autocast, GradScaler
import time

# 디바이스 설정
device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')
print(f"Using device: {device}")
if device.type == 'cuda':
    print(f"GPU: {torch.cuda.get_device_name(0)}")

# 모델 정의
class ResidualBlock(nn.Module):
    def __init__(self, dim):
        super().__init__()
        self.net = nn.Sequential(
            nn.Linear(dim, dim),
            nn.LayerNorm(dim),
            nn.GELU(),
            nn.Linear(dim, dim),
            nn.LayerNorm(dim),
        )

    def forward(self, x):
        return x + self.net(x)

class DeepModel(nn.Module):
    def __init__(self, input_dim=784, hidden_dim=512, num_classes=10, num_blocks=6):
        super().__init__()
        self.input_proj = nn.Linear(input_dim, hidden_dim)
        self.blocks = nn.Sequential(*[ResidualBlock(hidden_dim) for _ in range(num_blocks)])
        self.head = nn.Linear(hidden_dim, num_classes)

    def forward(self, x):
        x = self.input_proj(x)
        x = self.blocks(x)
        return self.head(x)

# 데이터 생성 (실제로는 데이터셋 로드)
X_train = torch.randn(50000, 784)
y_train = torch.randint(0, 10, (50000,))
train_dataset = TensorDataset(X_train, y_train)
train_loader = DataLoader(train_dataset, batch_size=256, shuffle=True,
                          num_workers=4, pin_memory=True)  # pin_memory 중요!

# 모델을 GPU로 이동
model = DeepModel().to(device)
criterion = nn.CrossEntropyLoss()
optimizer = optim.AdamW(model.parameters(), lr=1e-3, weight_decay=0.01)
scaler = GradScaler('cuda')

# 학습 루프 (Mixed Precision)
def train_epoch(model, loader, criterion, optimizer, scaler, device):
    model.train()
    total_loss = 0
    correct = 0
    total = 0

    for batch_idx, (data, target) in enumerate(loader):
        data, target = data.to(device, non_blocking=True), target.to(device, non_blocking=True)

        optimizer.zero_grad(set_to_none=True)  # set_to_none=True가 더 효율적

        with autocast('cuda'):
            output = model(data)
            loss = criterion(output, target)

        scaler.scale(loss).backward()
        scaler.step(optimizer)
        scaler.update()

        total_loss += loss.item()
        _, predicted = output.max(1)
        total += target.size(0)
        correct += predicted.eq(target).sum().item()

    return total_loss / len(loader), 100. * correct / total

# 학습 실행
num_epochs = 10
for epoch in range(num_epochs):
    start = time.time()
    loss, acc = train_epoch(model, train_loader, criterion, optimizer, scaler, device)
    elapsed = time.time() - start

    print(f"Epoch {epoch+1}/{num_epochs}: Loss={loss:.4f}, "
          f"Acc={acc:.2f}%, Time={elapsed:.2f}s")

    # GPU 메모리 상태 출력
    if device.type == 'cuda':
        alloc = torch.cuda.memory_allocated() / 1e9
        reserved = torch.cuda.memory_reserved() / 1e9
        print(f"  GPU Memory: {alloc:.2f} GB allocated, {reserved:.2f} GB reserved")

import torch
import torch.nn as nn

model = DeepModel().cuda()

# DataParallel 적용 (한 줄로 끝!)
if torch.cuda.device_count() > 1:
    print(f"Using {torch.cuda.device_count()} GPUs with DataParallel")
    model = nn.DataParallel(model)

# 이후 학습 코드는 동일
# 단점: GPU 0에 부하가 집중되는 불균형 문제

import torch
import torch.nn as nn
import torch.distributed as dist
from torch.nn.parallel import DistributedDataParallel as DDP
from torch.utils.data.distributed import DistributedSampler
import os

def setup(rank, world_size):
    os.environ['MASTER_ADDR'] = 'localhost'
    os.environ['MASTER_PORT'] = '12355'
    dist.init_process_group("nccl", rank=rank, world_size=world_size)
    torch.cuda.set_device(rank)

def cleanup():
    dist.destroy_process_group()

def train_ddp(rank, world_size):
    setup(rank, world_size)

    # 모델을 해당 GPU에 배치
    model = DeepModel().to(rank)
    ddp_model = DDP(model, device_ids=[rank])

    # 분산 데이터 로더
    train_dataset = TensorDataset(X_train, y_train)
    sampler = DistributedSampler(train_dataset, num_replicas=world_size, rank=rank)
    train_loader = DataLoader(train_dataset, batch_size=256,
                              sampler=sampler, num_workers=4, pin_memory=True)

    optimizer = optim.AdamW(ddp_model.parameters(), lr=1e-3)
    criterion = nn.CrossEntropyLoss()
    scaler = GradScaler('cuda')

    for epoch in range(10):
        sampler.set_epoch(epoch)  # 매 에폭마다 셔플링
        for data, target in train_loader:
            data = data.to(rank, non_blocking=True)
            target = target.to(rank, non_blocking=True)

            optimizer.zero_grad(set_to_none=True)
            with autocast('cuda'):
                output = ddp_model(data)
                loss = criterion(output, target)

            scaler.scale(loss).backward()
            scaler.step(optimizer)
            scaler.update()

        if rank == 0:
            print(f"Epoch {epoch+1} complete")

    cleanup()

# 실행
# torchrun --nproc_per_node=4 train_ddp.py
import torch.multiprocessing as mp
world_size = torch.cuda.device_count()
mp.spawn(train_ddp, args=(world_size,), nprocs=world_size, join=True)

import pycuda.autoinit
import pycuda.driver as cuda
from pycuda.compiler import SourceModule
import numpy as np

# Parallel Reduction 커널
mod = SourceModule("""
__global__ void parallel_reduce(float *input, float *output, int n) {
    extern __shared__ float sdata[];

    unsigned int tid = threadIdx.x;
    unsigned int i = blockIdx.x * blockDim.x * 2 + threadIdx.x;

    // 첫 번째 레벨 리덕션 (글로벌 로드 시 수행)
    float sum = 0.0f;
    if (i < n) sum += input[i];
    if (i + blockDim.x < n) sum += input[i + blockDim.x];
    sdata[tid] = sum;

    __syncthreads();

    // Shared Memory에서 리덕션
    for (unsigned int s = blockDim.x / 2; s > 32; s >>= 1) {
        if (tid < s) {
            sdata[tid] += sdata[tid + s];
        }
        __syncthreads();
    }

    // Warp 내에서는 동기화 불필요 (Warp-level reduction)
    if (tid < 32) {
        volatile float *smem = sdata;
        smem[tid] += smem[tid + 32];
        smem[tid] += smem[tid + 16];
        smem[tid] += smem[tid + 8];
        smem[tid] += smem[tid + 4];
        smem[tid] += smem[tid + 2];
        smem[tid] += smem[tid + 1];
    }

    if (tid == 0) {
        output[blockIdx.x] = sdata[0];
    }
}
""")

reduce_kernel = mod.get_function("parallel_reduce")

# 실행
n = 1 << 20  # 1M 원소
data = np.random.randn(n).astype(np.float32)

block_size = 256
grid_size = (n + block_size * 2 - 1) // (block_size * 2)

d_input = cuda.mem_alloc(data.nbytes)
d_output = cuda.mem_alloc(grid_size * 4)  # float32

cuda.memcpy_htod(d_input, data)

reduce_kernel(d_input, d_output, np.int32(n),
              block=(block_size, 1, 1), grid=(grid_size, 1),
              shared=block_size * 4)

# 부분 합계 가져오기
partial_sums = np.zeros(grid_size, dtype=np.float32)
cuda.memcpy_dtoh(partial_sums, d_output)
gpu_sum = partial_sums.sum()

# 검증
cpu_sum = data.sum()
print(f"GPU sum: {gpu_sum:.6f}")
print(f"CPU sum: {cpu_sum:.6f}")
print(f"Difference: {abs(gpu_sum - cpu_sum):.6f}")