150 लाइनों के C कोड से NumPy matrix multiplication को पीछे छोड़ा

 (salykova.github.io)
1 पॉइंट द्वारा GN⁺ 2024-07-05 | 1 टिप्पणियां | WhatsApp पर शेयर करें
  • NumPy में matrix multiplication बाहरी BLAS लाइब्रेरी पर निर्भर करती है, लेकिन यह implementation सिर्फ pure C और FMA3·AVX2 के सहारे single-thread और multi-thread performance को BLAS स्तर तक ले जाने का लक्ष्य रखती है
  • performance का केंद्र $C$ को छोटे blocks में बाँटना है, और 16×6 microkernel YMM registers के भीतर rank-1 update को बार-बार चलाकर memory access कम करता है
  • arbitrary size matrices में boundary handling आसानी से bottleneck बन सकती है, इसलिए masked store और zero-padding buffer को मिलाकर masked load की performance गिरावट से बचा गया है
  • cache reuse को k_c, m_c, n_c blocking से सुनिश्चित किया जाता है, और वास्तविक peak performance thread count, kernel size, tile size tuning पर बहुत निर्भर करती है
  • wider CPU support के लिए AVX-512 को शामिल नहीं किया गया, इसलिए AVX-512 CPU पर BLAS तेज हो सकता है, और OpenBLAS की तुलना भी AVX-512 बंद करके की गई है

implementation का लक्ष्य और comparison target

  • implementation code sgemm.c में उपलब्ध है, और यह modern processors पर multithreaded FP32 matrix multiplication को optimize करता है
  • NumPy matrix multiplication जैसी linear algebra operations के लिए बाहरी BLAS लाइब्रेरी पर निर्भर करता है
    • उदाहरण के तौर पर Intel MKL, Accelerate, BLIS, GotoBLAS, OpenBLAS हैं
    • OpenBLAS, GotoBLAS, BLIS को C/FORTRAN/Assembly में लिखा गया है, और इनमें CPU microarchitecture के हिसाब से manually optimized matrix multiplication implementations शामिल हैं
  • लक्ष्य यह है कि low-level assembly के बिना pure C में लिखा जाए और फिर भी नीचे की शर्तें पूरी हों
    • arbitrary matrix sizes पर काम करे
    • modern x86-64 processors पर चले
    • मौजूदा BLAS लाइब्रेरी से मुकाबला करे
    • code सरल हो और विस्तार करना आसान हो
  • संदर्भ सामग्री में Simon Boehm का Fast Multidimensional Matrix Multiplication on CPU from Scratch, Sergey Slotin का Matrix Multiplication, Geohot का Can you multiply a matrix?, और GotoBLAS·BLIS से जुड़े papers शामिल हैं

benchmark conditions और FLOPS calculation

  • test environment AMD Ryzen 7 9700X, 32GB DDR5 6000 MHz CL36, OpenBLAS 0.3.26, GCC 13.3, Ubuntu 24.04.1 LTS है
  • compile flags के रूप में -O3 -march=native -mno-avx512f -fopenmp का उपयोग किया गया है
  • fair comparison के लिए OpenBLAS install करते समय सही TARGET सेट करना और AVX-512 instructions को disable करना जरूरी है
    • Zen4/5 processors को make TARGET=ZEN से compile किया जाता है
    • नहीं तो OpenBLAS default रूप से AVX-512 instructions का उपयोग करता है
  • OpenBLAS FP32 matrix multiplication को cblas_sgemm API से चलाया गया है
  • benchmark square matrices पर किया गया है
    • m=n=k=200 से m=n=k=10000 तक 200 के interval पर evaluate किया गया है
    • matrix multiplication को n_iter बार दोहराया जाता है, और median execution time को performance measurement के लिए उपयोग किया जाता है
  • $M \times K$ matrix $A$ और $K \times N$ matrix $B$ को गुणा करने पर कुल operation count $2MNK$ FLOP होता है
    • performance को FLOPS=(2*m*n*k)/exec_time से calculate किया जाता है

theoretical limits और SIMD आधार

  • modern x86-64 CPU SIMD extension के जरिए कई data items को parallel में process करते हैं
  • मुख्य instructions AVX2 और FMA हैं
    • दोनों 256-bit YMM registers का उपयोग करते हैं
    • हर YMM register में 32-bit float के 8 मान आ सकते हैं
  • FMA instruction VFMADD231PS YMM1 = YMM2 * YMM3 + YMM1 रूप में packed single operation करती है
  • Ryzen 9700X पर fused multiply-add throughput 0.5 cycles/instruction है, यानी प्रति cycle 2 instructions
  • सैद्धांतिक रूप से Ryzen 9700X single core पर प्रति cycle 32 FLOP कर सकता है
    • calculation है 8 floats × 2(add+mul) × 2(1/TP)
    • 8 cores पर 4.7GHz sustained clock मानें तो multithreaded theoretical peak 1203 FLOPS अनुमानित होती है

basic implementation और microkernel

  • matrices को column-major क्रम में store किया जाता है
    • A[row][col] को C pointer में ptr[col*M + row] से access किया जाता है
  • सबसे सरल implementation $C$ की सभी rows और columns पर iterate करके हर element के लिए $A$ की row और $B$ की column का dot product निकालती है
  • high-performance implementation का core यह है कि $C$ को $m_R \times n_R$ submatrices में बाँटा जाए, और हर submatrix को efficiently calculate करने वाला microkernel इस्तेमाल किया जाए
  • kernel पहले $\bar{C}$ को registers में 0 से initialize करता है और फिर $K$ dimension के साथ iterate करता है
    • $\bar{A}$ के column vector और $\bar{B}$ के row vector को registers में लाता है
    • दोनों vectors का outer product निकालकर $\bar{C}$ accumulators में जोड़ता है
    • हर step एक rank-1 update है
  • यह तरीका naive approach के memory access $2K m_R n_R$ की तुलना में, registers में लाए जाने वाले elements की संख्या को $(m_R+n_R)K$ तक घटा देता है
  • AVX CPU में 16 YMM registers होते हैं, इसलिए kernel size को नीचे की constraint पूरी करनी होती है
    • $(m_R/8) \cdot n_R + m_R/8 + 1 \le 16$
    • $m_R$ को 8 का multiple होना चाहिए
  • सैद्धांतिक रूप से $m_R$ और $n_R$ जितने बड़े और बराबर होंगे, memory access reduction उतनी अधिक होगी, लेकिन वास्तविक Ryzen 9700X पर 16×6 kernel ने सबसे अच्छा performance दिखाया
  • implementation immintrin.h के intrinsics का उपयोग करती है
    • __m256 256-bit vector type है और YMM register contents को दर्शाता है
    • _mm256_loadu_ps से A का column vector load किया जाता है
    • _mm256_broadcast_ss से B के scalar value को 8-float vector में broadcast किया जाता है
    • _mm256_fmadd_ps से accumulators update होते हैं
    • _mm256_storeu_ps से result memory में store होता है
  • generated assembly में vfmadd231ps और vbroadcastss जैसे SIMD FMA instructions शामिल होते हैं

arbitrary size matrices के लिए padding

  • basic 16×6 kernel सीधे तब काम करता है जब $M$ और $N$ क्रमशः 16 और 6 के multiples हों
  • boundary region में यदि columns की संख्या $n$ 6 से कम हो, तो store loop को j < n तक ही चलाया जाता है
  • जब rows की संख्या $m$ 16 से कम हो, तब _mm256_storeu_ps एक बार में 8 elements store करता है, इसलिए masked store की जरूरत पड़ती है
    • _mm256_maskstore_ps सिर्फ उन्हीं memory elements को store करता है जिनके mask bit enabled हों
    • mask overlapping rows की संख्या $m$ के आधार पर बनाया जाता है
  • अगर boundary पर load भी _mm256_maskload_ps से किया जाए तो kernel performance काफी गिर सकती है
    • mask calculation के लिए अतिरिक्त instructions overhead बढ़ाते हैं
    • $n$ compile-time constant नहीं होने से compiler loop को efficiently unroll नहीं कर पाता
  • इसके बजाय, अगर $m \neq m_R$ हो तो $\bar{A}$ को buffer में copy करके 0 से pad किया जाता है, और अगर $n \neq n_R$ हो तो $\bar{B}$ को भी buffer में copy करके 0 से भर दिया जाता है
  • संबंधित implementation matmul_pad.h में है

cache blocking और data reuse

  • registers और DRAM के बीच CPU cache hierarchy होती है, और modern desktop CPU आमतौर पर L1, L2, L3 cache का उपयोग करते हैं
  • cache DRAM से तेज होती है, लेकिन उसकी capacity सीमित होती है, इसलिए पूरे $A$, $B$, $C$ को cache में रख पाना संभव नहीं है
  • matrices को छोटे blocks में बाँटकर cache में लाना और उसी data को कई rank-1 updates में reuse करना ही cache blocking या tiling है
  • single-thread cache blocking, BLIS structure जैसी 5-level loop form का उपयोग करता है
    • सबसे outer loop $N$ dimension के along $C_j$ और $B_j$ blocks बनाता है
    • अगला loop $K$ dimension के along $A_j$ और $B_p$ blocks बनाता है
    • $B_p$ को pack करके $\tilde{B}_p$ बनाया जाता है, और जरूरत पड़ने पर 0-padding देकर L3 cache reuse का लक्ष्य रखा जाता है
    • अगला loop $M$ dimension के along $C_i$ और $A_j$ blocks बनाता है, और $A_j$ को pack करके $\tilde{A}_j$ बनाया जाता है
    • अंतिम दो loops cache blocks को $m_R \times k_c$, $k_c \times n_R$ panels में बाँटकर kernel तक भेजते हैं
  • packed $\tilde{A}_j$ और $\tilde{B}_p$ अलग-अलग तरीके से store होते हैं
    • $\tilde{A}_j$ के internal panels column-major में store होते हैं
    • $\tilde{B}_p$ के internal panels row-major में store होते हैं
  • cache blocking parameters को CPU model की cache size के हिसाब से tune करना पड़ता है
    • $k_c \times n_c$ L3 cache भरने का starting point बनता है
    • $m_c \times k_c$ L2 cache भरने का starting point बनता है
    • $k_c \times n_R$ L1 cache भरने का starting point बनता है
  • व्यवहार में अक्सर theoretical values से बड़े values बेहतर performance देते हैं, और CPU cache placement को खुद manage करता है, इसलिए algorithm स्तर पर loops और access patterns को design करना जरूरी होता है
  • implementation matmul_cache.h में है

kernel micro-optimizations

  • __m256 C_buffer[6][2] जैसे array के रूप में accumulators define करने के बजाय, accumulator variables को explicitly unroll करके declare किया गया है
  • इस तरीके से GCC code को बेहतर optimize कर पाता है और register spilling से बचने में मदद मिलती है
  • mask calculation को भी vector instructions का उपयोग करने के लिए बदला गया है
    • mask[32] static array रखा गया है और _mm256_cvtepi8_epi32 तथा _mm_loadu_si64 का उपयोग किया जाता है
  • यह implementation matmul_micro.h में है

multithreading strategy

  • parallelization का लक्ष्य arithmetic operations और packing दोनों हैं
  • microkernel के बाहर वाली 5वीं, 4वीं, 3वीं loops cache block size unit पर iterate करती हैं
    • सभी threads को व्यस्त रखने के लिए iteration count thread count से कम नहीं होना चाहिए
    • input matrix dimensions लगभग thread count × cache block size या उससे बड़े होने चाहिए
  • Ryzen 9700X single-thread पर अच्छा performance देने वाले cache block sizes $n_c=1535$, $m_c=1024$ थे
    • सभी 8 cores का उपयोग करने के लिए कम से कम $\max(m_c,n_c) \times 8 = 1535 \times 8 = 12280$ आकार का dimension चाहिए
  • इसके उलट, अंतिम दो loops छोटे $m_R$, $n_R$ blocks पर iterate करती हैं, इसलिए वे parallelization के लिए उपयुक्त हैं
    • सामान्यतः $m_R$, $n_R$ 20 से छोटे होते हैं
    • अगर $m_c$, $n_c$ को core count के multiples के रूप में चुना जाए, तो काम को बराबर बाँटा जा सकता है
  • Ryzen 9700X पर दो inner loops को #pragma omp parallel for collapse(2) num_threads(NTHREADS) के साथ एकसाथ parallelize करने का तरीका सबसे अच्छा performance देता है
  • ज्यादा cores वाले processors, खासकर 16 cores से ऊपर, में nested parallelism और 2~3 loops parallelization पर विचार किया जा सकता है
  • $\tilde{A}$ और $\tilde{B}$ packing को भी OpenMP से parallelize किया गया है
    • pack_blockA mc को MR unit पर iterate करते हुए parallelize करता है
    • pack_blockB nc को NR unit पर iterate करते हुए parallelize करता है
  • multithreaded implementation में Ryzen 9700X पर अच्छा performance देने वाले parameters इस प्रकार थे
    • $m_c = m_R \times \text{number of threads} \times 5$
    • $n_c = n_R \times \text{number of threads} \times 50$
  • अंतिम multithreaded implementation matmul_parallel.h में है

संबंधित पढ़ाई

1 टिप्पणियां

 
GN⁺ 2024-07-05
Hacker News की राय

  • अगर इस लेख का सार यह है कि आम तौर पर performance की गुंजाइश बची रहती है, तो यह सुधार की सीमा को बल्कि कम करके आंकता है। ऐसा तब भी है जब matrix multiplication libraries में लगने वाली मेहनत अधिकतर software की तुलना में कहीं ज्यादा होती है
    अगर code पहले से बहुत ज्यादा optimized नहीं है, तो अक्सर बिना बहुत बड़ी मेहनत के मौजूदा code में 10~1000 गुना या उससे भी ज्यादा सुधार हो जाता है। मोटे तौर पर महत्व के क्रम में देखें तो सबसे अहम है कि algorithm का चयन सही है या नहीं और क्या काम को ही हटाया जा सकता है; kernel round trips या malloc जैसे भारी operations को घटा पाना भी बड़ा असर डालता है।
    vectorization के लिए explicit vector intrinsics अच्छे हैं, लेकिन कई बार data को structures की array के बजाय arrays/structures की array के रूप में फिर से organize कर देने भर से वही machine code निकल आता है। cache efficiency भी महत्वपूर्ण है, और parallel code में मामला तब और जटिल हो जाता है जब false sharing जैसी वजहों से per-thread data isolation नहीं हो पाता। अंत में, intrinsics या हाथ से लिखी assembly जैसी hardware-specific optimizations भी संभव हैं

    • network का असर भी नहीं छोड़ना चाहिए। एक बार मैंने पाया कि एक distributed query network के जरिए करीब 10 लाख rows ला रही थी और फिर join करके सिर्फ 5~10 rows बचा रही थी; इसे ठीक करके performance में सैकड़ों गुना सुधार मिला
      query को इस तरह बदला कि join remote server पर हो और network से सिर्फ 5~10 rows भेजी जाएं, तो यह तुरंत तेज हो गई। fixed overhead और latency हमेशा रहते हैं, लेकिन जरूरत से कहीं ज्यादा data network connection से भेजेंगे तो अंततः performance खराब होगी। latency के असर पर “It's the latency, stupid” भी पढ़ने लायक है: http://www.stuartcheshire.org/rants/latency.html
      कुल मिलाकर, ऊपर की बातों और उनके मोटे क्रम से सहमत हूं
    • “algorithm का चयन सही है या नहीं” असल में कुछ हद तक cargo cult बन गया है। “तेज” algorithm के वास्तविक constant factors इतने खराब हो सकते हैं कि ज्यादा काम करने वाला तरीका भी कई बार बेहतर performance दे देता है
      कई interviews में implementation धीमा क्यों है, यह reason करके, benchmark करके और ठीक करने का तरीका देखने के बजाय, “Google ऐसा करता है इसलिए” वाली obscure algorithms की रटंत quiz बन जाती है

  • आम coding patterns hardware के लिए पर्याप्त रूप से specialized नहीं होते, इसलिए काफी performance अधूरी छोड़ देते हैं। यह लेख इसका दिलचस्प उदाहरण है, और एक और classic demo “There's plenty of room at the top” है
    https://www.science.org/doi/10.1126/science.aam9744

  • इसे समझने के लिए BLIS repository के papers लगभग canonical reference हैं। मुझे समझ नहीं आता कि कोई क्यों सोचेगा कि optimized BLAS performance नहीं देगा; matrix पर्याप्त बड़ा हो तो CPU peak के 90% से ज्यादा की उम्मीद करनी चाहिए
    पिछली बार जब देखा था, serial OpenBLAS आम तौर पर MKL जैसा ही था, और BLAS basic linear algebra blocks के तौर पर matmul नहीं बल्कि GEMM implement करता है। आम benchmark framework के बजाय numpy इस्तेमाल करना भी समझ नहीं आता, और Zen पर AMD के BLAS, यानी BLIS-based implementation से compare करना चाहिए। BLIS में पहले OpenBLAS की तुलना में parallelization की कहानी बेहतर थी, और AMD BLIS में “छोटे” dimensions के लिए implementation switching भी है, लेकिन पता नहीं अभी OpenBLAS में है या नहीं
    micro-kernel vectorization के लिए SIMD intrinsics जरूरी नहीं हैं; अच्छा C compiler पूरी तरह vectorize कर देता है और loops भी unroll कर देता है। BLIS का pure C micro-kernel उचित block size पर Haswell के hand-optimized implementation की तुलना में 80% से ज्यादा performance देता है। फर्क शायद prefetching की वजह से होगा, हालांकि मैं ठीक से समझ नहीं पाया

    • SIMD intrinsics और manual loop unrolling निश्चित रूप से जरूरी हैं। सभी BLAS libraries loops को manually vectorize और unroll करती हैं, वजह यही है
      modern compilers भी auto-vectorization और loop unrolling को 100% success rate के साथ सही तरह नहीं कर पाते

  • लेख और implementation अच्छे लगते हैं, लेकिन जिज्ञासा है कि “रहस्य” क्या है। OpenBLAS इस सटीक समस्या के लिए दशकों से assembly+C में optimized होता आया है, तो इसे कैसे हराया जा सकता है
    caching वगैरह पर विस्तार से बात की गई है, लेकिन सोच रहा हूँ कि क्या BLAS इन चीज़ों का फायदा नहीं उठाता, या फिर यह किसी खास processor के लिए बेहतर tuned है

    • OpenBLAS किसी खास modern architecture के लिए इतना ज़्यादा optimized नहीं है। matrix भी इतनी बड़ी नहीं थी, और numpy में cffi overhead है
      performance gap average throughput की तुलना में peak throughput में कहीं ज़्यादा साफ दिखा, जबकि peak महत्वपूर्ण होने वाली applications बहुत कम होती हैं। दिखाए गए benchmark code में numpy वाली तरफ Python allocator से होकर जाता दिखता है और C implementation allocator से नहीं गुजरती, इसलिए microbenchmark error या mismatch की जाँच सबसे पहले यहीं करनी चाहिए। कई numpy routines in-place operations support करती हैं, इसलिए दोनों तरफ in-place versions के benchmarks साफ़ तौर पर देखने चाहिए लगते हैं
      numpy में lower implementation से independent चलने वाली bounds checking और error handling भी होती है, जो छोटी matrices में उसके pure Python lists से भी बहुत धीमा होने की वजह बनती है। कुछ हज़ार cycles का pure overhead जोड़ दें तो तेज़ बनाना मुश्किल हो जाता है
      यह implementation relevant cache को saturate करने की काफ़ी principled approach है, और एक अर्थ में obvious है, लेकिन ऐसे discussions में साफ़ engineering improvements को highlight करना वाजिब है। OpenBLAS में भी बहुत manpower लगा है, पर यह मानना मुश्किल है कि उन्होंने सब कुछ सोच लिया होगा। ठीक से समझाने के लिए दोनों codebases का गहरा analysis चाहिए
    • OpenBLAS को हराना न तो चौंकाने वाली बात है, न ही बिना मिसाल की। उदाहरण के लिए D भाषा की linear algebra library Mir ने भी कुछ साल पहले ऐसा किया था [1]
      C++ और C implementations के लिए metaprogramming approach [2], [3] देख सकते हैं। सच में हैरानी की बात यह है कि Matlab, Julia, Mojo जैसी कई modern languages अब भी 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 exact CPU model के हिसाब से compile करता है, इसलिए फायदा हो सकता है। numpy के अधिक generic और पुराने x86-64 target के लिए compile होने की संभावना ज़्यादा है
      Ryzen CPU पर -march=native शायद v4 इस्तेमाल करता है, और numpy शायद v1 या v2 को target करता होगा
      https://en.wikipedia.org/wiki/X86-64#Microarchitecture_level...
    • numpy 2.0 कई microarchitectures पर SIMD का बेहतर उपयोग करने के लिए Google highway integrate करता है, इसलिए numpy वाली comparison बेहतर हो जाएगी

  • लेख भी अच्छा है और benchmarks को आसानी से reproduce करने लायक बनाया गया है, यह भी शानदार है। मेरे 16-core Xeon W-2245 3.90GHz पर matmul.c ने 8192x8192 matrix multiplication gcc -O3 के साथ 1.41 सेकंड में, clang -O2 के साथ 1.47 सेकंड में की, और NumPy ने 1.07 सेकंड लिया
    मुझे लगता है AVX-512 kernel हो तो काफ़ी तेज़ होगा। performance कम लगने की एक और वजह OpenMP हो सकती है; अनुभव से, pthreads के साथ thread pool को explicitly manage करने पर overhead घट सकता है। CPU count hardcode करने के बजाय sysconf(_SC_NPROCESSORS_ONLN) इस्तेमाल करना भी बेहतर होगा

  • एक तरफ Python और दूसरी तरफ C रखकर अलग-अलग burden देने की वजह नहीं है। दोनों को C में लिखा जा सकता था, जहाँ एक BLAS library को call करे और दूसरा इस implementation को, ताकि apples-to-apples comparison हो सके

    • यहाँ Python से comparison करना सही है। आजकल ऐसी computation करने का सबसे popular तरीका numpy इस्तेमाल करने वाला Python ही है
      overhead बहुत बड़ा नहीं है, लेकिन जैसा इस thread में दूसरी जगह भी कहा गया है, सही तरीके से call करना महत्वपूर्ण है। naïve numpy code को tuned C code से भिड़ाना निश्चित रूप से fair comparison नहीं है

  • यह hot path नहीं है, लेकिन mask बनाने की अक्षमता, यानी bit_mask का इस्तेमाल, खटकता है। ज्यादा efficient तरीके के तौर पर {-1,-1,...,0,0,...} फॉर्म का global constant array बनाकर element offset 16-m, 8-m से load किया जा सकता है, या constant vector {0,1,2,3,4,...} की तुलना broadcast किए गए m और m-8 से की जा सकती है
    हालांकि यह matrix के सिर्फ एक column से जुड़ा है, और उसके बाद आने वाला maskload/maskstore loop कहीं ज्यादा समय लेता है, इसलिए यह बहुत मामूली nitpick है। खासकर store Zen 4 पर भी अभी धीमा है[1], और AVX-512 instruction सिर्फ इतना अलग है कि mask को mask register से लेता है, फिर भी 6 गुना तेज है। clang वैसे भी shift को auto-vectorize कर देता है, इसलिए यह मेरे सुझाव से शायद 2–3 गुना धीमा ही होगा
    [1]: https://uops.info/table.html?search=vmaskmovps&cb_lat=on&cb_...

    • मैं लेखक हूँ। C code optimization और intrinsics का इस्तेमाल सचमुच पहली बार कर रहा हूँ, इसलिए इस field का expert नहीं हूँ, लेकिन और सीखना चाहता हूँ
      नया perspective देने वाले feedback के लिए सच में आभारी हूँ। “constant global array बनाकर load करना” मुझे याद है कि test करने पर bit mask shift से थोड़ा धीमा लगा था, लेकिन पक्का करने के लिए फिर से test करूँगा। “constant vector {0, 1, 2, 3, 4, ...} की तुलना broadcast किए गए m और m-8 से करना” अच्छा idea है, इसे try करूँगा
    • global constant array बनाते समय elements को int8_t रखकर, load करते हुए byte को int32_t में sign-extend किया जा सकता है। _mm_loadu_si64 / _mm256_cvtepi8_epi32 combination memory operand का इस्तेमाल करने वाले single vpmovsxbd instruction में compile होगा
      ऐसा करने पर, alignas(32) से ठीक से align करने पर पूरा constant array एक cache line में आ जाता है। original use case में दो masks चाहिए, इसलिए दूसरा vpmovsxbd instruction निश्चित L1D cache hit होगा और अच्छी तरह fit बैठता है

  • jart का tinyBLAS कैसा रहेगा
    https://hacks.mozilla.org/2024/04/llamafiles-progress-four-m...
    और https://justine.lol/matmul/

    • कल Justine से काफी सक्रिय बातचीत हुई, और उस workstation पर यह implementation tinyBLAS से कम से कम 2 गुना तेज लगता है। पूरी discussion Mozilla AI Discord पर है: https://discord.com/invite/NSnjHmT5xY

  • benchmark के अलावा matrix multiplication को खुद multi-threaded बनाने की वजह क्या है। असल में क्या multiplication इस्तेमाल करने वाले algorithm की तरफ multi-threading करना ज्यादा फायदेमंद नहीं होगा

    • HPC में सच में आम तौर पर ऐसा ही किया जाता है। हालांकि parallel BLAS से replace करना भर कुछ खास तरह के R code में आसानी से मदद कर सकता है
      लेकिन HPC code में आम तौर पर GEMM bottleneck नहीं होता

  • अभी सिर्फ सरसरी तौर पर देखा है, लेकिन इस article में काफी details और explanation हैं। तेज matrix multiplication को architecture considerations ध्यान में रखकर कैसे implement किया जाता है, यह काफी अच्छी तरह समझाने वाला article लगता है, इसलिए reading list में डाल दिया है