1. 개요[편집]
| LAPACK Linear Algebra PACKage | |
|---|---|
| 첫 공개 | 1992 (사용자 안내서 SIAM, 3판 1999) |
| 언어 | Fortran 90 (초기 F77) |
| 전신 | EISPACK(1970s) · LINPACK(1979) |
| 대상 | 조밀·밴드 행렬 (희소는 대상 아님) |
| 기반 | BLAS 레벨 3 |
| 저장 방식 | 열 우선(column-major) + leading dimension |
| 라이선스 | 수정 BSD |
numpy.linalg.solve(A, b)한 줄 뒤에는 1992년에 이름이 정해진 포트란 서브루틴이 있다.
LAPACK(Linear Algebra PACKage)은 조밀 및 밴드 행렬에 대한 연립방정식 풀이, 최소제곱, 고유값 문제, 특이값 분해의 사실상 표준 구현체다. LU 분해·QR 분해·촐레스키 분해·특이값 분해·고유값 문제를 다루는 모든 상용·오픈소스 코드가 결국 이 API로 내려간다. MATLAB의 백슬래시, NumPy/SciPy의 linalg, R의 solve, 상용 유한요소법 코드의 조밀 블록 처리가 전부 여기서 만난다.
이 문서는 알고리즘이 아니라 라이브러리로서의 LAPACK을 다룬다. 각 분해의 수학은 위 개별 문서에 있고, 여기서는 왜 이런 구조로 짜여 있으며 사용자가 어디서 넘어지는지를 본다. 결론부터 말하면 LAPACK을 이해한다는 것은 대부분 성능이 어디서 나오는가(BLAS 레벨 3) 와 호출 규약이 어디서 배신하는가(열 우선·lda·info) 두 가지를 아는 일이다.
2. BLAS 레벨과 블록 알고리즘[편집]
LAPACK의 성능은 사실상 BLAS(Basic Linear Algebra Subprograms) 위에 올라탄 결과다. BLAS는 세 레벨로 나뉜다.
| 레벨 | 연산 | 연산량 | 데이터 | 대표 |
|---|---|---|---|---|
| 1 | 벡터-벡터 | daxpy, ddot | ||
| 2 | 행렬-벡터 | dgemv, dger | ||
| 3 | 행렬-행렬 | dgemm |
핵심은 마지막 줄이다. 레벨 3만이 데이터 한 바이트당 연산 수가 에 비례하므로, 캐시에 올린 블록을 여러 번 재사용해 부동소수점 유닛을 포화시킬 수 있다. 레벨 1·2는 아무리 잘 짜도 메모리 대역폭에서 막힌다.
족보를 짚으면 이렇다. 1970년대에 고유값 문제를 담당한 EISPACK과 1979년의 LINPACK이 먼저 있었고, 두 갈래를 하나의 규약으로 통합해 다시 쓴 것이 LAPACK이다. 전신인 LINPACK은 레벨 1 BLAS 기반이었고, 캐시 계층이 깊어진 1980년대 기계에서 성능이 무너졌다. LAPACK의 존재 이유가 여기 있다 — 같은 알고리즘을 블록 단위로 다시 짜서 작업량 대부분을 dgemm으로 밀어 넣는 것. 예컨대 dgetrf(LU)는 열 하나씩 소거하지 않고 열 패널을 처리한 뒤 나머지를 행렬-행렬 곱으로 갱신한다. 블록 크기 는 ILAENV가 기계별로 답한다.1
3. 세 계층과 이름 읽는 법[편집]
루틴은 세 계층으로 나뉜다.
- 드라이버 루틴 — 문제 하나를 통째로 푼다.
dgesv(연립방정식),dgels(최소제곱),dsyev(대칭 고유값),dgesvd(SVD). - 계산 루틴 — 한 단계만 한다.
dgetrf(LU 분해),dgetrs(삼각계 풀이),dgeqrf(QR 분해). - 보조 루틴 — 하우스홀더 반사자 생성
dlarfg, 기계 상수 조회dlamch등.
이름은 암호가 아니라 규약이다. XYYZZZ 꼴로, X는 정밀도(s 단정도, d 배정도, c/z 복소), YY는 행렬 종류(ge 일반, sy 대칭, po 대칭 양정부호, gb 일반 밴드, tr 삼각, he 에르미트), ZZZ는 계산 내용(sv 풀이, trf 분해, trs 대입, ev 고유값, qrf QR)이다.
dgesv= 배정도 + 일반 + 풀이dsyev= 배정도 + 대칭 + 고유값dpotrf= 배정도 + 양정부호 + 분해(= 촐레스키 분해)
드라이버에는 간이판과 확장판이 따로 있다. dgesv가 그냥 푸는 데 반해 dgesvx는 행·열 스케일링, 조건수 추정, 반복 개선, 전진·후진 오차 한계까지 돌려준다. 답만 필요한지 답의 신뢰도까지 필요한지에 따라 고르면 된다. 대칭 고유값 하나만 해도 dsyev(QR), dsyevd(분할 정복), dsyevr(MRRR), dsyevx(선택적 계산)가 다 있는데, 전부 필요하면 dsyevd, 일부 구간만 필요하면 dsyevr가 대체로 정답이다.
4. C·파이썬 사용자를 물어뜯는 지점[편집]
열 우선 저장. 포트란 관례를 그대로 쓴다. C의 2차원 배열은 행 우선이므로 그대로 넘기면 조용히 전치된 행렬을 푼다. 대칭 행렬로 시험하면 멀쩡히 통과했다가 비대칭에서 터지는 최악의 버그 패턴이 여기서 나온다.2
leading dimension(lda). 배열의 열 간 보폭이다. lda ≥ max(1, m)이면 되고, 큰 배열의 부분행렬을 복사 없이 넘길 수 있게 해 주는 장치다. 부분행렬을 넘기면서 lda를 부분행렬 크기로 준 순간 엉뚱한 메모리를 읽는다.
작업배열 조회. 다수의 루틴이 work 배열을 요구한다. lwork = -1로 한 번 부르면 실제 계산은 하지 않고 최적 크기를 work(1)에 실수로 돌려준다. 이걸 정수로 캐스팅해 할당한 뒤 진짜로 부르는 2단 호출이 관례다. 최소 크기만 주면 돌긴 하지만 언블록 경로로 떨어져 성능이 몇 배 손해다.
info 해석. info = 0이면 성공, info = -i면 번째 인자가 잘못됐다는 뜻이며, info > 0은 루틴마다 의미가 다르다. dgesv에서 info = i는 가 정확히 0이라 특이하다는 뜻이고, dsyev에서는 수렴하지 못한 비대각 원소 개수다. 양수 info를 무시하고 결과 배열을 그대로 쓰는 것이 실무에서 가장 흔한 사고다.
5. 생태계와 경계선[편집]
API는 하나지만 구현은 여럿이다. 넷립의 레퍼런스 LAPACK이 정의를 쥐고 있고, 실제 속도는 커널을 손으로 최적화한 구현이 낸다 — OpenBLAS, BLIS, 인텔 oneMKL, AMD AOCL, 애플 Accelerate. 분산 메모리는 2차원 블록 순환 분배를 쓰는 ScaLAPACK, 멀티코어 타일 알고리즘은 PLASMA, GPU 컴퓨팅 쪽은 MAGMA·cuSOLVER·rocSOLVER가 같은 이름 규약을 대체로 승계했다. 즉 인터페이스만 붙잡고 있으면 하드웨어를 갈아타도 코드가 산다.
경계선도 분명하다. LAPACK은 희소행렬을 다루지 않는다. 0을 전부 저장하고 채움(fill-in)을 그대로 받아들이는 설계라, 유한요소 강성행렬을 통째로 넘기면 메모리부터 터진다. 그쪽은 SuiteSparse(UMFPACK·CHOLMOD), MUMPS, SuperLU, PETSc의 영역이고, 희소 고유값은 ARPACK이다. 다만 이 도구들도 내부의 조밀 블록·초노드 연산에서는 결국 BLAS 3와 LAPACK을 부른다.
마지막으로 실무 결론 하나. 직접 짜지 마라. 가우스 소거법을 손으로 구현하는 것은 교육적으로 훌륭하지만, 30년 넘게 검증된 피벗 전략·스케일링·후진 안정성·예외값 처리를 재현할 수는 없다.3
6. 관련 문서[편집]
- LU 분해 · QR 분해 · 촐레스키 분해 · 특이값 분해
- 고유값 문제 · 최소자승법 · 조건수
- 희소행렬 · ARPACK · 반복법
- 병렬 컴퓨팅 · GPU 컴퓨팅 · 부동소수점 연산
- 가우스 소거법 · 기븐스 회전 · 강성행렬
7. Footnotes[편집]
-
그래서 “레퍼런스 LAPACK이 느리다”는 말은 절반만 맞다. 느린 건 대개 함께 링크된 레퍼런스 BLAS 쪽이고, 같은 LAPACK 소스에 OpenBLAS를 물리면 성능이 자릿수 단위로 달라진다. 파이썬 환경 두 개가 같은 코드에서 10배 차이 나는 미스터리는 거의 항상 여기서 나온다. ↩
-
C 표준 인터페이스 LAPACKE는
LAPACK_ROW_MAJOR를 받아 주지만, 내부적으로 전치 복사를 하거나 등가 문제로 바꿔 푸는 것이라 공짜가 아니다. 정공법은 애초에 열 우선으로 데이터를 잡는 것이고, NumPy에서order='F'를 쓰면 복사 한 번을 아낀다. 물론 그 한 번을 아끼려다 배열 레이아웃을 헷갈려 하루를 쓰는 게 국룰이다. ↩ -
유명한 반례로 LINPACK 벤치마크가 있다. TOP500 슈퍼컴퓨터 순위가 아직도 조밀 연립방정식 풀이 속도로 매겨지는 이유는, 이 문제가 실제 워크로드를 대표해서가 아니라 40년치 비교 가능한 숫자가 쌓여 있어서다. 그 관성이 싫어서 나온 것이 HPCG 벤치마크다. ↩