לדלג לתוכן

6.1 cuBLAS אלגברה לינארית מואצת הרצאה

בפרקים 0 עד 5 בנינו GPU מתוך אבני הבניין: הבנו את ה-SM (מעבד רב-זרמי - Streaming Multiprocessor) ואת ליבות ה-CUDA וה-Tensor (שיעורים 1.3, 1.5), כתבנו kernels משלנו וטיפסנו בפרויקט 3.5 מ-matmul naive של 3% מה-peak עד shared-memory tiling של 26%, ואפילו שם ראינו את התקרה המקצועית - cuBLAS - חולפת על פנינו פי 3.3 ומגיעה ל-86% מה-peak של FP32. בפרק הזה אנחנו עוברים לצד השני של המתרס: במקום לכתוב kernels, אנחנו קוראים ל-kernels שכתבו המומחים של NVIDIA. cuBLAS היא הספרייה שיושבת מתחת ל-PyTorch, ל-TensorFlow ולכל חישוב אלגברה לינארית צפוף על GPU, והיא נקודת ההשוואה ("reference ceiling") של כל פרויקט GEMM בקורס. בשיעור הזה נלמד להשתמש בה נכון: נתקין ונחווט, נבין את מנגנון בחירת ה-kernel האטום שלה, נכתוב cublasSgemm ל-FP32, נתמודד עם המלכוד מספר 1 של כל מי שנגע ב-cuBLAS - הפריסה ה-column-major מול row-major של C - ולבסוף נדליק את ה-Tensor Cores דרך cublasGemmEx ונמדוד את הזינוק ל-TF32 ול-FP16. את איך ה-Tensor Cores עובדות בפנים כבר יודעים מ-1.5; כאן לומדים איך לקרוא להן במשפט קוד אחד.

מהי cuBLAS - CUDA Basic Linear Algebra Subroutines

ראשי התיבות cuBLAS הם CUDA Basic Linear Algebra Subroutines - מימוש ה-GPU של NVIDIA לתקן הוותיק BLAS (ספריית תת-שגרות אלגברה לינארית בסיסית, בת עשרות שנים מעולם ה-Fortran המדעי). הרעיון פשוט: במקום שכל מפתח יכתוב ויכוונן ביד את פעולות הליבה של אלגברה לינארית - כפל מטריצה-וקטור, כפל מטריצה-מטריצה, פתרון מערכות - cuBLAS מספקת "kernels מכוונים היטב לפעולות אלגברה לינארית נפוצות", מוכנים לקריאה.

התקן BLAS מחולק לשלוש רמות, וכדאי להכיר אותן כי הן חוזרות בשמות הפונקציות:

רמה פעולה אופיינית Arithmetic Intensity דוגמה
BLAS-1 וקטור-וקטור נמוכה (memory-bound) axpy: y = a*x + y
BLAS-2 מטריצה-וקטור נמוכה-בינונית gemv: y = A*x
BLAS-3 מטריצה-מטריצה גבוהה (compute-bound) gemm: C = A*B

הכוכב של cuBLAS, ומוקד השיעור הזה, הוא GEMM - GEneral Matrix Multiply ברמת BLAS-3. זו הפעולה C = alpha * A*B + beta * C, ליבו של כל שכבה צפופה (dense / linear layer) ברשת נוירונים, וגם הפעולה שהarithmetic intensity הגבוהה שלה מאפשרת לרווֹת (saturate) את החומרה - בדיוק מה שראינו בפרויקט 3.5.

כיצד cuBLAS עובדת בפנים - בחירת kernel אטומה. נקודה מהותית שקובעת את כל אופן העבודה מולה: cuBLAS אינה kernel אחד. היא מכילה, כלשון ה-glossary, "מערך רחב של kernels, כל אחד מכוונן במיוחד לטיפוסי נתונים ספציפיים (למשל FP32, FP16), לגדלי מטריצות, ולדורות ארכיטקטורת ה-SM". לפעולה לוגית אחת (GEMM) יש עשרות מימושים קונקרטיים. בזמן ריצה cuBLAS משתמשת בהיוריסטיקות פנימיות (לא ידועות) כדי לבחור את ה-kernel המהיר ביותר ואת פרמטרי הlaunch האופטימליים עבור הבעיה שלכם. אתם לא בוחרים את ה-kernel - הספרייה בוחרת. זו "קופסה שחורה" (black box) כפי שתיארנו כבר ב-1.5: יעילה מאוד, אבל סגורה. אם תרצו לבנות GEMM מותאם בעצמכם, תרדו ל-CUTLASS ול-CuTe (שיעורים 6.3 ו-6.4); כאן אנחנו צרכנים.

המקום באקוסיסטם. cuBLAS היא "היסוד לרוב החישוב הנומרי בעל הביצועים הגבוהים על GPU של NVIDIA", ובפרט היא המנוע שמאיץ את פעולות הליבה של PyTorch. חלוקת העבודה עם אחותה cuDNN (שנפגוש בשיעור הבא): cuBLAS מטפלת באלגברה לינארית כללית/צפופה - הכפלים שבשכבות הצפופות - בעוד cuDNN מתמחה בפרימיטיבים של רשתות נוירונים (קונבולוציה, נורמליזציה, attention).

התקנה וחיווט - setup, handle ו-streams

cuBLAS מגיעה כחלק מ-CUDA Toolkit, כך שאין מה להתקין בנפרד. שני דברים בלבד דרושים בקוד:

  • כותרת: #include <cublas_v2.h>. הסיומת _v2 היא ה-API המודרני; אל תשתמשו ב-cublas.h הישן. הכותרת מגדירה את הטיפוסים והפונקציות שבהם נשתמש.
  • חיווט (linking): בעת הcompilation עם nvcc חובה לקשר לספרייה עם הדגל -lcublas.
nvcc -O2 -arch=sm_90a gemm.cu -o gemm -lcublas

עמוד השדרה של כל תוכנית cuBLAS הוא ה-handle (ידית) - אובייקט אטום מסוג cublasHandle_t שמחזיק את ההקשר (context) של הספרייה: הcaches הפנימיים, בחירות תצורה, וה-stream הפעיל. יוצרים אותו פעם אחת בתחילת התוכנית ומשמידים אותו בסוף:

cublasHandle_t handle;
CUBLAS_CHECK(cublasCreate(&handle));   // create once
// ... all GEMM calls ...
CUBLAS_CHECK(cublasDestroy(handle));   // release at the end

בדיקת שגיאות. בדיוק כפי שעטפנו כל קריאת Runtime API ב-CUDA_CHECK (שיעור 5.3), נעטוף כל קריאת cuBLAS ב-CUBLAS_CHECK. הפונקציות של cuBLAS אינן מחזירות cudaError_t אלא cublasStatus_t נפרד, ולכן צריך מאקרו משלו. מ-CUDA 11.4 קיימת cublasGetStatusString שממירה קוד סטטוס למחרוזת קריאה:

#define CUBLAS_CHECK(call)                                                    \
    do {                                                                      \
        cublasStatus_t st_ = (call);                                          \
        if (st_ != CUBLAS_STATUS_SUCCESS) {                                   \
            fprintf(stderr, "cuBLAS error at %s:%d\n  '%s'\n  -> %s\n",        \
                    __FILE__, __LINE__, #call, cublasGetStatusString(st_));   \
            exit(EXIT_FAILURE);                                               \
        }                                                                     \
    } while (0)

streams. כברירת מחדל cuBLAS רצה על ה-default stream. כדי לשלב GEMM בpipeline מקבילה (כמו זו שבנינו ב-5.3) קושרים את ה-handle ל-stream משלכם. כל הקריאות שאחרי כן יורצו על אותו stream:

cudaStream_t stream;
CUDA_CHECK(cudaStreamCreate(&stream));
CUBLAS_CHECK(cublasSetStream(handle, stream));   // every GEMM from here on runs on this stream

הערה למי שאין GPU מקומי: cuBLAS זמינה גם ב-Google Colab (T4 חינם) ובכל VM ענן עם GPU של NVIDIA. החליפו את sm_90a (H100) ב-sm_75 על T4 או sm_80 על A100.

פונקציית GEMM - cublasSgemm

הפונקציה הראשונה שנכתוב היא cublasSgemm - ה-S היא Single precision (FP32). כל משפחת ה-GEMM ב-cuBLAS עוקבת אחר אותה מוסכמת שמות: האות הראשונה היא הדיוק (S=FP32, D=FP64, H=FP16, C/Z=מרוכב), והשאר gemm. החתימה מממשת את הפעולה המלאה של BLAS-3:

   C  =  alpha * op(A) * op(B)  +  beta * C

כאשר op(X) הוא X או X^T (transpose) לפי דגל, ו-alpha, beta הם סקלרים (שימו לב: מצביעים אליהם, לא ערכים). החתימה המלאה:

cublasSgemm(handle,
            transa, transb,       // whether to transpose A, B
            m, n, k,              // problem dimensions: C is m x n, summed over k
            &alpha,               // pointer to scalar
            A, lda,               // matrix A and its leading dimension
            B, ldb,               // matrix B and its leading dimension
            &beta,                // pointer to scalar
            C, ldc);              // matrix C and its leading dimension

שלושה פרטים שקל לפספס ומפילים מתחילים:

  • הסקלרים הם מצביעים. &alpha ו-&beta הם const float*, לא float. כברירת מחדל הם נקראים מהhost (host pointer mode). לכפל טהור C = A*B נציב alpha = 1.0f, beta = 0.0f.
  • beta דורס או צובר. עם beta = 0 המטריצה C נכתבת מאפס (וגם לא נקראת - כך שאין צורך לאתחל אותה). עם beta = 1 מתקבל C += A*B, שימושי לצבירה.
  • הכל על ה-device. המצביעים A, B, C הם כתובות ב-global memory (זיכרון גלובלי) - תוצרי cudaMalloc, לא מצביעי host. cuBLAS אינה מעתיקה עבורכם; העתקת הקלט והפלט ב-cudaMemcpy באחריותכם.

מה שנשאר להסביר, וזה החלק הקשה, הוא מה בדיוק הם m, n, k, ה-transa/transb, ומעל הכל ה-lda/ldb/ldc - שלושתם נגזרים מהמלכוד הגדול של הפריסה בזיכרון. לשם אנחנו הולכים עכשיו.

המלכוד מספר 1 - column-major מול row-major

זהו, כלשון ה-glossary, "מקור השגיאות הנפוץ ביותר בשימוש ב-cuBLAS": פריסת הנתונים בזיכרון. cuBLAS, בירושה מ-Fortran ומ-BLAS הקלאסי, מצפה למטריצות בסדר column-major (עמודה-ראשית) - העמודות רצופות בזיכרון. C, C++ ו-NumPy (ברירת המחדל) הם row-major (שורה-ראשית) - השורות רצופות. זהו בדיוק ההפך.

   logical matrix         row-major (C/C++)         column-major (cuBLAS/Fortran)
   ┌         ┐
   │ 1  2  3 │           memory:                    memory:
   │ 4  5  6 │           1 2 3 4 5 6                 1 4 2 5 3 6
   └         ┘           (row after row)             (column after column)

leading dimension. לצד כל מטריצה מעבירים ל-cuBLAS פרמטר ld (leading dimension) - הצעד (stride) בזיכרון בין שתי עמודות עוקבות, כלומר מספר השורות של המטריצה כפי ש-cuBLAS רואה אותה (column-major). ה-ld מאפשר לעבוד על תת-מטריצה מתוך הקצאה גדולה יותר: אם המטריצה תופסת את כל הרוחב, ld שווה למספר השורות; אם היא חלון בתוך מערך רחב יותר, ld הוא רוחב המערך המלא. ld שגוי הוא באג קלאסי, ונחזור אליו.

תעלול ה-transpose - להאכיל row-major בלי להעביר

לפני שנרים ידיים ונעתיק (transpose) את המטריצות פיזית, יש זהות אלגברית שחוסכת את כל תנועת הנתונים. היא נשענת על עובדה אחת יפהפייה:

מטריצה row-major זהה סיבית-לסיבית ל-transpose שלה שנקראת column-major. הביטו בדוגמה למעלה: הזיכרון 1 2 3 4 5 6 הוא המטריצה 2x3 שלנו ב-row-major, אבל אם נקרא אותו column-major כמטריצה 3x2, נקבל את העמודות [1,2,3] ו-[4,5,6] - שהן בדיוק A^T. אותם בתים בדיוק, פרשנות אחרת. אין העתקה, אין העברה, רק שינוי נקודת מבט.

מכאן הזהות הבסיסית של אלגברה לינארית:

   C = A * B      <=>      C^T = B^T * A^T

נחבר את השניים. יש לנו בזיכרון מטריצות row-major: A בגודל M x K, B בגודל K x N, ואנחנו רוצים C = A*B בגודל M x N, גם הוא row-major. הטריק:

  1. buffer של A (row-major, M x K) הוא, ב-column-major, המטריצה A^T בגודל K x M. באותו אופן B הוא B^T בגודל N x K, ו-C הוא C^T בגודל N x M.
  2. בעולם ה-column-major שבו cuBLAS חיה, נבקש ממנה לחשב את C^T = B^T * A^T - כלומר מעבירים את B ואת A בסדר הפוך, בלי אף transpose מפורש (transa = transb = N).
  3. cuBLAS כותבת ל-buffer של C את C^T ב-column-major, שהוא בדיוק C שלנו כשקוראים אותו חזרה row-major. קסם, בלי תנועת נתונים.

הנה זה קונקרטית. cuBLAS מחשבת result(m x n) = op(A_cublas)(m x k) * op(B_cublas)(k x n), ואנחנו רוצים ש-result יהיה C^T בגודל N x M. לכן m = N, n = M, k = K, האופרנד הראשון הוא ה-B שלנו והשני הוא ה-A שלנו:

const float alpha = 1.0f, beta = 0.0f;
// dA: row-major M x K,  dB: row-major K x N,  dC: row-major M x N (all on device)
CUBLAS_CHECK(cublasSgemm(handle,
        CUBLAS_OP_N, CUBLAS_OP_N,
        N, M, K,          // m=N, n=M, k=K   (computing C^T)
        &alpha,
        dB, N,            // first operand = our B; ld = N (rows of B^T)
        dA, K,            // second operand = our A; ld = K (rows of A^T)
        &beta,
        dC, N));          // output C^T; ld = N (rows of C^T)

עבור המקרה הריבועי M = K = N (הנפוץ במדידות שלנו), כל שלושת ה-ld שווים ל-N וגם m = n = k = N. זה בדיוק הקוד ששילבנו בשלב 3 של פרויקט 3.5, ועכשיו אתם מבינים מדוע dB בא לפני dA.

הוכחה מספרית קטנה

בואו נאמת על מטריצות זעירות שאפשר לחשב ביד. A בגודל 2x3 ו-B בגודל 3x2, שתיהן row-major:

   A = ┌ 1 2 3 ┐   B = ┌  7  8 ┐        C = A*B = ┌  58   64 ┐
       └ 4 5 6 ┘       │  9 10 │                  └ 139  154 ┘
                       └ 11 12 ┘
   (A's buffer: 1 2 3 4 5 6 ; B's buffer: 7 8 9 10 11 12)

בעולם column-major: buffer של B נקרא כ-B^T בגודל 2x3, buffer של A כ-A^T בגודל 3x2. cuBLAS מחשבת B^T * A^T בגודל 2x2, שהוא C^T, וכותבת ל-buffer של C את הרצף 58 64 139 154. כשקוראים אותו חזרה row-major כמטריצה 2x2 מקבלים בדיוק את C, עם השורות 58 64 ו-139 154. שום transpose לא בוצע - רק החלפנו את סדר האופרנדים ובחרנו ld נכונים.

החלופה, כשבאמת רוצים transpose. אם מטריצה כן צריכה להיות מועברת לוגית, מעבירים CUBLAS_OP_T בדגל המתאים במקום CUBLAS_OP_N. אך שימו לב שזה שינוי לוגי בפעולה, לא פתרון לבעיית ה-row/column-major; לרוב תעדיפו את תעלול החלפת האופרנדים כי הוא חינמי לחלוטין.

דיוק מעורב ו-Tensor Cores - cublasGemmEx

עד כה cublasSgemm רץ על ה-CUDA Cores ב-FP32 טהור, ולכן התקרה שלו היא ~66.9 TFLOPS של ה-H100 (חישבנו זאת ב-3.5). אבל מ-1.5 אנחנו יודעים שרוב עוצמת ה-H100 חבויה ב-Tensor Cores (ליבות טנזור), שתפוקתן גבוהה בסדר גודל. כדי לרתום אותן משתמשים ב-cublasGemmEx - גרסת GEMM גמישה שמפרידה בין טיפוס הנתונים בזיכרון לבין טיפוס החישוב (compute type).

נזכיר את העיקרון מ-1.5: מכפילים בדיוק נמוך וזול, אבל צוברים בדיוק גבוה. cuBLAS חושפת זאת דרך שלושה פרמטרים חדשים - טיפוס לכל מטריצה (cudaDataType), טיפוס חישוב (cublasComputeType_t), ובחירת אלגוריתם. החתימה:

cublasGemmEx(handle, transa, transb, m, n, k,
             &alpha,
             A, Atype, lda,       // data type for each operand separately
             B, Btype, ldb,
             &beta,
             C, Ctype, ldc,
             computeType,          // internal accumulation precision
             algo);                // algorithm / kernel selection

מסלול TF32 - הכי קרוב לחינם. TF32 (TensorFloat-32) הוא הפורמט בן 19 הביט שראינו ב-1.5: מריצים קוד שנכתב ל-FP32, בלי לשנות את הנתונים בזיכרון (הם נשארים float), אך הכפל מתבצע על ה-Tensor Cores בדיוק חתוך. הקלט והפלט נשארים CUDA_R_32F, ורק טיפוס החישוב משתנה ל-CUBLAS_COMPUTE_32F_FAST_TF32:

CUBLAS_CHECK(cublasGemmEx(handle,
        CUBLAS_OP_N, CUBLAS_OP_N,
        N, M, K, &alpha,
        dB, CUDA_R_32F, N,
        dA, CUDA_R_32F, K,
        &beta,
        dC, CUDA_R_32F, N,
        CUBLAS_COMPUTE_32F_FAST_TF32,   // <-- turns on Tensor Cores in TF32
        CUBLAS_GEMM_DEFAULT));

שימו לב לנקודה חשובה: cuBLAS אינה משתמשת ב-TF32 כברירת מחדל. cublasSgemm ו-CUBLAS_COMPUTE_32F נותנים FP32 מדויק על CUDA Cores. כדי לקבל את TF32 חייבים לבחור זאת מפורשות (דרך compute type זה, או לחלופין cublasSetMathMode(handle, CUBLAS_TF32_TENSOR_OP_MATH)). זו בחירה מודעת של דיוק מול מהירות.

מסלול FP16 - התפוקה המקסימלית של Tensor Cores. כאן הנתונים עצמם ב-half precision. הקלט CUDA_R_16F (טיפוס __half, מ-<cuda_fp16.h>), הצבירה עדיין ב-FP32 עם CUBLAS_COMPUTE_32F כדי שהשגיאה לא תתפוצץ לאורך k:

// dAh, dBh, dCh are __half arrays on the device
CUBLAS_CHECK(cublasGemmEx(handle,
        CUBLAS_OP_N, CUBLAS_OP_N,
        N, M, K, &alpha,
        dBh, CUDA_R_16F, N,
        dAh, CUDA_R_16F, K,
        &beta,
        dCh, CUDA_R_16F, N,
        CUBLAS_COMPUTE_32F,             // FP32 accumulation over FP16 input
        CUBLAS_GEMM_DEFAULT));

מתחת למכסה, זו בדיוק ההוראה HMMA שפירקנו ב-1.5: קלט FP16, צבירה FP32. cuBLAS פשוט בוחרת עבורכם את ה-kernel שמנפיק אותה ומזין אליו את הtiles.

תקרת הייחוס - cuBLAS מול הקרנל שלנו

עכשיו נחבר הכל למספרים. זכרו מ-3.5 שכפל שתי מטריצות N x N הוא 2*N^3 פעולות; עבור N = 4096 זה 137.4 GFLOP של עבודה. אותה עבודה בדיוק רצה בכל מסלול; מה שמשתנה הוא התקרה שאליה משווים ואיזה חומרה עובדת. הנה תקרות ה-H100 SXM (צפיפות מלאה, ללא sparsity):

מסלול חומרה דיוק חישוב peak (H100 SXM) יחס ל-FP32
cublasSgemm CUDA Cores FP32 ~66.9 TFLOPS 1x
GemmEx + COMPUTE_32F_FAST_TF32 Tensor Cores TF32 ~494.7 TFLOPS ~7.4x
GemmEx + קלט FP16, COMPUTE_32F Tensor Cores FP16 ~989.4 TFLOPS ~14.8x

כבר מהתקרות רואים למה המעבר ל-Tensor Cores הוא הקפיצה הגדולה בקורס: עצם החלפת compute type מזיזה את הגג פי כמעט 7.4 (TF32) או 14.8 (FP16). כמובן, אף מסלול לא נוגע בתקרה שלו - חלק מהזמן הולך על טעינת הtiles, וככל שהחישוב מהיר יותר כך המשקל היחסי של תנועת הזיכרון גדל (מודל ה-Roofline מ-3.x). מספרים מייצגים על H100 עבור N = 4096, כדוגמה בלבד:

  path             time       throughput       % of that path's ceiling   speedup vs FP32
  ─────────────────────────────────────────────────────────────────────────────────────
  cublasSgemm      2.40 ms    57,300 GFLOP/s   ~86% of 66.9 TFLOPS        1.0x
  GemmEx TF32      0.43 ms   ~320,000 GFLOP/s  ~65% of 494.7 TFLOPS       ~5.6x
  GemmEx FP16      0.23 ms   ~600,000 GFLOP/s  ~61% of 989.4 TFLOPS       ~10.5x

והשוואה למה שכתבנו ביד ב-3.5: ה-kernel ה-tiled שלנו הגיע ל-17,180 GFLOP/s (26% מ-peak FP32). cuBLAS ב-FP32 מהיר ממנו פי 3.3, ו-cuBLAS ב-FP16 מהיר ממנו פי 35 בערך. זה בדיוק תפקידה של cuBLAS בקורס: היא ה"תקרת ייחוס" - המספר שאומר כמה "בשר" נשאר על העצם בכל אופטימיזציה שנכתוב בפרויקט 9.3.

מלכוד הדיוק בבדיקת נכונות. כשעוברים ל-TF32 ובעיקר ל-FP16, אסור להשוות לתוצאת CPU בסובלנות של FP32. הכפל ב-FP16 מאבד מנטיסה, והשגיאה היחסית המצטברת ב-GEMM גדולה בהרבה. סובלנות סבירה: rtol של כ-1e-3 ל-FP32, ~1e-2 ל-TF32, ו-~2e-2 ל-FP16. kernel שנכשל בבדיקה עם סובלנות FP32 אינו בהכרח שגוי - ייתכן שרק בחרתם דיוק נמוך יותר; התאימו את הסובלנות לדיוק שבחרתם.

מתי לא cuBLAS - הגבולות של הקופסה השחורה

cuBLAS מצוינת בדיוק כשהבעיה שלכם היא GEMM (או פעולת BLAS) סטנדרטית בטיפוס ובצורה נפוצים. אבל האטימות שלה היא גם המגבלה שלה, וכדאי לדעת מתי לצאת ממנה:

  • בחירת ה-kernel אינה בשליטתכם. מכיוון שההיוריסטיקות סגורות, הביצועים יכולים להשתנות בין טיפוסי נתונים, צורות מטריצה ודורות GPU באופן שאינכם שולטים בו ישירות. עבור צורות "מוזרות" (מטריצות דקות מאוד, ממדים לא-סטנדרטיים) לפעמים ה-kernel שנבחר אינו האופטימלי.
  • אין fusion. cuBLAS מחשבת GEMM ותו לא. אם אתם רוצים לאחד (fuse) את ה-GEMM עם bias, אקטיבציה, או עם GEMM נוסף לתוך kernel יחיד כדי לחסוך כתיבות ל-global memory, cuBLAS לא תעשה זאת. לשם כך יורדים ל-CUTLASS (שיעור 6.3), שבנויה בדיוק סביב mainloop ו-epilogue הניתנים לאיחוד.
  • לצורך שליטה עדינה יותר יש cublasLt. ספריית ה-Lightweight של cuBLAS (cublasLtMatmul) חושפת יותר בקרה על אלגוריתם, פריסה ו-epilogue מוגבל. היא הגשר בין ה-API הפשוט כאן לבין הכתיבה-מאפס של CUTLASS.

עבור רוב הצרכים, ולכל מי שרוצה "כפל מטריצות מהיר עכשיו", cuBLAS בקריאה אחת היא התשובה הנכונה - וזו בדיוק הסיבה ש-PyTorch נשענת עליה. את בניית ה-kernels מאפס, שם השליטה מלאה, נלמד בהמשך הפרק.

סיכום

  • cuBLAS (CUDA Basic Linear Algebra Subroutines) היא מימוש ה-GPU של NVIDIA לתקן BLAS; היא היסוד לחישוב הנומרי הצפוף על GPU והמנוע שמתחת ל-PyTorch, ומחלקת עבודה עם cuDNN (כללי מול פרימיטיבים של רשתות).
  • cuBLAS אינה kernel אחד אלא מאגר של עשרות kernels מכווננים לטיפוס נתונים, גודל מטריצה ודור SM; בזמן ריצה היוריסטיקות פנימיות (לא ידועות) בוחרות את ה-kernel - אתם צרכנים של "קופסה שחורה", לא מכווננים שלה.
  • ההתקנה מינימלית: #include <cublas_v2.h>, קישור עם -lcublas, יצירת cublasHandle_t עם cublasCreate/cublasDestroy, וקישור אופציונלי ל-stream עם cublasSetStream; עוטפים כל קריאה ב-CUBLAS_CHECK.
  • cublasSgemm מחשבת C = alpha*op(A)*op(B) + beta*C ב-FP32; הסקלרים מועברים כמצביעים, beta=0 דורס את C בלי לקרוא אותה, וכל המצביעים הם כתובות device.
  • המלכוד מספר 1: cuBLAS היא column-major בעוד C/C++/NumPy הם row-major, וכל מטריצה דורשת leading dimension (הצעד בין עמודות); זהו מקור הבאגים הנפוץ ביותר.
  • הזהות C = A*B <=> C^T = B^T*A^T, יחד עם העובדה ש-buffer row-major זהה ל-transpose שלו ב-column-major, מאפשרת להאכיל נתוני row-major ל-cuBLAS בלי כל העברה - פשוט מחליפים את סדר האופרנדים (dB לפני dA) ובוחרים ld נכונים.
  • cublasGemmEx מפרידה טיפוס נתונים מטיפוס חישוב ובכך פותחת את ה-Tensor Cores: CUBLAS_COMPUTE_32F_FAST_TF32 מריץ קוד FP32 ב-TF32 כמעט בחינם, וקלט FP16 עם CUBLAS_COMPUTE_32F מגיע לתפוקת השיא; TF32 אינו ברירת מחדל וחייב בחירה מפורשת.
  • על H100 SXM התקרות הן ~66.9 TFLOPS (FP32), ~494.7 TFLOPS (TF32) ו-~989.4 TFLOPS (FP16), כלומר המעבר ל-Tensor Cores מזיז את הגג פי 7.4 עד 14.8; ככל שהדיוק נמוך יותר, כך יש להרפות את סובלנות בדיקת הנכונות.
  • cuBLAS היא תקרת הייחוס של פרויקטי ה-GEMM (3.5 ו-9.3); כשצריך fusion, צורות לא-סטנדרטיות או שליטה מלאה יורדים ל-cublasLt ולאחר מכן ל-CUTLASS/CuTe (שיעורים 6.3 ו-6.4).