לדלג לתוכן

2.3 כתיבת kernel אינדוקס ו launch configuration תרגול

תרגול - כתיבת kernel - אינדוקס ו-launch configuration

בתרגול הזה תכתבו את כפל המטריצות (matmul) הנאיבי, ה-kernel שילווה את כל הקורס. תריצו אותו, תאמתו אותו מול רפרנס CPU, תמדדו את מהירותו ותחשבו כמה GFLOP/s הוצאתם מהכרטיס. ואז - החלק החשוב - תוכיחו במספרים, מול נקודת המפנה (ridge point) של ה-H100, שה-kernel הזה memory-bound ולעולם לא ינצל את ה-GPU. לבסוף תהפכו אותו ל-kernel נייד עם לולאת grid-stride שעובד לכל N ולכל תצורת launch, ותתנסו בבחירת גודל ה-block. עבדו לפי הסדר - כל תרגיל נשען על הקודם. הריצו על GPU אמיתי; אם אין לכם, השתמשו ב-Google Colab (T4 חינמי) או במכונת ענן, כפי שראינו בשיעור 0.2. ה-H100 הוא כרטיס הייחוס, אך כל תרגיל ירוץ על כל GPU של NVIDIA (החליפו את -arch).

הכנה

צרו קובץ matmul.cu עם השלד הבא. הוא כולל את מקרו CUDA_CHECK, רפרנס CPU, את ה-kernel הנאיבי, ומדידת זמן עם cudaEvent. חלק מהשורות מסומנות ב-// TODO - תשלימו אותן בתרגילים.

#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <cuda_runtime.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)

// CPU reference - simple matrix multiplication, row-major
void matmulCPU(const float* A, const float* B, float* C, int N) {
    for (int row = 0; row < N; row++)
        for (int col = 0; col < N; col++) {
            float sum = 0.0f;
            for (int k = 0; k < N; k++)
                sum += A[row * N + k] * B[k * N + col];
            C[row * N + col] = sum;
        }
}

__global__ void mm(const float* A, const float* B, float* C, int N) {
    // TODO (exercise 1): compute row, col from the two axes, bounds guard, k loop
}

int main(void) {
    int N = 1024;
    size_t bytes = (size_t)N * N * sizeof(float);

    float *h_A = (float*)malloc(bytes);
    float *h_B = (float*)malloc(bytes);
    float *h_C = (float*)malloc(bytes);
    float *h_ref = (float*)malloc(bytes);
    for (int i = 0; i < N * N; i++) { h_A[i] = 1.0f; h_B[i] = 2.0f; }

    float *d_A, *d_B, *d_C;
    CUDA_CHECK(cudaMalloc(&d_A, bytes));
    CUDA_CHECK(cudaMalloc(&d_B, bytes));
    CUDA_CHECK(cudaMalloc(&d_C, bytes));
    CUDA_CHECK(cudaMemcpy(d_A, h_A, bytes, cudaMemcpyHostToDevice));
    CUDA_CHECK(cudaMemcpy(d_B, h_B, bytes, cudaMemcpyHostToDevice));

    dim3 block(16, 16);
    dim3 grid((N + block.x - 1) / block.x, (N + block.y - 1) / block.y);

    // warm-up
    mm<<<grid, block>>>(d_A, d_B, d_C, N);
    CUDA_CHECK(cudaGetLastError());
    CUDA_CHECK(cudaDeviceSynchronize());

    // measuring the kernel's time
    cudaEvent_t s, e;
    CUDA_CHECK(cudaEventCreate(&s));
    CUDA_CHECK(cudaEventCreate(&e));
    CUDA_CHECK(cudaEventRecord(s));
    mm<<<grid, block>>>(d_A, d_B, d_C, N);
    CUDA_CHECK(cudaEventRecord(e));
    CUDA_CHECK(cudaEventSynchronize(e));
    float ms = 0.0f;
    CUDA_CHECK(cudaEventElapsedTime(&ms, s, e));

    CUDA_CHECK(cudaMemcpy(h_C, d_C, bytes, cudaMemcpyDeviceToHost));

    // TODO (exercise 2): CPU reference, correctness check, and GFLOP/s calculation
    printf("N = %d, kernel = %.3f ms\n", N, ms);

    CUDA_CHECK(cudaFree(d_A)); CUDA_CHECK(cudaFree(d_B)); CUDA_CHECK(cudaFree(d_C));
    free(h_A); free(h_B); free(h_C); free(h_ref);
    return 0;
}

הקומפילציה לאורך התרגול (החליפו sm_90a בארכיטקטורה שלכם - למשל sm_75 ל-T4):

nvcc -O2 -arch=sm_90a -o matmul matmul.cu

לבירור ה-compute capability של הכרטיס:

nvidia-smi --query-gpu=name,compute_cap --format=csv

תרגיל 1 - מימוש ה-kernel הנאיבי

  1. השלימו את גוף ה-kernel mm: חשבו row מציר ה-y ו-col מציר ה-x לפי הנוסחה blockIdx.* * blockDim.* + threadIdx.*.
  2. הוסיפו שומר גבולות if (row < N && col < N) והריצו לולאת k שסוכמת A[row*N + k] * B[k*N + col] לתוך sum, ואז כתבו C[row*N + col] = sum.
  3. קמפלו והריצו. ודאו שהתוכנית לא קורסת ושמודפס זמן ה-kernel.
  4. הריצו תחילה עם N = 4 (מטריצה קטנה) והדפיסו את h_C איבר-איבר. עם A מלא באחדות ו-B מלא בשתיים, כל איבר של C צריך להיות N * (1*2) = 2N. בדקו שקיבלתם 2*4 = 8 בכל איבר.

רמז: זהו בדיוק ה-kernel מההרצאה. שימו לב שהעמודה col באה מציר ה-x וה-שורה row מציר ה-y - זו המוסכמה הטבעית כי ציר x של ה-grid הוא הממד ה"מהיר" (הרציף בזיכרון), ובמטריצה row-major העמודות הן הממד הרציף. אם קיבלתם transpose או תוצאות שגויות בשוליים, בדקו את השיטוח row*N + col ואת שומר הגבולות.

תרגיל 2 - אימות מול CPU וחישוב GFLOP/s

  1. חזירו את N ל-1024. השלימו את בלוק ה-TODO השני: קראו ל-matmulCPU(h_A, h_B, h_ref, N), ואז חשבו את השגיאה המקסימלית בין h_C ל-h_ref.
  2. הדפיסו את max error. עם A=1, B=2 ו-N=1024, כל איבר של C הוא 1024 * 2 = 2048; ודאו שהשגיאה זעירה (סדר גודל של 1e-3 או פחות - נקודה צפה מצטברת).
  3. חשבו את מספר ה-FLOPs של כפל מטריצות: 2 * N^3 (כל אחד מ-N^2 איברי הפלט דורש N כפלים ו-N חיבורים).
  4. מזמן ה-kernel שמדדתם (ms), חשבו את הביצועים: GFLOP/s = (2.0 * N^3) / (ms * 1e-3) / 1e9. הדפיסו את התוצאה.
  5. רשמו את ה-GFLOP/s שקיבלתם. השוו לשיא ה-FP32 של הכרטיס (ב-H100 כ-66,900 GFLOP/s) - איזה אחוז זעום הוצאתם?

רמז: 2 * N^3 עבור N=1024 הוא 2 * 1024^3 = 2.147e9 FLOPs. אם ה-kernel לקח כ-7 ms, זה בערך 2.147e9 / 7e-3 / 1e9 ≈ 307 GFLOP/s - סביב חצי אחוז מהשיא. אם השגיאה גדולה, ודאו שה-CPU וה-GPU קוראים את אותן h_A,h_B. הריצו את החימום לפני התזמון (כבר בשלד) כדי לא למדוד את ה-JIT הראשוני.

תרגיל 3 - arithmetic intensity ונקודת המפנה

זה התרגיל התיאורטי המרכזי. תוכיחו במספרים למה ה-kernel memory-bound.

  1. ספרו את התעבורה בלולאה הפנימית: כמה קריאות global וכמה FLOPs בכל איטרציה של k? מכאן חשבו את הarithmetic intensity ב-FLOP לכל בית (זכרו: float = 4 בתים).
  2. חשבו את נקודת המפנה (ridge point) של ה-H100: שיא FP32 / רוחב פס HBM3 = 66.9 TFLOPS / 3.35 TB/s. כמה FLOP-לבית יוצא?
  3. השוו: פי כמה הarithmetic intensity של ה-kernel נמוכה מנקודת המפנה? האם הוא memory-bound או compute-bound?
  4. חשבו את התקרה הmemory-bound: arithmetic intensity * רוחב פס = (FLOP/בית) * 3.35e12. כמה TFLOPS זו התקרה, ואיזה אחוז מהשיא?
  5. השוו את התקרה שחישבתם ל-GFLOP/s שמדדתם בתרגיל 2. האם המדידה מתחת לתקרה כצפוי?

רמז: 2 קריאות של 4 בתים = 8 בתים, מול 2 FLOPs, נותן 0.25 FLOP/בית. נקודת המפנה היא 66.9/3.35 ≈ 20 FLOP/בית, כלומר ה-kernel רחוק פי 80 מתחתיה - memory-bound באופן מובהק. התקרה: 0.25 * 3.35e12 = 0.84e12 = 0.84 TFLOPS, כ-1.25% מהשיא. המדידה שלכם (סביב 0.3 TFLOPS) אמורה לשבת מתחת ל-0.84 (caches קולטים חלק מהקריאות ומקרבים, אך גם coalescing לא מושלם וגישות עמודה ב-B מורידות). את מודל ה-Roofline המלא נפתח בפרק על ביצועים.

תרגיל 4 - המרה ללולאת grid-stride

  1. כתבו kernel חדש mmStride עם לולאת grid-stride דו-ממדית: כל thread מדלג בצעד blockDim.x*gridDim.x בעמודות ו-blockDim.y*gridDim.y בשורות (ראו את התבנית בהרצאה).
  2. הריצו אותו עם grid קטן במכוון: dim3 grid(4, 4) בלבד (16 בלוקים) עם block(16,16), על N = 1024. זה 64x64 = 4096 threads בלבד מול 1024*1024 ≈ מיליון איברי פלט.
  3. אמתו מול הרפרנס שהתוצאה עדיין נכונה למרות שיש הרבה פחות threads מאיברים.
  4. הריצו את אותו kernel גם עם N = 1000 (לא כפולה של 16) ו-N = 1500, וודאו נכונות בשני המקרים.
  5. הריצו גם עם block יחיד: <<<1, dim3(16,16)>>>. עדיין נכון (רק איטי)?

רמז: שני התנאים row < N ו-col < N בלולאות משמשים גם כשומרי גבולות - אין צורך ב-if נפרד. אם התוצאה שגויה עם ה-grid הקטן, בדקו שהצעד הוא += stride (ולא ++) ושחישבתם stride מ-gridDim ולא מקבוע. היכולת להריץ block יחיד ולקבל תוצאה נכונה היא בדיוק מה שהופך את ה-grid-stride לנייד ולקל-לבדיקה.

תרגיל 5 - השפעת גודל ה-block

  1. חזרו ל-kernel הנאיבי mm עם N = 2048. הריצו אותו עם ארבע תצורות block שונות ומדדו GFLOP/s לכל אחת: dim3(8,8), dim3(16,16), dim3(32,8), dim3(32,32).
  2. בנו טבלה: block, מספר ה-threads ב-block, ו-GFLOP/s.
  3. איזו תצורה הכי מהירה? איזו הכי איטית ולמה? (שימו לב ש-8x8=64 הוא שני warps בלבד, ו-32x32=1024 הוא block מקסימלי.)
  4. הסבירו את התוצאה במונחי occupancy וארגון ה-warps: ציר x בגודל כפולה של 32 עוזר ל-coalescing.

רמז: 8x8=64 threads בדרך כלל איטי כי הבלוקים קטנים מדי (מגבלת ה-32 בלוקים ל-SM חוסמת תפוסה) ו-8 בציר x אינו כפולה של 32, אז threads warp נחתכים על פני שורות. 16x16 ו-32x8 בדרך כלל הכי טובים. כל הבדל כאן קטן יחסית כי כולנו memory-bound - אף תצורה לא תציל kernel נאיבי, מה שרק מחזק את המסקנה מתרגיל 3.

תרגיל 6 (בונוס) - כמה רחוקים מהתקרה האמיתית

  1. קמפלו את matmul.cu והריצו את התוכנית תחת הפרופיילר ncu על ה-kernel הנאיבי, וקראו את המדד של ניצול רוחב הפס לזיכרון:
ncu --set basic --kernel-name mm ./matmul
  1. חפשו בפלט את Memory Throughput (אחוז) ואת Compute (SM) Throughput. איזה מהם גבוה - הזיכרון או החישוב? האם זה מאשר את הניתוח מתרגיל 3?
  2. אם יש לכם cuBLAS זמין, קמפלו תוכנית קטנה שקוראת ל-cublasSgemm על אותה בעיה (N=2048) ומדדו את ה-GFLOP/s שלה. פי כמה היא מהירה מה-kernel הנאיבי שלכם?
  3. הסבירו את הפער: cuBLAS משתמש ב-tiling דרך shared memory וב-Tensor Cores. את ה-kernel המרוצף נכתוב בעצמנו בפרק 3.3, ונשלים את הדרך עד Tensor Cores בפרויקט 3.5.

רמז: ב-kernel הנאיבי, Memory Throughput יהיה גבוה יחסית (עשרות אחוזים) בעוד Compute Throughput נמוך - חתימה קלאסית של kernel memory-bound. cuBLAS על H100 יגיע לעשרות אלפי GFLOP/s, כלומר פי 100 ומעלה מהנאיבי. שימו לב ש-cuBLAS הוא column-major - יש לשים לב לסדר הארגומנטים או לחשב B*A כדי לקבל את התוצאה הנכונה בפריסת row-major.