150줄의 C 코드로 NumPy 행렬 곱셈 능가

 (salykova.github.io)
1P by GN⁺ | ★ favorite | 댓글 1개
  • NumPy의 행렬 곱셈은 외부 BLAS 라이브러리에 기대지만, 이 구현은 순수 C와 FMA3·AVX2만으로 단일 스레드와 멀티스레드 성능을 BLAS 수준까지 끌어올리는 것을 목표로 함
  • 성능의 중심은 $C$를 작은 블록으로 나누고 16×6 마이크로커널이 YMM 레지스터 안에서 rank-1 update를 반복해 메모리 접근을 줄이는 구조임
  • 임의 크기 행렬에서는 경계 처리가 병목이 되기 쉬워, 마스크 저장과 0 패딩 버퍼를 조합해 마스크 로드의 성능 저하를 피함
  • 캐시 재사용은 k_c, m_c, n_c 블로킹으로 확보하며, 실제 최고 성능은 스레드 수·커널 크기·타일 크기 튜닝에 크게 좌우됨
  • AVX-512는 더 넓은 CPU 지원을 위해 제외했으므로 AVX-512 CPU에서는 BLAS가 더 빠를 수 있고, OpenBLAS 비교도 AVX-512를 끈 조건에서 수행함

구현 목표와 비교 대상

  • 구현 코드는 sgemm.c에 공개되어 있으며, 최신 프로세서에서 멀티스레드 FP32 행렬 곱셈을 최적화함
  • NumPy는 행렬 곱셈 같은 선형대수 연산을 외부 BLAS 라이브러리에 의존함
    • 예시로 Intel MKL, Accelerate, BLIS, GotoBLAS, OpenBLAS가 있음
    • OpenBLAS, GotoBLAS, BLIS는 C/FORTRAN/Assembly로 작성되고, CPU 마이크로아키텍처별 수동 최적화 행렬 곱셈 구현을 포함함
  • 목표는 저수준 어셈블리 없이 순수 C로 작성하면서도 다음 조건을 만족하는 행렬 곱셈 구현임
    • 임의 행렬 크기에서 동작함
    • 최신 x86-64 프로세서에서 실행됨
    • 기존 BLAS 라이브러리와 경쟁함
    • 코드가 단순하고 확장하기 쉬움
  • 참고 자료는 Simon Boehm의 Fast Multidimensional Matrix Multiplication on CPU from Scratch, Sergey Slotin의 Matrix Multiplication, Geohot의 Can you multiply a matrix?, GotoBLAS·BLIS 관련 논문임

벤치마크 조건과 FLOPS 계산

  • 테스트 환경은 AMD Ryzen 7 9700X, 32GB DDR5 6000 MHz CL36, OpenBLAS 0.3.26, GCC 13.3, Ubuntu 24.04.1 LTS임
  • 컴파일 플래그는 -O3 -march=native -mno-avx512f -fopenmp를 사용함
  • 공정한 비교를 위해 OpenBLAS 설치 시 적절한 TARGET을 설정하고 AVX-512 명령어를 비활성화해야 함
    • Zen4/5 프로세서는 make TARGET=ZEN으로 컴파일함
    • 그렇지 않으면 OpenBLAS가 AVX-512 명령어를 기본 사용함
  • OpenBLAS FP32 행렬 곱셈은 cblas_sgemm API로 실행함
  • 벤치마크는 정사각 행렬을 대상으로 함
    • m=n=k=200부터 m=n=k=10000까지 200 간격으로 평가함
    • 행렬 곱셈을 n_iter번 반복하고, 중앙 실행 시간을 성능 측정에 사용함
  • $M \times K$ 행렬 $A$와 $K \times N$ 행렬 $B$를 곱하면 총 연산량은 $2MNK$ FLOP임
    • 성능은 FLOPS=(2*m*n*k)/exec_time으로 계산함

이론적 한계와 SIMD 기반

  • 최신 x86-64 CPU는 SIMD 확장으로 여러 데이터를 병렬 처리함
  • 주요 명령어는 AVX2와 FMA임
    • 둘 다 256비트 YMM 레지스터를 사용함
    • 각 YMM 레지스터는 32비트 float 8개를 담을 수 있음
  • FMA 명령어 VFMADD231PS는 YMM1 = YMM2 * YMM3 + YMM1 형태의 packed single 연산을 수행함
  • Ryzen 9700X에서 fused multiply-add 처리량은 0.5 cycles/instruction, 즉 사이클당 2개 명령어임
  • 이론적으로 Ryzen 9700X는 단일 코어에서 사이클당 32 FLOP를 수행할 수 있음
    • 계산식은 8 floats × 2(add+mul) × 2(1/TP)임
    • 8코어에서 4.7GHz 지속 클럭을 가정하면 멀티스레드 이론 피크는 1203 FLOPS로 추정됨

기본 구현과 마이크로커널

  • 행렬은 column-major 순서로 저장함
    • A[row][col]은 C 포인터에서 ptr[col*M + row]로 접근함
  • 가장 단순한 구현은 $C$의 모든 행과 열을 순회하며 각 원소마다 $A$의 행과 $B$의 열의 내적을 계산함
  • 고성능 구현의 핵심은 $C$를 $m_R \times n_R$ 부분 행렬로 나누고, 각 부분 행렬을 효율적으로 계산하는 마이크로커널임
  • 커널은 $\bar{C}$를 레지스터에 0으로 초기화한 뒤 $K$ 차원을 따라 반복함
    • $\bar{A}$의 열 벡터와 $\bar{B}$의 행 벡터를 레지스터로 가져옴
    • 두 벡터의 외적을 계산해 $\bar{C}$ 누산기에 더함
    • 각 단계는 rank-1 update임
  • 이 방식은 naive 방식의 메모리 접근량 $2K m_R n_R$과 비교해, 레지스터로 가져오는 원소 수를 $(m_R+n_R)K$로 줄임
  • AVX CPU에는 16개 YMM 레지스터가 있으므로 커널 크기는 다음 제약을 만족해야 함
    • $(m_R/8) \cdot n_R + m_R/8 + 1 \le 16$
    • $m_R$은 8의 배수여야 함
  • 이론적으로는 $m_R$과 $n_R$이 크고 같은 값일수록 메모리 접근 감소가 커지지만, 실제 Ryzen 9700X에서는 16×6 커널이 가장 좋은 성능을 보임
  • 구현은 immintrin.h의 intrinsic을 사용함
    • __m256은 256비트 벡터 타입이며 YMM 레지스터 내용을 나타냄
    • _mm256_loadu_ps로 A 열 벡터를 로드함
    • _mm256_broadcast_ss로 B의 스칼라 값을 8개 float 벡터로 브로드캐스트함
    • _mm256_fmadd_ps로 누산기를 갱신함
    • _mm256_storeu_ps로 결과를 메모리에 저장함
  • 생성된 어셈블리에는 vfmadd231ps와 vbroadcastss 같은 SIMD FMA 명령어가 포함됨

임의 크기 행렬을 위한 패딩

  • 기본 16×6 커널은 $M$과 $N$이 각각 16과 6의 배수일 때 바로 동작함
  • 경계 영역에서 열 수 $n$이 6보다 작으면 저장 루프를 j < n까지만 수행함
  • 행 수 $m$이 16보다 작을 때는 _mm256_storeu_ps가 8개 원소를 한 번에 저장하므로 마스크 저장이 필요함
    • _mm256_maskstore_ps는 마스크 비트가 켜진 원소만 메모리에 저장함
    • 마스크는 겹치는 행 수 $m$에 따라 생성함
  • 경계에서 로드까지 _mm256_maskload_ps로 처리하면 커널 성능이 크게 떨어질 수 있음
    • 마스크 계산 추가 명령어가 오버헤드를 만듦
    • $n$이 컴파일타임 상수가 아니어서 컴파일러가 루프를 효율적으로 언롤하기 어려움
  • 대신 $m \neq m_R$이면 $\bar{A}$를 버퍼에 복사해 0으로 패딩하고, $n \neq n_R$이면 $\bar{B}$도 버퍼에 복사해 0으로 채움
  • 관련 구현은 matmul_pad.h에 있음

캐시 블로킹과 데이터 재사용

  • 레지스터와 DRAM 사이에는 CPU 캐시 계층이 있으며, 최신 데스크톱 CPU는 보통 L1, L2, L3 캐시를 사용함
  • 캐시는 DRAM보다 빠르지만 용량이 제한되어, 전체 $A$, $B$, $C$를 모두 캐시에 담는 방식은 불가능함
  • 행렬을 작은 블록으로 나누어 캐시에 올리고 같은 데이터를 여러 rank-1 update에 재사용하는 방식이 캐시 블로킹 또는 타일링임
  • 단일 스레드 캐시 블로킹은 BLIS 구조와 유사한 5중 루프 형태임
    • 가장 바깥 루프는 $N$ 차원을 따라 $C_j$와 $B_j$ 블록을 만듦
    • 다음 루프는 $K$ 차원을 따라 $A_j$와 $B_p$ 블록을 만듦
    • $B_p$는 패킹되어 $\tilde{B}_p$가 되고, 필요한 경우 0으로 패딩되어 L3 캐시 재사용을 노림
    • 다음 루프는 $M$ 차원을 따라 $C_i$와 $A_j$ 블록을 만들고, $A_j$는 패킹되어 $\tilde{A}_j$가 됨
    • 마지막 두 루프는 캐시 블록을 $m_R \times k_c$, $k_c \times n_R$ 패널로 나누어 커널에 전달함
  • 패킹된 $\tilde{A}_j$와 $\tilde{B}_p$는 서로 다르게 저장됨
    • $\tilde{A}_j$ 내부 패널은 column-major로 저장됨
    • $\tilde{B}_p$ 내부 패널은 row-major로 저장됨
  • 캐시 블로킹 파라미터는 CPU 모델별 캐시 크기에 맞춰 조정해야 함
    • $k_c \times n_c$는 L3 캐시를 채우는 출발점이 됨
    • $m_c \times k_c$는 L2 캐시를 채우는 출발점이 됨
    • $k_c \times n_R$은 L1 캐시를 채우는 출발점이 됨
  • 실제로는 이론값보다 큰 값이 더 좋은 성능을 내는 경우가 많고, CPU가 캐시 배치를 자동 관리하므로 알고리듬 수준에서 루프와 접근 패턴을 설계해야 함
  • 구현은 matmul_cache.h에 있음

커널 미세 최적화

  • __m256 C_buffer[6][2]처럼 배열로 누산기를 정의하는 대신, 누산기 변수를 명시적으로 펼쳐 선언함
  • 이 방식은 GCC가 코드를 더 잘 최적화하고 레지스터 spilling을 피하는 데 도움을 줌
  • 마스크 계산도 벡터 명령어를 사용하도록 바꿈
    • mask[32] 정적 배열을 두고 _mm256_cvtepi8_epi32와 _mm_loadu_si64를 사용함
  • 해당 구현은 matmul_micro.h에 있음

멀티스레딩 전략

  • 병렬화 대상은 산술 연산과 패킹 모두임
  • 마이크로커널 바깥의 5번째, 4번째, 3번째 루프는 캐시 블록 크기 단위로 반복함
    • 모든 스레드를 바쁘게 유지하려면 반복 횟수가 스레드 수 이상이어야 함
    • 입력 행렬 차원은 대략 스레드 수 × 캐시 블록 크기 이상이어야 함
  • Ryzen 9700X 단일 스레드에서 좋은 성능을 보인 캐시 블록 크기는 $n_c=1535$, $m_c=1024$임
    • 8코어를 모두 활용하려면 최소 $\max(m_c,n_c) \times 8 = 1535 \times 8 = 12280$ 크기 차원이 필요함
  • 반대로 마지막 두 루프는 작은 $m_R$, $n_R$ 블록을 반복하므로 병렬화에 적합함
    • 일반적으로 $m_R$, $n_R$은 20보다 작음
    • $m_c$, $n_c$를 코어 수의 배수로 선택하면 작업을 균등하게 나눌 수 있음
  • Ryzen 9700X에서는 두 개의 내부 루프를 #pragma omp parallel for collapse(2) num_threads(NTHREADS)로 함께 병렬화하는 방식이 가장 좋은 성능을 냄
  • 많은 코어를 가진 프로세서, 특히 16코어 초과 환경에서는 중첩 병렬성과 2~3개 루프 병렬화를 고려할 수 있음
  • $\tilde{A}$와 $\tilde{B}$ 패킹도 OpenMP로 병렬화함
    • pack_blockA는 mc를 MR 단위로 순회하며 병렬화함
    • pack_blockB는 nc를 NR 단위로 순회하며 병렬화함
  • 멀티스레드 구현에서 Ryzen 9700X 기준 좋은 성능을 보인 파라미터는 다음과 같음
    • $m_c = m_R \times \text{number of threads} \times 5$
    • $n_c = n_R \times \text{number of threads} \times 50$
  • 최종 멀티스레드 구현은 matmul_parallel.h에 있음

함께 보면 좋은 글 β

GN⁺   [-]
Hacker News 의견들

  • 이 글의 요지가 보통 성능 여지가 남아 있다는 것이라면, 오히려 개선 폭을 과소평가한 편임. 행렬 곱셈 라이브러리에 들어가는 노력은 대부분의 소프트웨어보다 훨씬 큰데도 그렇다
    이미 강하게 최적화된 코드가 아니라면 큰 노력 없이도 기존 코드에서 10~1000배 이상 개선되는 일이 흔함. 대략 중요도 순으로 보면, 알고리즘 선택이 적절한지와 작업 자체를 없앨 수 있는지가 가장 중요하고, 커널 왕복이나 malloc 같은 무거운 작업을 줄일 수 있는지도 크다.
    벡터화는 명시적 벡터 내장 함수도 좋지만, 구조체 배열 대신 배열/구조체의 배열로 데이터를 재구성하는 것만으로도 같은 기계어가 나오는 경우가 많음. 캐시 효율도 중요하고, 병렬 코드에서는 거짓 공유처럼 스레드별 데이터 격리가 안 될 때 더 복잡해짐. 마지막으로 내장 함수나 손수 작성한 어셈블리처럼 하드웨어별 최적화도 가능하다

    • 네트워크 영향도 빼면 안 됨. 한 번은 분산 질의가 네트워크로 약 100만 행을 가져온 뒤 조인해서 5~10행만 남기는 걸 발견해 수백 배 성능 개선을 냈다
      조인이 원격 서버에서 일어나도록 질의를 바꾸고 네트워크로는 5~10행만 보내게 하니 바로 빨라졌음. 고정 오버헤드와 지연 시간은 늘 있지만, 필요한 것보다 훨씬 많은 데이터를 네트워크 연결로 보내면 결국 성능이 망가진다. 지연 시간의 영향을 다룬 “It's the latency, stupid”도 읽을 만함: http://www.stuartcheshire.org/rants/latency.html
      전체적으로는 위 고려사항과 대략적인 순서에 동의함
    • “알고리즘 선택이 적절한가”가 실제로는 카고 컬트가 되어버린 면이 있음. “더 빠른” 알고리즘이 실제 상수항이 끔찍해서, 일을 더 많이 하는 쪽이 오히려 성능이 나은 경우도 많다
      많은 면접이 구현이 왜 느린지 추론하고 벤치마크하며 고치는 방법을 보는 대신, “Google이 그렇게 하니까”식으로 obscure한 알고리즘 암기 퀴즈가 되어버림

  • 흔한 코딩 패턴은 하드웨어에 충분히 특화하지 않아 성능을 많이 남겨둔다. 이 글이 흥미로운 예이고, 또 다른 고전적인 시연으로는 “There's plenty of room at the top”이 있음
    https://www.science.org/doi/10.1126/science.aam9744

  • 이 내용을 이해하려면 BLIS 저장소의 논문들이 정본에 가깝다. 최적화된 BLAS가 성능이 안 나온다고 생각하는 이유를 모르겠고, 충분히 큰 행렬이라면 CPU 피크의 90% 이상을 기대해야 함
    마지막으로 봤을 때 직렬 OpenBLAS는 대체로 MKL과 비슷했고, BLAS는 기본 선형대수 블록으로 matmul이 아니라 GEMM을 구현한다. 보통 벤치마크 프레임워크 대신 numpy를 쓰는 것도 이해가 안 되고, Zen에서는 AMD의 BLAS, 즉 BLIS 기반 구현과 비교해야 한다고 봄. BLIS는 예전에 OpenBLAS보다 병렬화 쪽 이야기가 더 좋았고, AMD BLIS에는 “작은” 차원용 구현 전환도 있는데 현재 OpenBLAS에 있는지는 모르겠다
    마이크로 커널 벡터화에 SIMD 내장 함수가 꼭 필요한 건 아니며, 괜찮은 C 컴파일러는 완전히 벡터화하고 루프도 펼쳐준다. BLIS의 순수 C 마이크로 커널은 적절한 블록 크기에서 Haswell의 손수 최적화 구현 대비 80% 이상 성능을 냄. 차이는 아마 프리페치 때문일 텐데, 정확히 이해하진 못함

    • SIMD 내장 함수와 수동 루프 펼치기는 분명 필요함. 모든 BLAS 라이브러리가 루프를 수동으로 벡터화하고 펼치는 이유가 그거다
      최신 컴파일러도 자동 벡터화와 루프 펼치기를 100% 성공률로 제대로 해내지는 못함

  • 글과 구현은 좋아 보이지만 “비결”이 뭔지 궁금함. OpenBLAS는 이 정확한 문제를 위해 수십 년 동안 어셈블리+C로 최적화되어 왔는데, 어떻게 이길 수 있는 걸까
    캐싱 등을 자세히 다루는데 BLAS가 이런 걸 활용하지 않는 건지, 아니면 특정 프로세서에 더 잘 맞춘 건지 궁금하다

    • OpenBLAS가 특정 최신 아키텍처에 그렇게까지 최적화된 건 아님. 행렬도 그렇게 크지 않았고, numpy에는 cffi 오버헤드가 있다
      성능 차이는 평균 처리량보다 피크 처리량에서 훨씬 두드러졌는데, 피크가 중요한 애플리케이션은 거의 없다. 표시된 벤치마크 코드는 numpy 쪽은 Python 할당자를 지나가고 C 구현은 할당자를 안 지나는 듯해서, 마이크로벤치마크 오류나 불일치를 먼저 확인할 곳임. 많은 numpy 루틴은 제자리 연산을 지원하므로 양쪽 모두 제자리 버전 벤치마크를 명시적으로 봐야 할 것 같다
      numpy에는 하위 구현과 무관하게 실행되는 경계 검사와 오류 처리도 있어 작은 행렬에서는 순수 Python 리스트보다도 매우 느린 이유가 된다. 몇천 사이클의 순수 오버헤드를 더하면 빠르게 만들기 어렵다
      이 구현은 관련 캐시를 포화시키려는 꽤 원칙적인 접근이고, 어떤 의미에서는 뻔하지만 명확한 엔지니어링 개선은 이런 논의에서 강조할 가치가 있음. OpenBLAS도 많은 인력을 들였지만 모든 걸 다 생각했을 가능성은 낮다. 제대로 설명하려면 양쪽 코드 모두에 대한 깊은 분석이 필요함
    • OpenBLAS를 이기는 것은 놀랍지도 않고 전례가 없는 일도 아님. 예를 들어 D 언어의 선형대수 라이브러리 Mir도 몇 년 전에 그랬다 [1]
      C++와 C 구현은 메타프로그래밍 접근 [2], [3]을 보면 됨. 정말 놀라운 건 Matlab, Julia, Mojo 같은 많은 현대 언어가 아직도 OpenBLAS에 의존한다는 점인데, 물론 각자 이유가 있겠지
      [1] Numeric age for D: Mir GLAS is faster than OpenBLAS and Eigen (2016):
      http://blog.mir.dlang.io/glas/benchmark/openblas/2016/09/23/...
      [2] Vastly outperforming LAPACK with C++ metaprogramming (2018):
      https://wordsandbuttons.online/vastly_outperforming_lapack_w...
      [3] Outperforming LAPACK with C metaprogramming (2018):
      https://wordsandbuttons.online/outperforming_lapack_with_c_m...
    • -march=native가 정확한 CPU 모델에 맞춰 컴파일하므로 이점이 있을 수 있음. numpy는 더 범용적이고 오래된 x86-64 대상으로 컴파일됐을 가능성이 크다
      Ryzen CPU에서는 -march=native가 아마 v4를 쓰고, numpy는 v1이나 v2를 목표로 할 듯함
      https://en.wikipedia.org/wiki/X86-64#Microarchitecture_level...
    • numpy 2.0은 여러 마이크로아키텍처에서 SIMD를 더 잘 쓰기 위해 Google highway를 통합하므로, numpy 쪽 비교가 더 좋아질 것임

  • 글도 좋고 벤치마크를 쉽게 재현 가능하게 만든 점도 훌륭함. 내 16코어 Xeon W-2245 3.90GHz에서는 matmul.c가 8192x8192 행렬 곱셈을 gcc -O3로 1.41초, clang -O2로 1.47초에 수행했고, NumPy는 1.07초였다
    AVX-512 커널이면 훨씬 빨라질 거라고 봄. 성능이 아쉬운 또 다른 이유는 OpenMP일 수 있는데, 경험상 pthreads로 스레드 풀을 명시적으로 관리하면 오버헤드를 줄일 수 있다. CPU 개수도 하드코딩 대신 sysconf(_SC_NPROCESSORS_ONLN)을 쓰는 편이 좋음

  • 한쪽은 Python이고 다른 쪽은 C로 부담을 다르게 줄 이유가 없음. 양쪽 모두 C로 작성하고, 하나는 BLAS 라이브러리를 호출하고 다른 하나는 이 구현을 호출하는 식으로 사과 대 사과 비교를 할 수 있었을 것임

    • 여기서는 Python과 비교하는 게 맞다. 요즘 이런 계산을 수행하는 가장 대중적인 방식이 numpy를 쓰는 Python이기 때문임
      오버헤드가 아주 크진 않지만, 이 스레드의 다른 곳에서도 말했듯이 올바르게 호출하는 게 중요하다. 순진한 numpy 코드와 조정된 C 코드를 맞붙이는 건 분명 공정한 비교가 아님

  • 뜨거운 경로는 아니지만, 마스크 생성의 비효율성, 즉 bit_mask 사용이 거슬림. 더 효율적인 방법으로는 {-1,-1,...,0,0,...} 형태의 전역 상수 배열을 만들고 원소 오프셋 16-m, 8-m에서 로드하거나, 상수 벡터 {0,1,2,3,4,...}를 브로드캐스트된 m 및 m-8과 비교하는 방식이 있다
    다만 행렬의 한 열에만 해당하고, 뒤따르는 maskload/maskstore 루프가 훨씬 오래 걸리므로 아주 사소한 꼬집기임. 특히 저장은 Zen 4에서도 여전히 느리고[1], AVX-512 명령은 마스크를 마스크 레지스터에서 받는다는 차이만 있는데도 6배 빠르다. clang은 어차피 시프트를 자동 벡터화하니, 내 제안보다 2~3배 느린 정도일 듯함
    [1]: https://uops.info/table.html?search=vmaskmovps&cb_lat=on&cb_...

    • 글쓴이임. C 코드 최적화와 내장 함수 사용은 정말 처음이라 이 분야 전문가가 아니지만 더 배우고 싶다
      새 관점을 주는 피드백을 정말 고맙게 생각함. “상수 전역 배열을 만들고 로드하기”는 기억상 테스트했을 때 비트 마스크 시프트보다 조금 느렸던 것 같은데, 확실히 하려고 다시 테스트해보겠다. “상수 벡터 {0, 1, 2, 3, 4, ...}를 브로드캐스트된 m과 m-8과 비교”하는 방식은 좋은 아이디어라 시도해보겠음
    • 전역 상수 배열을 만들 때 원소를 int8_t로 두고, 로드하면서 바이트를 int32_t로 부호 확장할 수 있음. _mm_loadu_si64 / _mm256_cvtepi8_epi32 조합은 메모리 피연산자를 쓰는 단일 vpmovsxbd 명령으로 컴파일될 것이다
      이렇게 하면 alignas(32)로 제대로 정렬했을 때 전체 상수 배열이 캐시 라인 하나에 들어간다. 원문 사용 사례에는 마스크가 두 개 필요하므로 두 번째 vpmovsxbd 명령은 확실한 L1D 캐시 적중이 되어 잘 맞음

  • jart의 tinyBLAS는 어떨까
    https://hacks.mozilla.org/2024/04/llamafiles-progress-four-m...
    그리고 https://justine.lol/matmul/

    • 어제 Justine과 활발히 얘기했는데, 그 워크스테이션에서는 이 구현이 tinyBLAS보다 최소 2배 빠른 것 같음. 전체 논의는 Mozilla AI Discord에 있음: https://discord.com/invite/NSnjHmT5xY

  • 벤치마크 말고는 행렬 곱셈 자체를 다중 스레드화하는 이유가 뭘까. 실제로는 곱셈을 사용하는 알고리즘 쪽에서 다중 스레드를 쓰는 편이 더 유리하지 않을까

    • HPC에서는 실제로 보통 그렇게 함. 다만 병렬 BLAS로 교체하는 것만으로 특정한 종류의 R 코드는 쉽게 도움이 될 수 있다
      하지만 HPC 코드는 대개 GEMM이 병목이 아니다

  • 아직 훑어만 봤지만, 이 글은 세부사항과 설명이 많다. 빠른 행렬 곱셈이 아키텍처 고려사항을 반영해 어떻게 구현되는지 꽤 훌륭하게 설명하는 글처럼 보여서 읽을 목록에 넣어둠

답변달기