스텐실 연산

1. 개요

스텐실 연산 Stencil computation
정의	격자점마다 고정된 이웃 패턴의 가중합으로 새 값을 계산
대표	2D 5점·3D 7점 라플라시안, 3D 27점, 고차 25점 별
병목	거의 항상 메모리 대역폭
3D 7점 산술강도	약 0.33 flop/B (배정밀도, 쓰기 할당 포함)
최적화	레이어 조건·공간 블로킹·시간 블로킹·벡터화
쓰이는 곳	유한차분, LBM, 영상 필터, 셀룰러 오토마타

계산은 덧셈 몇 번인데, 데이터는 매번 메모리 끝에서 끌어와야 한다. 스텐실 최적화는 계산 최적화가 아니라 이삿짐 최적화다.

스텐실 연산(stencil computation)은 규칙 격자의 모든 점에서, 그 점을 중심으로 한 고정된 모양의 이웃 값들을 정해진 계수로 가중합해 새 값을 만드는 연산이다. 3차원 7점 라플라시안이라면

$$u^{\text{new}}_{i,j,k}=c_0\,u_{i,j,k}+c_1\big(u_{i\pm1,j,k}+u_{i,j\pm1,k}+u_{i,j,k\pm1}\big)$$

를 격자 전체에 적용한다. 유한차분법의 거의 모든 공간 연산자, 야코비 반복과 다중격자법의 평활기, 격자 볼츠만 방법의 스트리밍, 영상 블러·에지 검출, 셀룰러 오토마타가 모두 스텐실이다. 지진파 역산과 기상 모델의 실행 시간 대부분이 여기서 쓰인다.

수학적으로는 희소 행렬-벡터 곱의 특수한 경우이지만, 행렬을 저장하지 않고 패턴과 계수만 코드에 박는다는 점이 다르다. 그래서 인덱스 배열을 읽을 필요가 없고, 성능은 거의 전적으로 “격자 값을 메모리에서 얼마나 효율적으로 가져오느냐”로 정해진다. 이 문서는 그 데이터 이동을 줄이는 기법들을 다룬다. 성능 천장 계산 자체는 루프라인 모델 문서에 있다.

2. 왜 메모리 바운드인가

3D 7점 야코비 갱신 한 점(lattice update, LUP)에 드는 부동소수점 연산은 덧셈 67 번과 곱셈 12 번, 대략 8 flop 이다. 데이터는 이상적으로 원본 배열에서 한 번 읽고(8 B), 결과 배열에 한 번 쓰는데(8 B), 대부분의 CPU 캐시는 쓰기 전에 그 캐시 라인을 먼저 읽어 오므로(쓰기 할당) 8 B 가 더 든다. 합계 24 B/LUP, 산술강도 $8/24\approx0.33$ flop/B 다.

메모리 대역폭이 200 GB/s 인 소켓이라면 상한은

$$\frac{200\ \text{GB/s}}{24\ \text{B/LUP}}\approx 8.3\ \text{GLUP/s}\approx 67\ \text{GFLOP/s}$$

이고, 같은 소켓의 배정밀도 피크가 수 TFLOP/s 라면 피크의 몇 퍼센트에 불과하다. 스텐실 코드에서 “FLOP 효율 5%“는 실패가 아니라 정상이다. 진짜 질문은 대역폭 천장의 몇 %를 쓰고 있느냐다. 쓰기 전용 스트리밍 저장(non-temporal store)으로 쓰기 할당을 없애면 16 B/LUP 로 줄어 상한이 1.5 배가 된다.

문제는 이 24 B 가 이상적인 값이라는 것이다. 이웃 $u_{i,j,k\pm1}$ 은 한 평면 위아래에 있어서, 그 평면이 캐시에 남아 있지 않으면 같은 값을 메모리에서 세 번 읽게 된다.

3. 레이어 조건과 공간 블로킹

$k$ 바깥, $j$ 중간, $i$ 안쪽 순서로 훑는다고 하자. 점 $(i,j,k)$ 를 갱신할 때 필요한 $k-1$ 평면은 한 평면 전에, $k+1$ 평면은 지금 처음 읽힌다. 한 평면 뒤에 $k+1$ 평면이 다시 쓰이려면 원본 배열의 연속 세 평면이 캐시에 들어가야 한다. 이것을 레이어 조건(layer condition)이라 한다.

$$3\,N_iN_j\times8\ \text{B}\ \lesssim\ \tfrac12\,C_{\text{cache}}$$

(우변의 ½ 은 결과 배열과 충돌을 고려한 경험적 안전 계수다.) $N_i=N_j=512$ 면 좌변이 6.3 MB 라서, 마지막 단계 캐시(L3)에는 들어가도 코어별 L2(1~2 MB)에는 안 들어간다. 레이어 조건이 깨지고 행 단위 재사용만 남으면 원본을 평면 셋에서 따로 읽어 40 B/LUP 근처로 뛴다 — 같은 코드가 격자 크기만 바꿨는데 40% 느려지는 이유다.

해결책은 공간 블로킹이다. $j$ 방향을 폭 $b_j$ 의 띠로 잘라, 띠마다 $k$ 를 끝까지 훑는다. 그러면 조건이 $3N_ib_j\times8$ B 로 줄어 $b_j=64$ 면 약 790 KB 가 되어 L2 에 들어간다. 블록 크기는 캐시 크기·스레드 수·하드웨어 프리페처와 얽혀 있어 해석적으로 딱 떨어지지 않으므로, 몇 후보를 재 보고 고르는 자동 튜닝이 표준 관행이다(데이타 등 2008). 레이어 조건 분석을 캐시 계층마다 더 정밀하게 하려면 ECM 모델을 쓴다 — 루프라인이 “메모리 하나”만 보는 반면, ECM 은 L1↔L2↔L3↔메모리 각 단계 전송 시간을 더해 단일 코어 성능과 포화 코어 수까지 예측한다.

4. 시간 블로킹 — 천장을 뚫는 유일한 방법

공간 블로킹은 이상적인 24 B 로 돌아오게 할 뿐, 그 아래로는 못 내려간다. 한 시간 스텝마다 격자 전체를 한 번은 메모리에서 읽고 써야 하기 때문이다. 그 벽을 넘으려면 한 블록을 캐시에 올린 김에 여러 시간 스텝을 진행해야 한다. 이게 시간 블로킹이다.

문제는 의존성이다. 스텝 $t+1$ 의 점은 스텝 $t$ 의 이웃을 필요로 하므로, 블록 경계에서는 매 스텝 유효 영역이 스텐실 반경만큼 줄어든다. 처리 방식에 따라

• 중첩 타일링: 경계 영역을 겹치게 잡아 중복 계산한다. 단순하지만 스텝 수가 늘면 중복이 커진다.
• 웨이브프론트: 시간-공간 평면에서 비스듬한 띠를 따라 진행해, 여러 스레드가 서로 다른 시간 스텝을 파이프라인처럼 맡는다.
• 다이아몬드·사다리꼴 타일링: 시공간 타일을 다이아몬드 모양으로 잘라 중복 없이 병렬성을 확보한다. 다면체(polyhedral) 컴파일러가 자동으로 찾는 모양이기도 하다.

$T$ 스텝을 묶으면 메모리 트래픽이 이론상 $1/T$ 로 줄지만, 대가로 코드가 복잡해지고 경계 조건·소스 항·출력 시점처럼 스텝마다 끼어드는 일과 충돌한다. 그래서 실제 응용 코드에서는 생각보다 드물고, 순수 커널 벤치마크나 도메인 특화 언어 안에서 주로 쓰인다.

5. 병렬화 · GPU · DSL

분산 메모리. 격자를 영역 분할법으로 나누면 각 프로세스는 경계 바깥에 스텐실 반경만큼의 **고스트 셀(할로)**을 두고, 매 스텝 이웃과 MPI 로 교환한다. 교환량은 표면적, 계산량은 부피에 비례하므로 블록이 클수록 유리하다. 통신과 내부 점 계산을 겹치는 것(경계 먼저 계산 → 비동기 송신 → 내부 계산 → 수신 대기)이 표준 기법이다.

공유 메모리. OpenMP 로 $k$ 또는 블록 루프를 나누면 간단하지만, 대역폭 포화로 소켓의 일부 코어만으로 성능이 멈춘다. 초기화 루프도 같은 방식으로 병렬화해야 NUMA 배치가 맞는다.

벡터화. 안쪽 $i$ 루프는 연속 메모리라 SIMD 에 잘 맞지만, $u_{i\pm1}$ 접근이 정렬되지 않은 로드를 만든다. 정렬·별칭 정보를 컴파일러에 줘야 자동 벡터화가 되고, 고차 스텐실에서는 레지스터 내 셔플로 이웃을 만드는 기법도 쓴다.

GPU. 스레드 블록이 $(i,j)$ 타일을 맡아 공유 메모리에 올리고, $k$ 방향으로는 각 스레드가 위·현재·아래 값을 레지스터에 들고 미끄러지며 진행하는 2.5 차원 블로킹이 고전이다(미치케비치우스 2009). 지진파 모델링의 8차 정확도(25점) 스텐실이 이 방식의 대표 사례다.

도메인 특화 언어. 최적화가 하드웨어마다 달라지니, “무엇을 계산할지”와 “어떻게 배치할지”를 분리하는 언어들이 나왔다. 영상처리의 Halide 는 알고리즘과 스케줄(타일·벡터화·병렬화)을 따로 적고, 지구물리의 Devito 는 기호 PDE 에서 최적화된 C/GPU 스텐실 코드를 생성한다. 사람이 손으로 짠 커널과 비슷한 성능을 내면서 이식성을 얻는 것이 목표다.¹

6. 실무 감각

• B/LUP 로 생각하라. 측정한 성능(LUP/s)에 이론 B/LUP 를 곱해 STREAM 대역폭과 비교하면, 지금 코드가 레이어 조건을 지키는지 바로 보인다. 측정 바이트가 24 보다 한참 크면 블로킹부터 한다.²
• 배열 패딩. 격자 폭이 2 의 거듭제곱이면 이웃 평면들이 같은 캐시 세트로 몰려 충돌 미스가 난다. 512 대신 512+8 같은 패딩이 공짜 성능 향상을 주곤 한다.
• 연산 융합. 스텐실 하나 하고 배열 저장, 다음 스텐실 하고 저장… 대신 한 루프에서 여러 연산을 연달아 하면 중간 배열 왕복이 사라진다. 메모리 바운드에서는 계산을 더 해서라도 트래픽을 줄이는 게 이득이다.
• 정밀도. 단정밀도로 바꾸면 바이트가 절반이라 성능이 거의 2 배다. 정확도가 허락하는지는 혼합 정밀도 판단의 몫이다.
• 가우스-자이델은 다르다. 갱신한 값을 바로 쓰는 제자리 스텐실은 의존성 때문에 단순 병렬화가 안 된다. 적흑(red-black) 순서로 바꾸면 병렬화되지만 메모리 접근 효율은 떨어진다.³

7. 관련 문서

• 루프라인 모델 · ECM 모델 · 메모리 계층
• OpenMP · MPI · SIMD · 자동 벡터화 · GPU 컴퓨팅
• 유한차분법 · 다중격자법 · 격자 볼츠만 방법 · 반복법
• 영역 분할법 · 병렬 컴퓨팅 · 혼합 정밀도

8. Footnotes

1. DSL 이 성공하는 조건은 “사용자가 스케줄을 몰라도 된다”가 아니라 “스케줄을 바꿔도 알고리즘을 안 건드려도 된다”다. 성능 엔지니어와 물리학자가 같은 파일을 놓고 싸우지 않게 해 준다는 게 진짜 기능. ↩

2. 스텐실 튜닝 경험담의 국룰: 처음 재면 대역폭의 30%, 레이어 조건 맞추면 70%, 스트리밍 저장 넣으면 85%, 그 이상은 시간 블로킹인데 거기서부터는 코드가 아무도 못 읽게 된다. ↩

3. 반복법 교과서에서는 가우스-자이델이 야코비보다 두 배 빨리 수렴한다고 배우지만, 현대 하드웨어에서는 야코비가 병렬화·벡터화가 잘 돼서 벽시계 시간으로 이기는 경우가 많다. 수렴률과 실행 시간은 다른 이야기다. ↩