לדלג לתוכן

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

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

בתרגול הזה תלמדו להשתמש ב-cuBLAS נכון, על העור שלכם. תתחילו מכפל מטריצות row-major מול cuBLAS ה-column-major, תפתרו את מלכוד הפריסה עם תעלול החלפת האופרנדים ותאמתו מול ייחוס CPU; תשוו את התפוקה של cuBLAS ל-kernel ה-tiled שכתבתם בפרויקט 3.5 ותמדדו את הפער מהתקרה של H100; תדליקו את ה-Tensor Cores דרך cublasGemmEx ב-TF32 וב-FP16 ותמדדו את הזינוק; ולבסוף - תשברו בכוונה את ה-leading dimension כדי לראות איך נראית התקלה, ואז תתקנו אותה. עבדו לפי הסדר; כל תרגיל נשען על קודמו. כל התרגילים רצים על GPU של NVIDIA; אם אין לכם כרטיס מקומי השתמשו ב-Google Colab (T4 חינם) והחליפו את דגל הארכיטקטורה.

הכנה

צרו קובץ עבודה cublas_gemm.cu. בראשו שימו את הכותרות, שני מאקרו הבדיקה, ומימוש ייחוס וכלי מדידה. חלק מהתשתית (מד-זמן, בודק נכונות, מחשבון GFLOP/s) זהה לזו של פרויקט 3.5 - העתיקו אותה משם אם כתבתם.

#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <cuda_runtime.h>
#include <cuda_fp16.h>
#include <cublas_v2.h>

#define CUDA_CHECK(call)                                                       \
    do {                                                                       \
        cudaError_t err_ = (call);                                             \
        if (err_ != cudaSuccess) {                                             \
            fprintf(stderr, "CUDA error at %s:%d\n  '%s'\n  -> %s\n",          \
                    __FILE__, __LINE__, #call, cudaGetErrorString(err_));      \
            exit(EXIT_FAILURE);                                                \
        }                                                                      \
    } while (0)

#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)

static const double H100_FP32_PEAK_GFLOPS = 66900.0;   // FP32 CUDA Cores

// CPU reference (row-major): C = A * B, accumulates in double
void matmulCPU(const float* A, const float* B, float* C, int N) {
    for (int r = 0; r < N; r++)
        for (int c = 0; c < N; c++) {
            double s = 0.0;
            for (int k = 0; k < N; k++) s += (double)A[r*N+k] * (double)B[k*N+c];
            C[r*N+c] = (float)s;
        }
}

// correctness check with relative tolerance
bool checkResult(const float* ref, const float* got, int N, float rtol) {
    double maxRel = 0.0;
    for (int i = 0; i < N*N; i++) {
        double denom = fmax(fabs((double)ref[i]), 1e-6);
        double rel = fabs((double)got[i] - (double)ref[i]) / denom;
        if (rel > maxRel) maxRel = rel;
    }
    printf("  max relative error = %.3e (tolerance %.1e)  -> %s\n",
           maxRel, rtol, (maxRel <= rtol ? "PASS" : "FAIL"));
    return maxRel <= rtol;
}

void report(const char* name, int N, float ms, double peak) {
    double flops  = 2.0 * (double)N * N * N;
    double gflops = flops / (ms / 1000.0) / 1e9;
    printf("  %-14s N=%d  %8.3f ms  %10.1f GFLOP/s  %5.2f%% peak\n",
           name, N, ms, gflops, 100.0 * gflops / peak);
}

הcompilation (על H100; החליפו ל-sm_75 על T4, sm_80 על A100):

nvcc -O2 -arch=sm_90a cublas_gemm.cu -o cublas_gemm -lcublas

תרגיל 1 - כפל row-major נכון עם cublasSgemm

מטרתכם: להכפיל שתי מטריצות row-major ולקבל תוצאה נכונה מ-cuBLAS ה-column-major, ולהוכיח זאת מול ייחוס CPU.

  1. הקצו והאתחלו על ה-host שלוש מטריצות float ריבועיות N x N (התחילו קטן, N = 512, כדי שהייחוס על ה-CPU יסתיים מהר). מלאו את A ו-B בערכים אקראיים בטווח [-1, 1].
  2. העתיקו את A ו-B ל-dA, dB ב-device עם cudaMemcpy; הקצו dC.
  3. צרו cublasHandle_t עם cublasCreate. הגדירו alpha = 1.0f, beta = 0.0f.
  4. קראו ל-cublasSgemm בעזרת תעלול החלפת האופרנדים: העבירו CUBLAS_OP_N, CUBLAS_OP_N, ממדים N, N, N, ואת dB לפני dA, עם ld = N לכל אחת. העתיקו את dC חזרה ל-host.
  5. הריצו matmulCPU על אותם A, B והשוו עם checkResult בסובלנות 1e-3. ודאו PASS.
  6. הניסוי המלמד: קראו שוב ל-cublasSgemm, אבל הפעם העבירו את dA לפני dB (הסדר ה"טבעי" והשגוי). אל תשנו כלום אחר. הריצו את checkResult שוב. מה קרה, והאם cuBLAS החזירה שגיאה?

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

תרגיל 2 - cuBLAS מול הקרנל שלכם והפער מהתקרה

עכשיו נמדוד: כמה מהר cuBLAS לעומת ה-kernel ה-tiled שכתבתם ב-3.5, וכמה רחוקים שנינו מתקרת ה-H100.

  1. הגדילו ל-N = 4096 (גדול מספיק שה-kernel, ולא הoverhead, ישלוט בזמן). הוסיפו מד-זמן מבוסס cudaEvent עם warmup אחד ו-20 חזרות ממוצעות (כמו ב-3.5 וב-5.3).
  2. מדדו את cublasSgemm ודווחו זמן, GFLOP/s ו-%-peak מול H100_FP32_PEAK_GFLOPS.
  3. הביאו את ה-kernel mmTiled מפרויקט 3.5 (העתיקו את הפונקציה ואת הlaunch). מדדו גם אותו, על אותה N ואותם נתונים.
  4. אמתו ש-שניהם עוברים checkResult מול הייחוס (ל-N קטן) לפני שאתם מדווחים תפוקה. אף פעם אל תדווחו GFLOP/s של kernel שלא עבר נכונות.
  5. בנו טבלה: לכל אחד מהשניים - זמן, GFLOP/s, %-peak, ואת היחס זמן_tiled / זמן_cuBLAS. כתבו משפט: איזה אחוז מהתקרה כל אחד תפס, ומהו הפער.

רמז: מדדו רק את זמן ה-kernel (בין cudaEventRecord לפני ואחרי הlaunch), בלי ה-cudaMemcpy, כדי לבודד את איכות ה-kernel. ה-warmup הכרחי כי ההרצה הראשונה של cuBLAS כוללת אתחול פנימי ובחירת kernel. צפו לכך ש-cuBLAS יגיע ל-80%+ מהתקרה בעוד ה-tiled יישאר סביב 25% - הפער הוא register tiling, double-buffering ו-Tensor Cores שעוד לא עשיתם.

תרגיל 3 - הדלקת Tensor Cores עם cublasGemmEx

תחליפו את מנוע החישוב מ-CUDA Cores ל-Tensor Cores ותמדדו את הזינוק, קודם ב-TF32 ואז ב-FP16.

  1. מסלול TF32: החליפו את הקריאה ל-cublasGemmEx עם קלט ופלט CUDA_R_32F, טיפוס חישוב CUBLAS_COMPUTE_32F_FAST_TF32 ואלגוריתם CUBLAS_GEMM_DEFAULT. שמרו על אותו תעלול החלפת אופרנדים ואותם ld. מדדו זמן ותפוקה, והפעם דווחו %-peak מול תקרת TF32 של H100 (~494700 GFLOP/s).
  2. אמתו נכונות מול הייחוס, אך הרפו את הסובלנות ל-1e-2 (TF32 חותך מנטיסה).
  3. מסלול FP16: הקצו מערכי __half ב-device (dAh, dBh, dCh). כתבו kernel זעיר או לולאת host שממירה float -> __half (עם __float2half). קראו ל-cublasGemmEx עם קלט/פלט CUDA_R_16F וטיפוס חישוב CUBLAS_COMPUTE_32F (צבירת FP32 מעל קלט FP16). מדדו, ודווחו %-peak מול תקרת FP16 (~989400 GFLOP/s).
  4. אמתו נכונות מול הייחוס, אך הרפו עוד את הסובלנות ל-2e-2 (FP16 מאבד עוד מנטיסה). המירו את dCh חזרה ל-float לפני ההשוואה.
  5. בנו טבלת סיכום של שלושת המסלולים (FP32, TF32, FP16): זמן, GFLOP/s, %-peak של אותו מסלול, ו-speedup מול ה-FP32 של תרגיל 2. כתבו משפט על הקשר בין ירידת הדיוק לעליית התפוקה.

רמז: CUBLAS_COMPUTE_32F_FAST_TF32 הוא כל ההבדל בין FP32 טהור ל-TF32 - אין צורך לשנות את הנתונים בזיכרון, רק את טיפוס החישוב. שימו לב שהם משווים מול תקרות שונות: אותה עבודה 2*N^3, אבל הגג זז. אם ה-speedup ב-FP16 קטן מפי 2, בדקו ש-N גדול מספיק (4096) ושאתם באמת פוגעים ב-Tensor Cores - kernel קטן מדי נשלט בlaunch overhead.

תרגיל 4 - לשבור את ה-leading dimension ולתקן

עכשיו נבין את מלכוד ה-ld מבפנים - קודם נשבור אותו בכוונה, נראה איך נראית התקלה, ואז נתקן.

  1. חזרו לקריאה התקינה של cublasSgemm מתרגיל 1 (עם checkResult שעובר). זה נקודת ההתחלה הבריאה.
  2. שבירה א (שגיאת ערך): שנו את ה-ld של האופרנד הראשון (dB) מ-N ל-N/2. הריצו. מכיוון שעטפתם ב-CUBLAS_CHECK, שימו לב איזה קוד סטטוס חוזר ומתי התוכנית נעצרת.
  3. שבירה ב (תוצאה שגויה בשקט): החזירו את ה-ld ל-N, אבל השתמשו במטריצות לא ריבועיות כדי לחשוף את הבאג האמיתי של ld. הגדירו M=256, K=512, N=128 (row-major: A בגודל M x K, B בגודל K x N, C בגודל M x N). קראו ל-cublasSgemm עם התעלול הנכון (m=N, n=M, k=K, האופרנד הראשון dB עם ld=N, השני dA עם ld=K, הפלט dC עם ld=N), ואמתו PASS מול ייחוס CPU מלבני.
  4. עכשיו קלקלו רק את ה-ld של dA מ-K ל-M (טעות נפוצה: לחשוב ש-ld הוא מספר השורות הלוגי במקום הצעד בזיכרון column-major). הריצו. האם קיבלתם שגיאה, תוצאה שגויה, או קריסה? הסבירו למה התוצאה שונה בין שבירה א לשבירה ב.
  5. תקנו חזרה לכל ה-ld הנכונים ואמתו PASS. סכמו: מהו בדיוק ה-ld שכל אופרנד דורש בתעלול החלפת האופרנדים, ולמה.

רמז: ld קטן מדי מ-max(1,m) הוא שגיאת ארגומנט שנתפסת מיד - CUBLAS_STATUS_INVALID_VALUE. ld שגוי אך "חוקי" גורם ל-cuBLAS לפרש את הצעד בין עמודות לא נכון, וכך היא קוראת מהמקום הלא נכון בזיכרון - שגיאה מספרית בשקט, בלי סטטוס שגיאה. בתעלול החלפת האופרנדים ל-C(M x N)=A(M x K)*B(K x N) row-major: dB דורש ld=N, dA דורש ld=K, dC דורש ld=N.

תרגיל 5 (בונוס) - לראות את ה-HMMA ב-SASS ואת בחירת ה-kernel

אמתו חזותית שה-Tensor Cores באמת נדלקו, ושהיוריסטיקת בחירת ה-kernel של cuBLAS עובדת.

  1. פרפלו את שלושת המסלולים (FP32, TF32, FP16) עם Nsight Systems: nsys profile -o gemm --stats=true ./cublas_gemm. אתרו בטבלת cuda_gpu_kern_sum את שמות ה-kernels שנבחרו - שימו לב שהם שונים בין המסלולים, כי cuBLAS בחרה kernel שונה לכל דיוק.
  2. פרפלו kernel יחיד עם Nsight Compute: ncu --set full -k regex:".*gemm.*" ./cublas_gemm. חפשו במטריקות את שיעור הניצול של ה-Tensor (למשל sm__pipe_tensor_op_hmma... או המדד "Tensor Core utilization") והשוו בין מסלול FP32 (שלא אמור לגעת ב-Tensor Cores) למסלול FP16 (שאמור להיות עמוס בהם).
  3. הריצו את מסלול ה-FP16 עם שני גדלי מטריצה שונים מאוד (למשל N=256 ו-N=8192) ובדקו ב-nsys האם cuBLAS בחרה אותו kernel לשניהם או kernels שונים. הסבירו את הממצא במונחי ההיוריסטיקה האטומה של cuBLAS.

רמז: את מנגנון ה-HMMA פירקנו ב-1.5; כאן אתם רואים אותו נבחר אוטומטית. ב-nsys stats השם של kernel של cuBLAS ארוך ומקודד (כולל את הטיפוס, את צורת הtile ואת דור ה-SM) - זה בדיוק ה"מאגר של עשרות kernels" מהשיעור. אם באותו דיוק אתם רואים kernel שונה לגדלים שונים, זו ההיוריסטיקה בפעולה.