ScaLAPACK

1. 개요

ScaLAPACK Scalable LAPACK
대상	분산 메모리(MPI) 위의 밀집·띠 행렬 선형대수
계층	BLACS(통신) → PBLAS(병렬 BLAS) → ScaLAPACK(드라이버)
데이터 분배	2차원 블록 순환(block-cyclic)
이름 규약	LAPACK 이름 앞에 P: PDGESV, PZHEEVD
언어	Fortran 77/90 + C (BLACS), 열 우선 저장
후계·대안	SLATE, ELPA(고유값), Elemental, DPLASMA, cuSOLVERMp

행렬이 노드 한 대 메모리에 안 들어간다. 그럼 여러 대에 나눠 담고, 다 같이 갈아넣는다.

ScaLAPACK(Scalable LAPACK)은 LAPACK의 밀집 행렬 알고리즘(연립방정식·최소제곱·고유값·특이값)을 분산 메모리 병렬 컴퓨터로 옮긴 라이브러리로, 행렬을 2차원 프로세스 격자에 블록 순환 방식으로 흩어 놓고 각 프로세스가 자기 조각에 로컬 BLAS를 돌리며 MPI로 통신한다. 1990년대 초 테네시 대학·오크리지 국립연구소·UC 버클리 등이 개발했고, 1997년 SIAM 사용자 안내서로 정리됐다.

LAPACK이 “공유 메모리 한 대에서 캐시를 잘 쓰는 법”이었다면 ScaLAPACK은 “통신을 덜 하면서 모든 프로세스가 끝까지 놀지 않게 하는 법”이다. 30년 가까이 지난 지금도 밀도범함수이론 코드(Quantum ESPRESSO·VASP·CP2K)의 부분공간 대각화, 희소 직접 솔버 MUMPS의 루트 노드 분해, 경계요소법의 밀집 시스템이 이 위에서 돈다. 화려하진 않지만 HPC 소프트웨어 스택 바닥에 깔린 지반이다.

2. 세 층 구조

• BLACS(Basic Linear Algebra Communication Subprograms): 프로세스를 $P_r\times P_c$ 격자로 묶는 문맥(context)을 만들고, 행·열 방향 방송과 점대점 송수신을 행렬 블록 단위로 제공한다. 오늘날은 MPI 위의 얇은 껍데기다.
• PBLAS(Parallel BLAS): 분산된 행렬에 대한 BLAS 레벨 1~3. PDGEMM 이 분산 행렬곱이다. 내부적으로 각 프로세스의 로컬 BLAS 호출 + BLACS 통신으로 짜여 있다.
• ScaLAPACK 본체: LAPACK과 같은 이름 체계의 드라이버·계산 루틴. PDGETRF(LU), PDPOTRF(촐레스키 분해), PDGEQRF(QR 분해), PDSYEV·PDSYEVD·PDSYEVX·PDSYEVR(대칭 고유값), PDGESVD(특이값 분해) 등.

LAPACK 이름 읽는 법(D=배정밀도 실수, GE=일반, TRF=분해 …)은 그대로 승계되므로 LAPACK을 알면 ScaLAPACK 루틴 목록은 이미 반쯤 읽은 셈이다.

3. 2차원 블록 순환 분배 — 이 라이브러리의 심장

행렬을 $MB\times NB$ 크기 블록으로 자르고, 블록 행 $I_b$ 를 프로세스 행 $I_b \bmod P_r$ 에, 블록 열 $J_b$ 를 프로세스 열 $J_b\bmod P_c$ 에 돌려 가며 배정한다. 0부터 세는 전역 행 번호 $i$ 의 주인과 로컬 위치는

$$p_r = \Big\lfloor \tfrac{i}{MB}\Big\rfloor \bmod P_r,\qquad i_{\text{loc}} = \Big\lfloor \tfrac{i}{MB\,P_r}\Big\rfloor MB + (i \bmod MB)$$

(시작 프로세스 RSRC=0 일 때). 각 프로세스가 가진 로컬 행 수는 NUMROC 함수가 계산한다. $N=1000$, $NB=64$ 에서 직접 세 보면 블록은 64짜리 15개와 40짜리 1개이고,

프로세스 수	로컬 행 수
3	360, 320, 320
4	256, 256, 256, 232

로 나뉜다. 마지막 자투리 블록 하나만큼만 불균형이 생긴다.

왜 순환인가. LU·촐레스키·QR은 왼쪽 위에서 오른쪽 아래로 진행하며 남은 부분행렬이 계속 줄어든다. 행렬을 큰 덩어리로 나눠 주면(블록 분배) 앞쪽 덩어리를 가진 프로세스는 초반에 일을 끝내고 남은 시간 내내 논다. 순환 분배는 각 프로세스가 행렬 전역에 고루 흩어진 조각을 갖게 해 끝까지 일거리가 남게 만든다.

왜 블록인가. 블록이 1이면(순수 순환) 부하 균형은 최고지만 로컬 연산이 BLAS-1·2 로 쪼그라들고 메시지가 잘아진다. 블록을 키우면 BLAS-3 효율과 메시지 수는 좋아지지만 부하 불균형이 늘어난다.

간단한 1차원 모형으로 이 줄다리기를 재 봤다. $n=1024$ 열, 프로세스 8개, 우측 갱신 작업만 셈하고 매 단계 가장 바쁜 프로세스 시간을 합산한 병렬 효율(이상적 = 1):

분배	블록 크기	효율
블록 순환	1	0.995
블록 순환	8	0.961
블록 순환	32	0.867
블록 순환	64	0.777
순수 블록 (프로세스당 128열)	—	0.670

통신 비용과 BLAS 효율을 뺀 모형이라 블록 1이 이기지만, 실제로는 그 둘이 작은 블록을 응징한다. 결국 현실의 최적은 “블록당 BLAS-3가 제 속도를 내는 최소 크기” 이고, 실무 기본값이 32~256 사이(흔히 64)에 몰리는 이유다. 위 표가 알려 주는 진짜 교훈은 $n/(NB\cdot P)$, 즉 각 프로세스가 받는 블록 순환 횟수가 너무 적으면 순환의 이득이 사라진다는 것이다. $NB=64$ 면 각 프로세스가 고작 2바퀴분만 받는다.

4. 격자 모양과 통신량

$P$ 개 프로세스를 $P_r\times P_c$ 로 나누는 방식도 성능을 좌우한다. 2차원 분배에서 LU 한 번의 프로세스당 통신량은 대략 $O(N^2/\sqrt{P})$ 단어로, 1차원(행 또는 열만) 분배의 $O(N^2)$ 보다 $\sqrt{P}$ 배 적다. 이게 2차원 격자를 쓰는 근본 이유다. 경험칙은 정사각형에 가깝게, LU처럼 열 방향 피벗 탐색이 잦은 연산은 열 쪽을 약간 많게 잡는 것이다. 64 프로세스면 8×8, 128이면 8×16. 소수 개수의 프로세스(예: 7)를 주면 1×7 격자밖에 안 돼 1차원으로 퇴화하니, 작업 배치 때 프로세스 수를 합성수로 맞추는 것이 국룰이다.¹

메모리는 프로세스당 $N^2/P$ 원소라 프로세스를 늘리는 만큼 더 큰 행렬을 담는다. 이 덕분에 ScaLAPACK은 약한 확장성(문제 크기를 프로세스 수에 비례해 키우기)에서 특히 강하다. 반대로 고정 크기 문제를 프로세스만 늘리면(강한 확장성) 로컬 블록이 작아져 결국 통신이 계산을 덮는다.

5. 쓰는 법 — 그리고 물어뜯는 지점

전형적인 호출 순서:

1. BLACS_GET / BLACS_GRIDINIT 으로 $P_r\times P_c$ 격자 문맥 생성
2. NUMROC 로 로컬 행·열 수 계산, 로컬 배열 할당
3. DESCINIT 으로 배열 기술자(descriptor) 9개 정수 작성 — 행렬 종류, 문맥, 전역 크기 M·N, 블록 MB·NB, 시작 프로세스 RSRC·CSRC, 로컬 선행 차원 LLD
4. 로컬 조각 채우기 → PDGESV 등 호출 → BLACS_GRIDEXIT

현장에서 자주 터지는 곳:

• 열 우선 저장. Fortran 규약이다. C/C++ 에서 행 우선 배열을 그대로 넘기면 전치를 푼 셈이 된다. LAPACK 문서의 같은 함정이 분산판에서 더 고약하게 재현된다 — 블록 경계까지 꼬인다.
• LLD 와 로컬 크기. 기술자의 LLD 를 전역 크기로 적거나 NUMROC 결과보다 작게 주면 조용히 메모리를 밟는다. LLD ≥ max(1, 로컬 행 수).
• 작업공간 질의. LWORK=-1 로 한 번 불러 필요한 작업공간 크기를 받은 뒤 할당하는 LAPACK 관례가 그대로다. 고유값 루틴은 프로세스마다 필요량이 달라 최댓값을 써야 한다.
• 입출력 병목. 행렬 전체를 랭크 0에서 읽어 흩뿌리는 코드는 $N$ 이 커지면 랭크 0 메모리가 먼저 터진다. 각 프로세스가 자기 블록만 읽게 하거나 PDGEMR2D(서로 다른 분배 사이 재분배)를 쓴다.
• 고유값 루틴 선택. 고유벡터가 다 필요하면 분할정복 PDSYEVD 나 MRRR 기반 PDSYEVR 이 PDSYEV(QR 반복)보다 대개 훨씬 빠르다. 일부만 필요하면 PDSYEVX·PDSYEVR 의 범위 지정. 이분법+역반복의 PDSYEVX 는 고유값이 뭉쳐 있으면 재직교화 때문에 느려지고 직교성이 나빠질 수 있다. 원리는 분할 정복 고유값 알고리즘·역반복법 문서 참고.

6. 생태계 — 노병은 죽지 않는다

ScaLAPACK의 약점은 설계 시대에서 온다. 프로세스 단위 동기적 진행(fork-join 스타일)이라 코어가 수십 개인 노드 안에서 스레드를 잘 못 살리고, GPU 컴퓨팅을 고려하지 않았다. 그래서 후계자들이 나왔다.

• ELPA: 독일 계열이 만든 대칭/에르미트 고유값 전용 라이브러리. 2단계 삼중대각화로 ScaLAPACK보다 빠르고 GPU를 지원해, DFT 코드들이 대각화 백엔드로 많이 갈아탔다. 기술자·분배는 ScaLAPACK 것을 그대로 받는다.
• SLATE: 미국 엑사스케일 컴퓨팅 프로젝트(ECP)가 “ScaLAPACK 대체”를 명시 목표로 만든 C++ 라이브러리. 타일 기반 작업 스케줄링과 GPU 가속, ScaLAPACK 호환 API 를 제공한다.
• Elemental, DPLASMA: 원소 순환 분배·작업 그래프 런타임 같은 다른 설계를 택한 분산 밀집 선형대수.
• cuSOLVERMp 같은 벤더 라이브러리: 다중 GPU·다중 노드용.

그럼에도 ScaLAPACK이 남는 이유는 인터페이스다. 블록 순환 분배 + 9칸 기술자는 이 바닥의 공용어가 됐고, 후계자들도 그 형식으로 행렬을 받는다. 2010년대 이후에도 2.x 판이 나오며 MRRR 대칭 고유값 풀이기와 개선된 비대칭 고유값 루틴이 추가됐다. 새 코드를 짤 때는 SLATE·ELPA 를 먼저 보되, ScaLAPACK 기술자 형식은 알아 둬야 남의 코드를 읽을 수 있다.²

7. 관련 문서

• LAPACK · BLAS · MPI · 병렬 컴퓨팅 · GPU 컴퓨팅
• LU 분해 · 촐레스키 분해 · QR 분해 · 특이값 분해 · 고유값 문제
• 분할 정복 고유값 알고리즘 · 역반복법 · 부하 분산
• MUMPS · 희소 직접 솔버 · PETSc
• 밀도범함수이론 · Quantum ESPRESSO · VASP · ELPA · SLATE

8. Footnotes

1. 클러스터 대기열에 프로세스 97개를 던졌다가 97×1 격자로 하루를 날리는 일이 실제로 있다. 소수는 정수론에서는 아름답지만 프로세스 격자에서는 킹받는 존재다. ↩

2. 사용자 안내서 표지의 저자 목록이 열세 명이다. 이름만 봐도 1990년대 미국 수치선형대수 학계의 절반이 모였다는 걸 알 수 있다. 동가라·뎀멜·해머링 같은 이름이 LAPACK 표지에도 나온다는 건 덤. ↩