BLAS

1. 개요

BLAS Basic Linear Algebra Subprograms
정체	밀집 선형대수 기본 연산의 인터페이스 표준
레벨 1	벡터-벡터 (1979, 로슨 외)
레벨 2	행렬-벡터 (1988, 동가라 외)
레벨 3	행렬-행렬 (1990, 동가라 외)
대표 구현	OpenBLAS, oneMKL, BLIS, Accelerate, cuBLAS
기본 배치	열 우선(column-major) + 선행 차원 LDA
얹혀 사는 것	LAPACK, ScaLAPACK, MUMPS, NumPy/SciPy, MATLAB

dgemm 한 줄이 이 세상 부동소수점 연산의 꽤 큰 몫을 태운다. 딥러닝도 결국 곱셈이다.

BLAS(Basic Linear Algebra Subprograms)는 벡터 덧셈·내적, 행렬-벡터 곱, 행렬-행렬 곱처럼 밀집 선형대수의 가장 밑바닥 연산을 정해진 이름과 인자 순서로 호출하게 만든 인터페이스 표준이다. 1979년 로슨·핸슨·킨케이드·크로그가 포트란용 레벨 1을 내놓았고, 1988년 레벨 2, 1990년 레벨 3이 동가라 등에 의해 추가됐다.

핵심은 BLAS가 라이브러리가 아니라 계약이라는 점이다. Netlib의 참조 구현은 “이렇게 동작하면 된다”는 정답지일 뿐 느리고, 실제 성능은 각 CPU·GPU 제조사와 오픈소스 프로젝트가 같은 이름의 함수를 기계에 맞춰 갈아 넣은 구현에서 나온다. 그래서 LAPACK·ScaLAPACK·MUMPS 같은 상위 라이브러리는 BLAS만 부르도록 짜 두고, 링크할 때 OpenBLAS를 붙이느냐 MKL을 붙이느냐로 속도가 몇 배씩 갈린다. 레벨 구분과 LAPACK의 블록 알고리즘은 LAPACK 문서에 정리돼 있으니, 여기서는 “왜 레벨 3만 빠른가”를 숫자로 확인하고, 그 빠른 dgemm이 실제로 어떻게 만들어지는지에 집중한다.

2. 이름 읽는 법과 호출 규약

BLAS 함수명은 LAPACK과 같은 문법을 쓴다. 앞 글자가 자료형, 가운데가 행렬 구조, 끝이 연산이다.

조각	뜻	예
첫 글자	s 단정밀도 실수, d 배정밀도 실수, c·z 복소	dgemm, zgemm
구조	ge 일반, sy 대칭, he 에르미트, tr 삼각, gb 띠	dsymm, dtrsm
연산	mm 행렬곱, mv 행렬-벡터, sv 삼각 풀이, rk 랭크-$k$ 갱신	dsyrk, dtrsv

가장 유명한 dgemm은 $C \leftarrow \alpha\,\mathrm{op}(A)\,\mathrm{op}(B) + \beta C$ 를 계산한다. $\mathrm{op}$ 는 전치 여부를 TRANSA·TRANSB 문자로 고른다. 레벨 3의 dtrsm은 삼각행렬 여러 우변 풀이, dsyrk는 $C\leftarrow \alpha AA^T+\beta C$ 로 촐레스키 분해의 몸통을 이룬다.

호출 규약에서 사람을 물어뜯는 지점은 정해져 있다.

• 열 우선 배치. 포트란 출신이라 $A_{ij}$ 가 A[i + j*LDA] 에 있다. C·NumPy의 행 우선 배열을 그대로 넘기면 전치된 행렬을 계산한다. 국룰 우회법은 $C^T = B^T A^T$ 로 인자 순서를 바꿔 넘기는 것이고, CBLAS는 CblasRowMajor 플래그로 이걸 대신 해 준다.
• 선행 차원 LDA. 큰 배열의 부분 블록을 복사 없이 넘기려고 “열 사이 간격”을 따로 받는다. 블록 알고리즘이 패널을 잘라 넘길 수 있는 이유가 이것이다. LDA 를 행 수로 착각해 블록을 넘기면 조용히 엉뚱한 원소를 읽는다.
• 정수 폭. LP64(32비트 정수)와 ILP64(64비트 정수) 빌드가 따로 있다. 원소가 $2^{31}$ 을 넘는 행렬을 LP64에 넘기면 크기 인자가 넘친다. 반대로 두 빌드를 한 프로세스에 섞어 링크하면 심볼이 꼬인다.¹

3. 왜 레벨 3만 빠른가 — 연산 강도 실측

연산 하나가 메모리에서 읽어 오는 바이트당 몇 번의 부동소수점 연산을 하는지를 연산 강도라 한다. 배정밀도에서 대략

연산	flop	읽는 원소	flop/원소
daxpy ($y\leftarrow \alpha x+y$)	$2n$	$3n$	$2/3$
dgemv ($y\leftarrow Ax$)	$2n^2$	$n^2$	$2$
dgemm ($C\leftarrow AB$)	$2n^3$	$3n^2$ (이상적)	$2n/3$

레벨 1·2는 $n$ 을 키워도 원소당 연산 수가 상수라 메모리 대역폭에 묶이고, 레벨 3만 $n$ 에 비례해 커져서 연산 유닛을 포화시킬 수 있다. 루프라인 모델로 말하면 레벨 1·2는 대역폭 경사면 위에, 레벨 3은 지붕 아래에 산다.

이 문서를 쓰며 4코어 Xeon 컨테이너(OpenBLAS 0.3.31, 배정밀도)에서 직접 재 봤다. 수치는 기계·부하에 따라 흔들리니 비율만 보라.

방식	크기	GFLOP/s
C 삼중 루프 ijk (-O3 -march=native)	$1024^2$	0.86
C 삼중 루프 ikj (안쪽 루프 연속 접근)	$1024^2$	5.7
OpenBLAS dgemm, 1스레드	$1024^2$	37
OpenBLAS dgemm, 4스레드	$2048^2$	약 170
OpenBLAS dgemv, 4스레드	$4096^2$	19

교과서 그대로의 ijk 루프와 튜닝된 BLAS 사이가 약 40배(단일 스레드)다. 루프 순서만 ikj 로 바꿔 $B$ 를 행 방향으로 연속 접근하게 해도 7배가 나오는데, 그래도 BLAS의 6분의 1이다. 같은 4스레드에서 dgemv 가 dgemm 의 약 9분의 1인 것은 위 표의 연산 강도 차이가 그대로 드러난 것이다. 행렬 4096×4096 하나가 128MB라 매 호출마다 메모리에서 통째로 끌어와야 하고, 원소 하나당 곱셈-덧셈 한 번밖에 못 한다.

4. 빠른 dgemm 은 이렇게 만든다 — 고토의 해부도

2000년대 중반 고토 가즈시게의 GotoBLAS가 당시 벤더 라이브러리를 이기며 정리된 설계가 지금 OpenBLAS·BLIS의 뼈대다.² 고토와 판 더 헤이인(2008)의 요지는 **“캐시 계층마다 한 행렬씩 붙박아 두고, 레지스터에서 작은 외적을 반복한다”**이다.

1. 바깥 두 루프 — $B$ 패널을 L3에. $B$ 를 $k_c\times n_c$ 패널로 자르고, 연속 메모리로 패킹(복사·재배치)한다. 이 패널이 L3 캐시에 머문다.
2. 가운데 루프 — $A$ 블록을 L2에. $A$ 를 $m_c\times k_c$ 블록으로 잘라 역시 패킹한다. L2에 들어가는 크기로 고른다.
3. 마이크로커널 — 레지스터에서 외적. $C$ 의 $m_r\times n_r$ 조각(예: BLIS 의 AVX2 배정밀도 커널은 $6\times 8$)을 레지스터에 통째로 들고, $A$ 의 $m_r$ 길이 열 조각과 $B$ 의 $n_r$ 길이 행 조각의 외적을 $k_c$ 번 누적한다. 한 단계마다 $m_r+n_r$ 개를 읽고 $m_r n_r$ 번 FMA를 한다.

패킹이 의외의 주인공이다. 복사에 $O(n^2)$ 이 들지만 $O(n^3)$ 연산 내내 캐시 라인과 TLB 미스 없이 순차 접근하게 해 주니 남는 장사다. 그리고 마이크로커널 하나만 기계마다 어셈블리·인트린직으로 짜고 나머지 루프는 공통 C 코드로 두면 된다는 게 BLIS의 발견이다. 새 아키텍처에 BLAS 전체(수백 개 함수)를 다시 짤 필요 없이 커널 몇 개만 갈아 끼우면 레벨 3 전부가 따라온다. 그 이전 세대인 ATLAS(휘틀리·동가라, 1998)는 블록 크기를 설치 시점에 자동 탐색(오토튜닝)으로 정했고, 지금도 그 발상은 GPU 커널 튜닝에 살아 있다.

레벨 3 연산 대부분은 dgemm 으로 환원된다. dtrsm 은 작은 대각 블록을 풀고 나머지를 dgemm 으로 갱신하고, dsyrk 는 대칭성 때문에 반만 계산할 뿐 구조가 같다. 그래서 BLAS 구현의 품질 = dgemm 마이크로커널의 품질이라고 해도 크게 틀리지 않다.

5. 구현 생태계와 고르는 법

구현	특징
Netlib 참조 BLAS	정답지. 단순 포트란 루프라 느리다. 결과 검증용
OpenBLAS	GotoBLAS2 계승. 오픈소스 기본값, NumPy 휠에 기본 동봉
Intel oneMKL	인텔 CPU에서 최상급. 희소·FFT까지 포함한 종합 패키지
BLIS / AMD AOCL-BLIS	위 해부도를 프레임워크로 만든 것. AMD의 공식 BLAS
Apple Accelerate	애플 실리콘의 행렬 보조 유닛까지 쓴다
cuBLAS / rocBLAS	GPU판. 텐서 코어로 저정밀 GEMM을 수십~수백 TFLOP/s로

GPU 쪽은 같은 이름 규약을 쓰되 메모리가 장치에 있어야 한다. 작은 행렬 수천 개를 따로따로 부르면 커널 실행 오버헤드에 묻히므로 배치 GEMM(gemmBatched, gemmStridedBatched) 같은 확장이 생겼다. 저정밀 GEMM을 반복 정제로 배정밀도까지 끌어올리는 기법은 혼합 정밀도 문서에 있다.

6. 실전 지뢰

• 스레드 과잉 구독. OpenBLAS가 코어 수만큼 스레드를 띄우는데, 바깥에서 이미 MPI 랭크나 OpenMP 스레드로 코어를 다 쓰고 있으면 스레드가 코어 수의 제곱만큼 생겨 오히려 느려진다. MPI 코드에서는 OPENBLAS_NUM_THREADS=1(또는 OMP_NUM_THREADS, MKL_NUM_THREADS)이 국룰이다.
• 비트 재현성. 스레드 수나 CPU 명령어 집합이 바뀌면 누적 순서가 바뀌고, 부동소수점 연산은 결합법칙이 안 성립하니 마지막 몇 비트가 달라진다. 회귀 테스트가 == 로 비교하고 있으면 “어제는 됐는데”가 터진다. MKL은 이를 위해 조건부 재현성(CNR) 모드를 따로 둔다.
• 작은 행렬. $n\lesssim 16$ 정도에서는 호출 오버헤드와 패킹 비용이 계산을 넘는다. 유한요소 요소 행렬처럼 작은 곱이 수백만 번 반복되면 BLAS 대신 컴파일 시점 크기가 고정된 전용 커널(Eigen 고정 크기, LIBXSMM 등)이 낫다.
• 희소 행렬은 BLAS가 아니다. BLAS는 밀집 전용이다. 희소행렬 곱·풀이는 별도 라이브러리(Sparse BLAS, SuiteSparse 등)의 몫이고, 희소 직접 솔버는 내부에서 밀집 프론트 행렬을 만든 뒤에야 BLAS 레벨 3을 부른다.
• 슈트라센은 안 쓴다. 슈트라센 알고리즘은 점근적으로 $O(n^{2.807})$ 이지만, 오차 한계가 원소별이 아닌 노름 단위로 약해지고 메모리 접근이 복잡해 주류 BLAS의 dgemm 은 고전적 $2n^3$ 알고리즘을 쓴다.³

7. 시뮬레이션에서의 의미

시뮬레이션 코드에서 BLAS를 직접 부르는 일은 생각보다 적다. 대신 간접적으로 거의 모든 곳에 있다. LU 분해·QR 분해·고유값 풀이는 LAPACK을 통해, 다중 프론트 희소 솔버는 프론트 행렬 처리로, 스펙트럴 방법의 미분 행렬 곱, 밀도범함수 계산의 부분공간 대각화, 신경망 대리 모델 학습은 GEMM 그 자체로 BLAS를 태운다. 그래서 “코드 하나도 안 고쳤는데 링크만 바꿨더니 3배 빨라졌다”는 일이 실제로 일어난다.⁴

반대로 성능 분석에서 첫 질문도 여기서 나온다. 프로파일러를 돌렸을 때 시간이 dgemm 에 몰려 있으면 남은 개선 여지는 “더 좋은 BLAS”나 “더 적은 곱”뿐이고, daxpy·ddot 같은 레벨 1에 몰려 있으면 알고리즘을 블록화해 레벨 3으로 끌어올리는 것이 먼저다. GMRES 같은 크릴로프 방법이 내적을 묶어 처리하는 블록·s-단계 변형으로 가는 것도 같은 이유다.

8. 관련 문서

• LAPACK · ScaLAPACK · MUMPS · SuiteSparse
• 부동소수점 연산 · 혼합 정밀도 · GPU 컴퓨팅 · 병렬 컴퓨팅
• LU 분해 · 촐레스키 분해 · QR 분해 · 희소행렬
• 루프라인 모델 · 슈트라센 알고리즘 · SIMD

9. Footnotes

1. NumPy가 OpenBLAS를 ILP64로 동봉하면서 심볼 이름 끝에 64_ 접미사를 붙이는 이유가 이것이다. 다른 확장 모듈이 LP64 OpenBLAS를 따로 들고 오면 한 프로세스에 BLAS 두 벌이 공존하는 기묘한 상황이 된다. 그리고 둘 다 코어 수만큼 스레드를 띄운다. ↩

2. 고토는 일본 출신 엔지니어로, 텍사스 대학 TACC에서 일하던 시절 손으로 짠 커널이 벤더 라이브러리를 이겨 이름을 날렸다. “TLB 미스를 줄여라”가 그의 트레이드마크 통찰인데, 당시 다들 캐시만 보고 있었다. ↩

3. 반대로 “무조건 고전 알고리즘이어야 한다”는 규정이 표준에 있는 것은 아니다. 그냥 다들 안 쓸 뿐이다. 슈트라센을 구현한 실험적 BLAS는 있고, 특정 크기에서 실제로 빨라지기도 한다. ↩

4. 그리고 반대 방향으로도 일어난다. 클러스터 모듈을 바꿨더니 참조 BLAS가 링크돼 있어서 하룻밤 돌 계산이 일주일이 되는 사태. ldd 로 어느 libblas 가 붙었는지 확인하는 습관은 현생을 지켜 준다. ↩