선형대수 연산의 핵심인 행렬 곱셈에서, 고성능 수학 라이브러리가 단순 반복문 구현 대비 얼마나 효율적인지 측정해보겠습니다. LAPACK의 DGEMM 루틴과 기본적인 삼중 루프 구현을 다양한 차원에서 비교합니다.
DGEMM 래퍼 설계
Fortran 인터페이스를 C++에서 직접 호출하기 위해 타입 안전한 래퍼를 구성합니다. Fortran의 열 우선 저장 방식과 C의 행 우선 저장 방식 차이를 고려하여 파라미터 순서를 조정합니다.
#include <iostream>
#include <vector>
#include <chrono>
#include <random>
// Fortran BLAS 인터페이스
extern "C" {
void dgemm_(const char* TRANSA, const char* TRANSB,
const int* M, const int* N, const int* K,
const double* ALPHA, const double* A, const int* LDA,
const double* B, const int* LDB,
const double* BETA, double* C, const int* LDC);
}
// LAPACK DGEMM 래퍼: C = A × B
void optimized_gemm(int dim, const double* matA, const double* matB, double* result) {
const char no_transpose = 'N';
const int n = dim;
const double scale_a = 1.0;
const double scale_c = 0.0;
// Fortran은 열 우선이므로, C에서 B^T × A^T로 전달하면
// 결과적으로 C++ 행렬 곱셈 형태가 됨
dgemm_(&no_transpose, &no_transpose,
&n, &n, &n,
&scale_a, matB, &n, matA, &n,
&scale_c, result, &n);
}기준 구현: 단순 중 루프
캐시 고려 없이 직접 작성한 참조 구현입니다. 행 우선 접근 패턴으로 메모리 지역성이 낮습니다.
void naive_multiply(int dim, const double* srcA, const double* srcB, double* dst) {
for (int row = 0; row < dim; ++row) {
for (int col = 0; col < dim; ++col) {
double accum = 0.0;
const int row_offset = row * dim;
for (int inner = 0; inner < dim; ++inner) {
accum += srcA[row_offset + inner] * srcB[inner * dim + col];
}
dst[row_offset + col] = accum;
}
}
}성능 측정 프레임워크
고해상도 타이머를 사용하여 다양한 차원에서 실행 시간을 기록합니다. 콜드 스타트 효과를 줄이기 위해 워밍업 실행은 생략하고, 각 크기별 단일 측정값을 수집합니다.
int main() {
const std::vector<int> test_sizes = {10, 20, 30, 40, 50,
100, 200, 500, 800, 1000};
std::vector<double> blas_times;
std::vector<double> naive_times;
std::mt19937_64 rng(42);
std::uniform_real_distribution<double> dist(0.0, 1.0);
for (int n : test_sizes) {
const int elements = n * n;
// 메모리 할당
double* buffer_a = new double[elements];
double* buffer_b = new double[elements];
double* buffer_c = new double[elements];
// 난수 초기화
for (int i = 0; i < elements; ++i) {
buffer_a[i] = dist(rng);
buffer_b[i] = dist(rng);
}
// OpenBLAS/LAPACK 측정
auto t_start = std::chrono::high_resolution_clock::now();
optimized_gemm(n, buffer_a, buffer_b, buffer_c);
auto t_end = std::chrono::high_resolution_clock::now();
double t_blas = std::chrono::duration<double>(t_end - t_start).count();
blas_times.push_back(t_blas);
// 수작업 구현 측정
t_start = std::chrono::high_resolution_clock::now();
naive_multiply(n, buffer_a, buffer_b, buffer_c);
t_end = std::chrono::high_resolution_clock::now();
double t_naive = std::chrono::duration<double>(t_end - t_start).count();
naive_times.push_back(t_naive);
std::cout << "n=" << n << ": BLAS=" << t_blas
<< "s, Naive=" << t_naive << "s, "
<< "Ratio=" << t_naive/t_blas << "x\n";
delete[] buffer_a;
delete[] buffer_b;
delete[] buffer_c;
}
return 0;
}컴파일 및 실행
Ubuntu/Debian 환경에서 필요한 패키지 설치:
sudo apt-get install libopenblas-dev liblapack-dev컴파일 명령:
g++ -O2 -o benchmark gemm_benchmark.cpp -lopenblas
./benchmark결과 분석
측정된 실행 시간(초) 및 가속비:
| 차원(n) | OpenBLAS DGEMM | 단순 구현 | 가속비 |
|---|---|---|---|
| 10 | 3.8×10⁻⁵ | 7.0×10⁻⁵ | 1.8× |
| 20 | 3.1×10⁻⁵ | 1.1×10⁻⁴ | 3.6× |
| 30 | 7.8×10⁻⁵ | 4.1×10⁻⁴ | 5.2× |
| 40 | 2.3×10⁻⁴ | 8.0×10⁻⁴ | 3.4× |
| 50 | 3.1×10⁻⁴ | 1.5×10⁻³ | 4.8× |
| 100 | 2.7×10⁻³ | 5.5×10⁻³ | 2.0× |
| 200 | 3.4×10⁻³ | 1.9×10⁻² | 5.5× |
| 500 | 5.1×10⁻² | 3.2×10⁻¹ | 6.4× |
| 800 | 1.9×10⁻¹ | 1.6×10⁰ | 8.0× |
| 1000 | 4.1×10⁻¹ | 2.9×10⁰ | 7.3× |
차원이 커질수록 성능 격차가 확대되며, n ≥ 500 구간에서 6-8배 가속을 달성합니다. 작은 차원에서는 함수 호출 오버헤드와 캐시 적중률 차이로 변동성이 존재합니다.
최적화 기법 분석
BLAS 구현이 우수한 성능을 내는 핵심 요인:
- 블록 분할 (Tiling): L1/L2 캐시에 맞춘 서브행렬 단위 연산으로 메모리 대역폭 활용 극대화
- SIMD 벡터화: AVX2/AVX-512 명령어로 4-8개 배정도 부동소수점 동시 처리
- 패킹 (Packing): 연속 메모리 접근을 위한 행렬 재배열로 TLB 미스 감소
- 마이크로커널 최적화: 레지스터 활용도를 극대화하는 어셈블리 수준 루틴
단순 삼중 루프는 이러한 최적화가 전혀 적용되지 않아, 메모리 계층 구조를 효과적으로 활용하지 못합니다.