Chapter 8

스텐실: 미분방정식 풀이의 GPU 처세술

3D 격자에서 시간을 굴리는, 면적/부피의 싸움

8.1 배경: PDE의 이산화와 스텐실 패턴

유체역학 시뮬레이션, 열전도, 음향 전파, 자기장 — 자연 현상을 기술하는 편미분방정식(PDE)을 컴퓨터로 풀려면 연속 공간을 격자(grid)로 잘라 미분을 차분으로 바꿔야 한다. 이를 유한 차분(finite difference)이라 부른다.

가장 흔한 예가 라플라시안(Laplacian) ∇²u를 격자에서 근사하는 7-point stencil이다. 격자점 (i,j,k)에서의 라플라시안은 자기 자신과 6개 이웃의 가중합으로 계산된다.

∇²u[i,j,k] ≈ (u[i±1,j,k] + u[i,j±1,k] + u[i,j,k±1] - 6·u[i,j,k]) / h²

이런 패턴 — 격자점 하나의 새 값이 자기와 인접한 몇 개 점의 가중합으로 결정되는 — 을 stencil이라 부른다. 2D에서는 5-point stencil(자신 + 상하좌우 4)가, 3D에서는 7-point(자신 + 6면 이웃)나 27-point(자신 + 모든 인접 입체) 같은 변형이 흔하다.

컨볼루션과의 차이

구조적으로 stencil은 컨볼루션의 사촌이다. 둘 다 출력 한 점이 입력 윈도우의 가중합이다. 그러나 사용 양상이 사뭇 다르다.

반복. 컨볼루션은 한 번의 sweep으로 끝난다. stencil은 시간 step마다 같은 격자를 update하며 수백~수만 번 반복한다.
차원. 컨볼루션은 보통 2D 이미지가 주력. stencil은 3D 물리 격자가 일상적이다. 이게 의외로 큰 차이를 만든다.
정밀도. 시뮬레이션이 시간으로 누적되므로 작은 오차도 발산할 수 있다. 보통 double precision을 쓰고, 산술 강도가 그만큼 더 낮아진다.
가중치 sparsity. 7-point stencil의 25-point 윈도우 중 6개만 nonzero. 컨볼루션 필터는 빽빽한 반면 stencil 필터는 듬성듬성하다. 그래서 stencil은 가중치 메모리 자체보다 입력 메모리 트래픽이 지배적이다.

8.2 기본 병렬 스텐실

매핑은 익숙하다. 출력 격자점 1개당 스레드 1개. 3D 그리드 차원의 블록과 그리드를 띄운다.

__global__ void laplacian7pt(const double* in, double* out,
                              int Nx, int Ny, int Nz, double cself, double cnb) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    int j = blockIdx.y * blockDim.y + threadIdx.y;
    int k = blockIdx.z * blockDim.z + threadIdx.z;

    if (i < 1 || i >= Nx - 1 || j < 1 || j >= Ny - 1 || k < 1 || k >= Nz - 1)
        return;   // 경계는 별도 처리

    int idx = (k * Ny + j) * Nx + i;
    int ip1 = idx + 1,    im1 = idx - 1;
    int jp1 = idx + Nx,   jm1 = idx - Nx;
    int kp1 = idx + Nx*Ny, km1 = idx - Nx*Ny;

    out[idx] = cself * in[idx]
             + cnb * (in[ip1] + in[im1] + in[jp1] + in[jm1] + in[kp1] + in[km1]);
}

분석. 출력 1점당 입력 7개를 읽고(56 byte, double 8byte × 7), 7번의 곱셈과 6번의 덧셈, 즉 13 FLOP. 산술 강도는 13 / 56 ≈ 0.23 FLOP/B. 7장의 컨볼루션과 비슷하게 메모리 바운드 그 자체다.

다행히 이웃한 출력점은 입력을 광범위하게 공유한다. (i,j,k)와 (i+1,j,k)는 in[i,j,k]를 둘 다 읽는다. 한 격자선을 따라 보면 출력 1점당 진짜로 새로 읽어야 할 입력은 한두 개뿐이다. 타일링과 코어스닝이 통할 풍경이다.

8.3 스텐실용 shared memory 타일링

2D 컨볼루션과 똑같이 입력 타일을 공유 메모리에 협력 로드해 보자. 3D 출력 타일이 O × O × O라면 필요한 입력 타일은 가장자리에 r=1만큼씩 늘려 (O+2) × (O+2) × (O+2)다.

타일링 효율을 계산하자. 출력 타일에서 글로벌 트래픽 절감 비율은

효율 = O³ × 7 / (O + 2)³

O=8이면 512 × 7 / 1000 = 3.58배. O=16이면 4096 × 7 / 5832 ≈ 4.92배. 비교를 위해 7장 2D 컨볼루션 r=1, O=16의 효율은 256 × 9 / 324 ≈ 7.1배였다. 3D는 같은 O에서도 효율이 더 낮다.

함정

왜 3D가 불리한가? 면적 vs 부피의 싸움이다. halo는 부피의 표면적에 해당한다. 출력 부피는 O³로 자라고 입력 부피는 (O+2)³로 자라는데, 이 비율은 차원이 높을수록 1에 더 가까워지는 게 아니라 halo 비중이 더 커진다. 한 변 O가 작으면 표면적이 부피와 비슷한 차수가 된다. 게다가 공유 메모리 한도가 부피로 소모되니 O를 크게 키우기도 어렵다. O=16 double 입력 타일이면 18³ × 8 = 46.6 KB — 거의 SM의 공유 메모리 절반이다. 한 SM에 한 블록밖에 못 띄운다.

#define TILE 8
__global__ void laplacian7pt_tiled(const double* in, double* out,
                                    int Nx, int Ny, int Nz,
                                    double cself, double cnb) {
    __shared__ double s[TILE + 2][TILE + 2][TILE + 2];

    int tx = threadIdx.x, ty = threadIdx.y, tz = threadIdx.z;
    int i = blockIdx.x * TILE + tx;
    int j = blockIdx.y * TILE + ty;
    int k = blockIdx.z * TILE + tz;

    // 한 스레드가 한 점만 로드할 수도 있고, 가장자리 처리를 위해
    // 협력 로드 스킴을 짜야 한다. 여기서는 단순화를 위해
    // 블록 차원을 (TILE+2)^3 으로 잡았다고 가정.
    int gi = blockIdx.x * TILE + tx - 1;
    int gj = blockIdx.y * TILE + ty - 1;
    int gk = blockIdx.z * TILE + tz - 1;

    s[tz][ty][tx] = (gi >= 0 && gi < Nx && gj >= 0 && gj < Ny && gk >= 0 && gk < Nz)
                  ? in[(gk * Ny + gj) * Nx + gi] : 0.0;
    __syncthreads();

    if (tx >= 1 && tx < TILE + 1 && ty >= 1 && ty < TILE + 1 && tz >= 1 && tz < TILE + 1
        && gi < Nx - 1 && gj < Ny - 1 && gk < Nz - 1
        && gi > 0    && gj > 0    && gk > 0) {
        double v = cself * s[tz][ty][tx]
                 + cnb * (s[tz][ty][tx+1] + s[tz][ty][tx-1]
                        + s[tz][ty+1][tx] + s[tz][ty-1][tx]
                        + s[tz+1][ty][tx] + s[tz-1][ty][tx]);
        out[(gk * Ny + gj) * Nx + gi] = v;
    }
}

이 커널은 잘 작동하지만 위에서 본 효율 한계와 공유 메모리 부담 때문에 실제 성능이 기대만큼 안 나온다. 더 나은 전략이 필요하다.

8.4 스레드 코어스닝: z축 슬랩 스위프

3D 큐브 타일이 비효율적이라면 한 차원만 큐브를 쓰지 말자. 아이디어: xy 평면만 2D 타일로 잡고, z축은 한 스레드가 여러 슬랩을 sweep한다. 한 블록은 이제 (O × O × 1) 모양이고, 한 스레드가 z=0,1,2,...,Nz-1까지 한 칸씩 이동하며 출력을 누적한다.

이 방식의 장점은 명료하다. (1) 공유 메모리 부담이 부피에서 면적으로 줄어든다. O=16이면 (O+2)² = 324개의 double = 2.6 KB. 큐브 타일의 46 KB보다 한 자릿수 적다. (2) z 방향 재사용이 레지스터로 들어간다. 한 슬랩의 가운데 평면을 처리할 때 이전 평면(z-1)과 다음 평면(z+1)이 필요한데, 다음 슬랩으로 넘어가면 현재 평면이 새 z-1, 다음 z+1만 새로 읽으면 된다. (3) 면적/부피 비율이 좋아진다. halo는 xy 평면의 가장자리만 있고 z 방향은 두 평면(앞/뒤)뿐.

z 방향 스위프 ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ iter k=2: prev = plane[1] curr = plane[2] next = plane[3] iter k=3: prev = plane[2] curr = plane[3] next = plane[4] ← prev/curr 슬라이딩 iter k=4: prev = plane[3] curr = plane[4] next = plane[5] 공유 메모리: 한 평면 (O+2)² 만 잡고 매 iter마다 갱신 레지스터: 각 스레드가 prev, curr, next의 자기 위치 값을 보관

그림 8.1 — 코어스닝된 z 슬랩 sweep의 동작

__global__ void laplacian7pt_coarsened(const double* in, double* out,
                                        int Nx, int Ny, int Nz,
                                        double cself, double cnb) {
    __shared__ double s_curr[TILE + 2][TILE + 2];
    __shared__ double s_next[TILE + 2][TILE + 2];

    int tx = threadIdx.x, ty = threadIdx.y;
    int i = blockIdx.x * TILE + tx - 1;
    int j = blockIdx.y * TILE + ty - 1;

    bool inX = (i >= 0 && i < Nx);
    bool inY = (j >= 0 && j < Ny);

    // k=0 평면을 prev로, k=1 평면을 curr로 초기화
    double prev = (inX && inY) ? in[(0 * Ny + j) * Nx + i] : 0.0;
    s_curr[ty][tx] = (inX && inY && Nz > 1) ? in[(1 * Ny + j) * Nx + i] : 0.0;
    __syncthreads();

    for (int k = 1; k < Nz - 1; ++k) {
        // next 평면 로드
        s_next[ty][tx] = (inX && inY) ? in[((k + 1) * Ny + j) * Nx + i] : 0.0;
        __syncthreads();

        // 안쪽 스레드만 출력 계산
        if (tx >= 1 && tx < TILE + 1 && ty >= 1 && ty < TILE + 1
            && i > 0 && i < Nx - 1 && j > 0 && j < Ny - 1) {
            double v = cself * s_curr[ty][tx]
                     + cnb * (s_curr[ty][tx+1] + s_curr[ty][tx-1]
                            + s_curr[ty+1][tx] + s_curr[ty-1][tx]
                            + prev + s_next[ty][tx]);
            out[(k * Ny + j) * Nx + i] = v;
        }
        __syncthreads();

        // 슬라이딩: prev←curr, curr←next, next는 다음 iter에 새로 읽음
        prev = s_curr[ty][tx];
        s_curr[ty][tx] = s_next[ty][tx];
        __syncthreads();
    }
}

한 스레드가 Nz-2번의 z slab를 sweep하므로 그리드의 z 차원이 사라졌다. blockDim도 2D, 그리드도 2D다. 글로벌 트래픽 측면에서는 한 평면을 정확히 한 번씩 읽으면 되므로 N×N×N 격자의 입력 트래픽이 N³ × 8 byte로 이론 하한에 거의 닿는다.

8.5 레지스터 타일링: prev를 공유 메모리에서 빼기

위 코드를 한 단계 더 다듬을 수 있다. 8.4의 prev는 사실 평면 전체를 공유 메모리에 두고 있는 게 아니라 각 스레드가 자기 위치의 값만 있으면 된다. 라플라시안 계산에서 prev 평면은 자기 위치 한 점만 쓰이지, 옆 점은 안 쓰이기 때문이다. (xy 방향 이웃은 curr 평면에서, z 방향 이웃은 prev/next의 자기 위치에서.)

그러므로 prev와 next도 공유 메모리에 두지 말고 레지스터에 두자. 위 코드는 이미 prev를 레지스터(double prev)에 두고 있어서 그 부분은 OK다. next도 같은 방식으로 처리하면 공유 메모리는 curr 평면 하나만 들고 가면 된다.

__global__ void laplacian7pt_regTiled(const double* in, double* out,
                                       int Nx, int Ny, int Nz,
                                       double cself, double cnb) {
    __shared__ double s_curr[TILE + 2][TILE + 2];

    int tx = threadIdx.x, ty = threadIdx.y;
    int i = blockIdx.x * TILE + tx - 1;
    int j = blockIdx.y * TILE + ty - 1;
    bool inX = (i >= 0 && i < Nx);
    bool inY = (j >= 0 && j < Ny);

    double prev = (inX && inY) ? in[(0 * Ny + j) * Nx + i] : 0.0;
    double curr = (inX && inY && Nz > 1) ? in[(1 * Ny + j) * Nx + i] : 0.0;
    double next;

    s_curr[ty][tx] = curr;
    __syncthreads();

    for (int k = 1; k < Nz - 1; ++k) {
        next = (inX && inY) ? in[((k + 1) * Ny + j) * Nx + i] : 0.0;

        if (tx >= 1 && tx < TILE + 1 && ty >= 1 && ty < TILE + 1
            && i > 0 && i < Nx - 1 && j > 0 && j < Ny - 1) {
            double v = cself * curr
                     + cnb * (s_curr[ty][tx+1] + s_curr[ty][tx-1]
                            + s_curr[ty+1][tx] + s_curr[ty-1][tx]
                            + prev + next);
            out[(k * Ny + j) * Nx + i] = v;
        }

        // 레지스터 슬라이딩 윈도우
        prev = curr;
        curr = next;
        __syncthreads();
        s_curr[ty][tx] = curr;
        __syncthreads();
    }
}

공유 메모리 사용량이 (TILE+2)² × 8 바이트로 줄어들었다 (TILE=16이면 약 2.6KB). 한 SM에 동시 상주하는 블록 수가 늘어 occupancy가 회복된다. 레지스터는 prev·curr·next 3개씩 늘었지만 보통 GPU는 이 정도는 여유가 있다. 이 패턴 — z 방향 이웃을 레지스터에 슬라이딩 윈도우로 보관 — 은 stencil 코드의 거의 표준 idiom이다.

잘 된 사례

3D stencil은 "차원이 늘 때 어디에 데이터를 두느냐"에서 성능이 결정된다. 좋은 stencil 커널의 패턴: (1) 가장 큰 차원을 sweep으로 변환, (2) 평면 한 장만 공유 메모리에, (3) sweep 차원의 이웃은 레지스터에. 이 셋이 맞물리면 트래픽은 이론 하한에, 공유 메모리는 작은 면적에 묶인다.

8.6 정리

스텐실은 컨볼루션의 사촌이지만 3D 격자와 시간 반복이라는 두 가지 어려움을 더 안고 있다. 단순 큐브 타일링은 면적/부피 비율 때문에 효율이 낮고 공유 메모리 부담이 크다. 해법은 한 차원을 sweep으로 바꿔서 한 스레드가 여러 슬랩을 처리하는 코어스닝, 그리고 sweep 차원의 이웃을 레지스터에 슬라이딩으로 잡는 레지스터 타일링이다. 이 두 기법을 합치면 공유 메모리는 평면 한 장으로 줄고, 글로벌 트래픽은 격자 한 번 읽기 수준에 도달한다.

실제 PDE 시뮬레이션에서는 시간 step이 수만 번 반복되므로 한 step의 작은 비효율이 누적되어 큰 차이를 만든다. 8장의 코어스닝 + 레지스터 타일링은 단순 버전 대비 보통 3~5배의 속도 향상을 가져온다. 다음 장은 또 다른 종류의 redundancy 패턴, 히스토그램과 atomic 연산을 다룬다. 거기서는 메모리가 아니라 경합이 새로운 적이다.

이 챕터에서 챙길 것

stencil = 격자점이 자신과 인접 점들의 가중합으로 update되는 패턴. PDE 풀이의 기본 도구.
컨볼루션과 닮았지만 3D, 반복 sweep, double precision이 일반적이라 더 까다롭다.
3D 큐브 타일링은 halo 부피 비중이 커서 효율이 낮다 (O=16, r=1에서 약 5배 절감에 그침).
해법은 z 슬랩 sweep 코어스닝: 한 스레드가 z 차원을 sweep하면 공유 메모리는 평면 1장으로, 효율은 이론 하한에 가까워진다.
레지스터 타일링: prev/curr/next의 자기 위치 값만 레지스터에 슬라이딩 윈도우로 보관해 공유 메모리 압력을 더 줄인다.
잘 짠 3D stencil 패턴: 큰 차원 sweep → 공유 메모리는 평면 → sweep 차원 이웃은 레지스터.