לדלג לתוכן

3.3 Shared memory ו tiling תרגול

תרגול - Shared memory ו-tiling

בתרגול הזה תהפכו את ה-matmul הנאיבי מ-2.3 ל-kernel מרוצף (tiled) מהיר, ותמדדו את השיפור במו עיניכם. תממשו את ה-matmul עם shared memory, תמדדו GFLOP/s ותשוו לנאיבי, תסרקו ערכי TILE_WIDTH שונים, תמירו את ה-kernel ל-shared memory דינמי, ולבסוף תזריקו ותתקנו את הבאג הקלאסי של מחסום __syncthreads() חסר. עבדו לפי הסדר - כל תרגיל נשען על הקודם. הריצו הכל על GPU אמיתי; אם אין לכם, השתמשו ב-Google Colab (T4 חינמי) או במכונת ענן, כפי שראינו בשיעור 0.2. ה-H100 הוא כרטיס הייחוס שלנו, אך כל תרגיל כאן ירוץ על כל GPU של NVIDIA (החליפו sm_90a ב-sm_75 ל-T4).

הכנה

צרו קובץ matmul.cu והתחילו מהשלד הבא, הכולל את מקרו CUDA_CHECK, את ה-kernel הנאיבי מ-2.3, בדיקת נכונות מול ה-CPU, ומדידת זמן עם cudaEvent. את ה-kernel המרוצף mmTiled תשלימו בתרגיל 1.

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

#define TILE_WIDTH 16

__global__ void mmNaive(const float* A, const float* B, float* C, int N) {
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    int col = blockIdx.x * blockDim.x + threadIdx.x;
    if (row < N && col < N) {
        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;
    }
}

// exercise 1: complete the tiled kernel here.
__global__ void mmTiled(const float* A, const float* B, float* C, int N) {
    // ...
}

// simple correctness check on the CPU (slow - use a small N for verification only)
void mmHost(const float* A, const float* B, float* C, int N) {
    for (int r = 0; r < N; r++)
        for (int c = 0; c < N; c++) {
            float s = 0.0f;
            for (int k = 0; k < N; k++) s += A[r * N + k] * B[k * N + c];
            C[r * N + c] = s;
        }
}

int main(int argc, char** argv) {
    int N = (argc > 1) ? atoi(argv[1]) : 1024;
    size_t bytes = (size_t)N * N * sizeof(float);

    float *h_A = (float*)malloc(bytes), *h_B = (float*)malloc(bytes);
    float *h_C = (float*)malloc(bytes), *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(TILE_WIDTH, TILE_WIDTH);
    dim3 grid((N + TILE_WIDTH - 1) / TILE_WIDTH,
              (N + TILE_WIDTH - 1) / TILE_WIDTH);

    cudaEvent_t start, stop;
    CUDA_CHECK(cudaEventCreate(&start));
    CUDA_CHECK(cudaEventCreate(&stop));

    // warm-up + measurement (replace mmNaive with mmTiled once you complete it)
    mmNaive<<<grid, block>>>(d_A, d_B, d_C, N);
    CUDA_CHECK(cudaDeviceSynchronize());

    CUDA_CHECK(cudaEventRecord(start));
    mmNaive<<<grid, block>>>(d_A, d_B, d_C, N);
    CUDA_CHECK(cudaEventRecord(stop));
    CUDA_CHECK(cudaEventSynchronize(stop));
    float ms = 0.0f;
    CUDA_CHECK(cudaEventElapsedTime(&ms, start, stop));

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

    double gflops = (2.0 * N * N * N) / (ms / 1e3) / 1e9;
    printf("N = %d, time = %.3f ms, %.1f GFLOP/s\n", N, ms, gflops);
    printf("C[0] = %.1f (expected %.1f)\n", h_C[0], 2.0f * N);

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

הקומפילציה לאורך כל התרגול:

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

שימו לב לנוסחת ה-GFLOP/s: matmul מרובע בגודל N מבצע 2*N^3 FLOPs (כפל וחיבור לכל אחד מ-N^3 הצירופים). למטריצות שכולן 1 ו-2, כל איבר של C שווה בדיוק 2*N, מה שנותן בדיקת שפיות מיידית.

תרגיל 1 - מימוש ה-matmul המרוצף ומדידת המהירות

  1. השלימו את mmTiled לפי ההרצאה: שני מערכי __shared__ בגודל TILE_WIDTH x TILE_WIDTH, לולאה על השלבים, טעינה שיתופית עם ריפוד אפסים בשוליים, שני מחסומי __syncthreads(), וכתיבה מוגנת ב-if (row < N && col < N).
  2. הריצו קודם עם N = 512 והשוו את C[0] שהתקבל לערך הצפוי 2*N. אם הוא נכון, המימוש תקין.
  3. החליפו את שתי הקריאות ל-mmNaive ב-main (החימום והמדידה) ל-mmTiled, וקמפלו מחדש.
  4. הריצו את שתי הגרסאות (הנאיבי והמרוצף) עבור N = 2048 ואז N = 4096, ורשמו את ה-GFLOP/s של כל אחת. חשבו את מקדם ההאצה (speedup).

רמז: הטעינה השיתופית היא As[ty][tx] = (row < N && aCol < N) ? A[row*N + aCol] : 0.0f; כאשר aCol = t*TILE_WIDTH + tx, ובאופן דומה Bs[ty][tx] עם bRow = t*TILE_WIDTH + ty. הלולאה הפנימית היא for (int k = 0; k < TILE_WIDTH; k++) sum += As[ty][k] * Bs[k][tx];. אם קיבלתם תוצאה שגויה רק בקצוות, כנראה שכחתם את ריפוד האפסים לtiles חלקיים.

תרגיל 2 - אימות נכונות מול ה-CPU

  1. הוסיפו קריאה ל-mmHost (הפניה שכבר בשלד) עבור N קטן, למשל N = 256, ומלאו את h_A, h_B בערכים אקראיים במקום קבועים (h_A[i] = (float)rand() / RAND_MAX;).
  2. השוו איבר-איבר את h_C (מה-GPU) ל-h_ref (מה-CPU) וחשבו את השגיאה המקסימלית: maxErr = fmax(maxErr, fabs(h_C[i] - h_ref[i]));.
  3. ודאו שהשגיאה זעירה (בסדר גודל של 1e-3 או פחות עבור float), לא אפס מוחלט.
  4. הסבירו לעצמכם למה השגיאה אינה אפס מדויק גם כשהאלגוריתם נכון.

רמז: ב-float (32 סיביות) הסדר שבו מצטברים המחצבים משפיע על תוצאת העיגול, וה-GPU צובר בסדר שונה מה-CPU. שגיאה יחסית קטנה (ולא אפס) היא הסימן לנכונות, לא לבאג. אל תשוו ב-== בין תוצאות נקודה-צפה - השוו סף (tolerance).

תרגיל 3 - סריקת TILE_WIDTH

  1. הפכו את TILE_WIDTH לפרמטר שקל לשנות (ערכו את ה-#define וקמפלו מחדש בכל ערך), וסרקו את הערכים 8, 16, 32.
  2. לכל ערך, הריצו עבור N = 4096 ורשמו את ה-GFLOP/s. בנו טבלה: TILE_WIDTH, מספר threads ב-block, בתי shared לכל block, GFLOP/s.
  3. חשבו לכל ערך את הarithmetic intensity התיאורטית (FLOPs לכל קריאת global) - היא אמורה לצאת TILE_WIDTH.
  4. איזה ערך היה הכי מהיר על הכרטיס שלכם? האם 32 תמיד מנצח את 16? הסבירו.

רמז: TILE_WIDTH = 8 נותן block של 64 threads (שני warps) - arithmetic intensity נמוכה ומעט ניצול-חוזר. TILE_WIDTH = 32 נותן block של 1024 threads (מקסימום) ו-8 KiB shared. לפעמים 32 מנצח ולפעמים 16, כי 32 צורך יותר משאבים ומוריד occupancy (ראו 8.1). בתי ה-shared הם 2 * TILE_WIDTH * TILE_WIDTH * 4.

תרגיל 4 - המרה ל-shared memory דינמי

  1. כתבו גרסה mmTiledDyn(const float* A, const float* B, float* C, int N, int tile) שבה גודל הtile מגיע כפרמטר בזמן ריצה, וה-shared memory מוכרז extern __shared__ float smem[];.
  2. חלקו את smem לשני חצאים: float* As = smem; ו-float* Bs = smem + tile*tile;. השתמשו באינדוקס משוטח As[ty*tile + tx] במקום As[ty][tx].
  3. השיקו עם הפרמטר השלישי בסוגריים המשולשים: mmTiledDyn<<<grid, block, 2*tile*tile*sizeof(float)>>>(d_A, d_B, d_C, N, tile); כאשר גם block בגודל tile x tile.
  4. ודאו שהתוצאה זהה לגרסה הסטטית, והשוו את ה-GFLOP/s. האם יש הבדל ביצועים ניכר?

רמז: הפרמטר השלישי של <<< >>> הוא מספר הבתים של ה-shared הדינמי, לא מספר האיברים - אל תשכחו את sizeof(float). אם קיבלתם invalid argument בlaunch, כנראה ביקשתם יותר shared ממה שמותר; לחריגה מ-48 KiB צריך cudaFuncSetAttribute(mmTiledDyn, cudaFuncAttributeMaxDynamicSharedMemorySize, bytes). הביצועים אמורים להיות כמעט זהים לסטטי - האינדוקס הידני אינו יקר יותר.

תרגיל 5 - הזרקת הבאג הקלאסי ותיקונו

  1. קחו עותק של mmTiled התקין וקראו לו mmTiledBug. מחקו את המחסום השני __syncthreads() (זה שאחרי הלולאה הפנימית), והשאירו את הראשון.
  2. הריצו עבור N = 1024, כמה פעמים ברצף. ודאו שאתם משתמשים באימות מול ה-CPU מתרגיל 2. תעדו: האם התוצאה שגויה? האם היא אותה שגיאה בכל הרצה?
  3. עכשיו מחקו במקום זאת את המחסום הראשון (והחזירו את השני). הריצו שוב ותעדו את ההתנהגות.
  4. החזירו את שני המחסומים ואמתו שהתוצאה שוב נכונה ויציבה.
  5. הסבירו במילים שלכם: איזו סכנה (hazard) כל מחסום מונע, ולמה השמטתו נותנת תוצאה לא דטרמיניסטית.

רמז: השמטת המחסום השני יוצרת סכנת כתיבה-לפני-קריאה (write-after-read): warp מהיר דורס את הtile בשלב הבא בעוד warp איטי עדיין קורא ממנו. התלות בתזמון ה-warps היא שהופכת את השגיאה ללא-יציבה בין הרצות. ל-N קטן שנכנס בtile יחיד (N <= TILE_WIDTH) יש רק שלב אחד, ואז הבאג עלול לא להתגלות - השתמשו ב-N גדול בהרבה מ-TILE_WIDTH כדי לכפות כמה שלבים.

תרגיל 6 (בונוס) - פרופיילינג ואחוז מהשיא

  1. הריצו את הגרסה המרוצפת תחת ncu והשוו לנאיבי, למשל:
    ncu --set full --kernel-name-base mangled -k "mm.*" ./matmul 2048
  2. השוו בין הנאיבי למרוצף את המדדים: DRAM Throughput, Compute (SM) Throughput, ו-Achieved Occupancy. ראו איך הנאיבי memory-bound (DRAM גבוה, SM נמוך) בעוד המרוצף מזיז את המשקל לכיוון החישוב.
  3. חפשו בדוח את מדד ה-Shared Memory Bank Conflicts. האם יש התנגשויות בגרסה שלכם?
  4. חשבו את אחוז השיא: קחו את ה-GFLOP/s הטוב ביותר שלכם וחלקו ב-66.9 TFLOPS (שיא ה-FP32 של H100, או השיא של הכרטיס שלכם). כמה רחוקים אתם מ-cuBLAS?

רמז: ncu דורש הרשאות (על ריצה מקומית ייתכן sudo או הגדרת NVreg_RestrictProfilingToAdminUsers=0). את הטיפול המלא בהתנגשויות בנקים ובריפוד [TILE_WIDTH + 1] נעשה בפרק 8.4, ואת האיחוד של כל שכבות האופטימיזציה (מרוצף -> ריצוף register -> Tensor Cores) - בפרויקט 3.5. אחוז של כמה אחוזים בודדים מהשיא הוא צפוי ל-kernel פשוט על CUDA Cores; cuBLAS מגיע לעשרות אחוזים כי הוא משתמש ב-Tensor Cores וריצוף רב-שכבתי.