3.5 פרויקט כפל מטריצות tiled פרויקט
פרויקט - כפל מטריצות tiled¶
זהו פרויקט הביניים המסכם של הקורס עד כה. בפרקים 0 עד 3 אספנו את כל אבני הבניין: כתבנו kernel ראשון והבנו את מודל ה-host/device (פרק 0), פירקנו את היררכיית ה-threads, ה-blocks וה-warps (פרק 2), ולמדנו את היררכיית הזיכרון - registers, shared memory ו-global memory - ואת ההבדל הביצועי העצום ביניהם (פרק 3). עכשיו נחבר הכל למשימה אחת קלאסית שהיא גם ה-"Drosophila" של תכנות GPU: כפל מטריצות (matrix multiply, בקיצור matmul או SGEMM עבור single-precision). נבנה אותו בשלבים - naive, ואז coalesced, ואז tiled עם shared memory - נמדוד כל שלב במדדים אמיתיים (זמן, GFLOP/s, אחוז מה-peak), ונשווה לתקרה: הספרייה המקצועית cuBLAS. המטרה איננה לנצח את cuBLAS (אי אפשר, וזה בסדר), אלא להסביר כל קפיצת ביצועים במושגים שלמדנו: הarithmetic intensity, coalescing (איחוד גישות), וניצול-חוזר (reuse) של נתונים דרך shared memory. את שני השלבים המתקדמים שנשארים - register tiling ו-Tensor Cores - נבנה בפרויקט 9.3, שם נטפס את שאר הסולם עד קרוב ל-peak האמיתי.
מה נבנה - the GEMM ladder¶
הפרויקט בנוי כסולם (ladder) של גרסאות ל-kernel, כל אחת מתקנת את צוואר הבקבוק של קודמתה. בכל שלב מודדים, רושמים בטבלה, ומסבירים. הנה מפת הדרך:
stage kernel what it fixes bottleneck after
─────────────────────────────────────────────────────────────────────────
stage 0 naive matmul correctness vs CPU global-memory bound,
(baseline) non-coalesced accesses
stage 1 coalesced thread->data mapping still memory-bound,
so accesses are coalesced every FLOP reads from global
stage 2 tiled (shared mem) staging through shared mem tile-overhead bound
-> reuse, less traffic (one tile per thread)
─────────────────────────────────────────────────────────────────────────
ceiling cuBLAS SGEMM reference ceiling (professional; target 9.3)
─────────────────────────────────────────────────────────────────────────
[9.3] register tiling + each thread computes an output tile, Tensor Cores
Tensor Cores then matmuls on the Tensor Core
מטריצות ההרצה יהיו ריבועיות בגודל N x N מסוג float (FP32). נעבוד עם N = 4096 כערך ברירת המחדל למדידה (גדול מספיק כדי שה-kernel, ולא הoverhead, ישלוט בזמן), אך הקוד יעבוד לכל N. ה-GPU הרץ לדוגמה הוא H100 (SXM), כרגיל בקורס.
הרקע המספרי - the metrics we report¶
לפני שכותבים שורת קוד, צריך להסכים איך מודדים. שלושה מספרים ילוו כל שלב.
1. מספר פעולות הנקודה הצפה (FLOPs). כפל של שתי מטריצות N x N דורש, לכל אחד מ-N^2 איברי הפלט, N כפלים ו-N חיבורים - כלומר 2N פעולות. סך הכל:
עבור N = 4096 יוצא 2 * 4096^3 = 2 * 68,719,476,736 ≈ 1.374 x 10^11 פעולות, כלומר בערך 137.4 GFLOP של עבודה חישובית שקבועה בכל השלבים. מה שמשתנה בין השלבים הוא רק כמה מהר אנחנו מבצעים את אותה כמות עבודה.
2. תפוקה - GFLOP/s. מחלקים את העבודה בזמן:
3. אחוז מה-peak. ה-H100 מגדיר תקרה תיאורטית לחישוב FP32 שאינו על Tensor Cores. נגזור אותה: 132 SMs, ובכל SM 128 ליבות CUDA (CUDA Cores) ל-FP32, כל אחת מבצעת פעולת FFMA (fused multiply-add) שהיא 2 FLOP למחזור שעון, בתדר boost של כ-1.98 GHz:
כל התוצאות שלנו יימדדו כאחוז מ-66,900 GFLOP/s. שימו לב לנקודה חינוכית: cuBLAS ב-FP32 מגיע בדרך כלל ל-80% ומעלה מהתקרה הזו, בעוד ה-kernels שנכתוב בעצמנו יישארו נמוך בהרבה - וזה בדיוק העניין. הפער הוא מה שנסביר.
הarithmetic intensity. המושג המרכזי מפרק 3: היחס בין פעולות חישוב לבין בתים שנקראים מ-global memory. ה-kernel ה-naive מבצע בערך FLOP אחד לכל קריאה מ-global memory (שני loads ל-FMA שמייצר ~2 FLOP). כפי שהזהיר אותנו הגלוסארי: "לעולם לא תנצלו את כל ה-GPU כך, כי רוחב הפס האריתמטי של ה-CUDA Cores גבוה בהרבה מרוחב הפס של הזיכרון בין ה-GPU RAM ל-SMs". כל הסולם הוא, במובן אחד, מסע להעלאת הarithmetic intensity.
תשתית ההתחלה - starter scaffolding¶
הקוד הבא הוא השלד שתתחילו ממנו. הוא כולל את מקרו ה-CUDA_CHECK, מימוש ייחוס (reference) על ה-CPU, בודק נכונות עם סובלנות (tolerance), מד-זמן מבוסס cudaEvent, ומחשבון GFLOP/s ו-%-peak. שמרו אותו כ-matmul.cu והשלימו את ה-kernels שלב אחר שלב.
#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)
// the theoretical peak of H100 for non-Tensor-Core FP32, in GFLOP/s
static const double H100_FP32_PEAK_GFLOPS = 66900.0;
// ------- reference implementation on the CPU (row-major): C = A * B -------
// intentionally slow; run it once for a small N to verify correctness.
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 sum = 0.0; // accumulate in double for reference precision
for (int k = 0; k < N; k++)
sum += (double)A[r * N + k] * (double)B[k * N + c];
C[r * N + c] = (float)sum;
}
}
// ------- correctness check with relative tolerance -------
// matrix multiply accumulates rounding error; never compare exact equality.
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 r = ref[i], g = got[i];
double denom = fmax(fabs(r), 1e-6);
double rel = fabs(g - r) / denom;
if (rel > maxRel) maxRel = rel;
}
printf(" max relative error = %.3e (tolerance %.1e)\n", maxRel, rtol);
return maxRel <= rtol;
}
// ------- metrics calculator -------
void report(const char* name, int N, float ms) {
double flops = 2.0 * (double)N * (double)N * (double)N;
double gflops = flops / (ms / 1000.0) / 1e9;
double pct = 100.0 * gflops / H100_FP32_PEAK_GFLOPS;
printf(" %-16s N=%d %8.3f ms %9.1f GFLOP/s %5.2f%% peak\n",
name, N, ms, gflops, pct);
}
// ------- cudaEvent-based timer -------
// run one warmup (not measured), then REPS runs, and return the average time per kernel.
template <typename LaunchFn>
float timeKernel(LaunchFn launch, int reps = 20) {
cudaEvent_t start, stop;
CUDA_CHECK(cudaEventCreate(&start));
CUDA_CHECK(cudaEventCreate(&stop));
launch(); // warmup: takes JIT/caches out of the picture
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaEventRecord(start));
for (int r = 0; r < reps; r++) launch();
CUDA_CHECK(cudaEventRecord(stop));
CUDA_CHECK(cudaEventSynchronize(stop)); // blocks until the GPU has finished
float ms = 0.0f;
CUDA_CHECK(cudaEventElapsedTime(&ms, start, stop));
CUDA_CHECK(cudaEventDestroy(start));
CUDA_CHECK(cudaEventDestroy(stop));
return ms / reps; // average time per run
}
שלוש הערות על התשתית שחשוב להפנים:
- מודדים את זמן ה-kernel בלבד, לא את ההעברות. ה-
cudaEventמתעד נקודות זמן על ה-stream של ה-GPU, ולכן מודד את זמן ה-kernel נטו, בלי ה-cudaMemcpyשל הקלט והפלט. זו המדידה הנכונה כדי לבודד את איכות ה-kernel. כפי שראינו בשיעור 0.3, אם היינו כוללים את ה-PCIe, היינו מודדים את ה-bus ולא את החישוב. - warmup הכרחי. ההרצה הראשונה של kernel כוללת JIT של PTX ל-SASS, אכלוס caches וקיבוע תדר. מריצים אותה פעם אחת מחוץ למדידה, אחרת ההרצה הראשונה מזהמת את הממוצע.
- צוברים ב-double בייחוס ובודקים בסובלנות יחסית. ה-kernel על ה-GPU צובר ב-FP32, ולכן ייתן תוצאה מעט שונה מ-CPU שצובר ב-double. השוואת שוויון מדויק תיכשל תמיד; לכן
rtolשל כ-1e-3היא הנכונה עבור SGEMM.
שלב 0 - the naive baseline¶
נתחיל מה-kernel ה-naive שכתבנו בשיעור 2.3: כל thread אחראי לאיבר פלט אחד C[row][col], ומחשב אותו בלולאה על k. המיפוי ה"אינטואיטיבי" הראשון שרובנו כותבים משייך את ציר ה-x של ה-thread לשורה:
__global__ void mmNaive(const float* A, const float* B, float* C, int N) {
int row = blockIdx.x * blockDim.x + threadIdx.x; // x -> row
int col = blockIdx.y * blockDim.y + threadIdx.y; // y -> column
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;
}
}
הlaunch, עם blocks דו-ממדיים בגודל 16 x 16 = 256 threads:
dim3 block(16, 16);
dim3 grid((N + 15) / 16, (N + 15) / 16);
auto launch = [&]{ mmNaive<<<grid, block>>>(dA, dB, dC, N); };
float ms = timeKernel(launch);
report("naive", N, ms);
קריטריון קבלה לשלב 0: ה-kernel עובר את checkResult מול matmulCPU (ל-N קטן, למשל 512, כדי שהייחוס על ה-CPU יסתיים בזמן סביר), ואתם מודדים ורושמים GFLOP/s ו-%-peak עבור N = 4096. אין דרישת ביצועים - זהו בסיס ההשוואה בלבד.
מספרים מייצגים על H100 (ישתנו אצלכם, אך סדר הגודל יציב):
בערך 3% מה-peak. למה כה נמוך? המיפוי x -> row הוא האסון: בתוך warp (קבוצת 32 threads) ה-threadIdx.x רץ, ולכן threads סמוכים מקבלים שורות סמוכות ואותה עמודה. הכתיבה C[row * N + col] של threads סמוכים קופצת בצעדים של N איברים בזיכרון - גישה לא-מאוחדת (uncoalesced) - וכך גם הקריאה A[row * N + k]. כל transaction של הזיכרון מביא 32 בתים שמתוכם משתמשים ב-4 בלבד. את זה מתקנים בשלב 1.
שלב 1 - global memory coalescing¶
התיקון קטן להפליא בקוד ועצום בביצועים: מחליפים את המיפוי כך שציר ה-x יזין את העמודה. עכשיו threads סמוכים (ב-threadIdx.x) מקבלים עמודות סמוכות ואותה שורה:
__global__ void mmCoalesced(const float* A, const float* B, float* C, int N) {
int col = blockIdx.x * blockDim.x + threadIdx.x; // x -> column
int row = blockIdx.y * blockDim.y + threadIdx.y; // y -> row
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;
}
}
עכשיו נתבונן בגישות של warp שלם (32 threads עם col רץ, row קבוע), לכל ערך k בלולאה:
A[row*N + k] -> same address for all 32 threads = broadcast (efficient)
B[k*N + col] -> 32 contiguous addresses = coalesced (one transaction)
C[row*N+col] -> 32 contiguous addresses = coalesced (one transaction)
שלוש הגישות אידיאליות: אחת broadcast ושתיים מאוחדות. הקוד זהה כמעט לחלוטין ל-naive - רק החלפנו איזה ציר מזין מה - אך כל transaction של הזיכרון מנוצל עכשיו במלואו.
קריטריון קבלה לשלב 1: אותו checkResult עובר, ואתם מראים שיפור של כמה-מונים ב-GFLOP/s לעומת שלב 0. רשמו את היחס זמן_naive / זמן_coalesced.
פי 3.4 מהר יותר, מתיקון של מיפוי אינדקסים בלבד, בלי לגעת בכמות החישוב. זו ההדגמה הנקייה ביותר בקורס לכך שcoalescing הוא לא אופטימיזציה מתקדמת - הוא תנאי סף. ובכל זאת, 10% מה-peak עדיין רחוק מהתקרה. הסיבה: ה-kernel נשאר memory-bound. לכל FMA אנחנו עדיין קוראים מ-global memory; הarithmetic intensity לא השתנתה, רק ניצול רוחב הפס. את כמות התעבורה נחתוך בשלב 2.
שלב 2 - shared memory tiling¶
זה השלב שבו מיישמים את הרעיון המרכזי של פרק 3, ואת התבנית שהראה הגלוסארי: staging דרך shared memory (זיכרון משותף) כדי למקסם את מספר הפעולות האריתמטיות לכל טעינה מ-global memory. במקום שכל thread יקרא בעצמו את כל שורת A ועמודת B, ה-block משתף פעולה: הוא טוען tile של A וtile של B ל-shared memory, כל הthreads משתמשים בהם שוב ושוב, ואז עוברים לtile הבא.
#define TILE 32 // 32x32 tile = 1024 threads per block
__global__ void mmTiled(const float* A, const float* B, float* C, int N) {
__shared__ float As[TILE][TILE]; // tile of A in shared memory
__shared__ float Bs[TILE][TILE]; // tile of B in shared memory
int tx = threadIdx.x, ty = threadIdx.y;
int row = blockIdx.y * TILE + ty; // output row for this thread
int col = blockIdx.x * TILE + tx; // output column for this thread
float sum = 0.0f;
int numTiles = (N + TILE - 1) / TILE;
for (int t = 0; t < numTiles; t++) {
int aCol = t * TILE + tx; // coalesced load of tile A
int bRow = t * TILE + ty; // coalesced load of tile B
As[ty][tx] = (row < N && aCol < N) ? A[row * N + aCol] : 0.0f;
Bs[ty][tx] = (bRow < N && col < N) ? B[bRow * N + col ] : 0.0f;
__syncthreads(); // barrier: the tile is fully loaded
for (int k = 0; k < TILE; k++) // every FMA reads from shared, not global
sum += As[ty][k] * Bs[k][tx];
__syncthreads(); // barrier: before overwriting the tile
}
if (row < N && col < N) C[row * N + col] = sum;
}
הlaunch עם blocks בגודל 32 x 32:
dim3 block(TILE, TILE); // 1024 threads - the maximum per block
dim3 grid((N + TILE - 1) / TILE, (N + TILE - 1) / TILE);
auto launch = [&]{ mmTiled<<<grid, block>>>(dA, dB, dC, N); };
מדוע זה עובד - החשבון של הarithmetic intensity. בלי tiling, כל FMA דרש load מ-global memory. עם tile בגודל TILE = 32, טוענים כל איבר של A ושל B ל-shared memory פעם אחת, ואז משתמשים בו ב-TILE = 32 פעולות FMA לפני שנזרק. חתכנו את תעבורת ה-global memory פי TILE. שני המחסומים __syncthreads() הכרחיים: הראשון מבטיח שהtile נטען לגמרי לפני שמישהו קורא ממנו, והשני מבטיח שכולם סיימו לחשב לפני שהtile הבא דורס אותו.
staging of one tile (TILE=32):
global loads: 2 * 32 * 32 = 2048 elements (tile A + tile B)
FMAs produced: 32 * 32 * 32 = 32768 operations
-> ~32 FLOP per load from global, vs ~1 in naive
קריטריון קבלה לשלב 2: checkResult עובר, ואתם מראים שיפור של כמה-מונים לעומת שלב 1 (coalesced). מדדו גם עם TILE = 16 ועם TILE = 32 והשוו.
פי 2.5 מהר יותר מ-coalesced, ופי 8.5 מ-naive. הגענו ל-26% מה-peak - עדיין לא התקרה, כי כל thread מחשב עדיין רק איבר פלט אחד, וה-kernel חסום בoverhead של טעינת tiles ומחסומים יחסית לחישוב. את הצעד הבא - שכל thread יחשב tile שלם של פלט מתוך ה-registers שלו (register tiling) - נעשה ב-9.3.
שלב 3 - the reference ceiling: cuBLAS¶
עכשיו נמדוד את התקרה. cuBLAS היא ספריית ה-BLAS המקצועית של NVIDIA, כתובה ומכווננת ידנית ברמת ה-SASS לכל ארכיטקטורה. הפונקציה cublasSgemm היא SGEMM - Single-precision GEneral Matrix Multiply. נשווה אליה כדי לדעת כמה "בשר" נשאר על העצם.
יש כאן מלכוד שחייבים לזכור: cuBLAS היא column-major (עמודה-ראשית), בעוד המטריצות שלנו ב-C הן row-major (שורה-ראשית). מטריצה row-major M שנקראת column-major היא בעצם M^T. הטריק הסטנדרטי: כדי לחשב C = A * B בשמירה row-major, מבקשים מ-cuBLAS לחשב, בעולם ה-column-major שלה, את C^T = B^T * A^T - כלומר מעבירים את B ו-A בסדר הפוך ובלי transpose:
#include <cublas_v2.h>
cublasHandle_t handle;
cublasCreate(&handle);
const float alpha = 1.0f, beta = 0.0f;
// computes (in the column-major world) dC = dB * dA, which is (A*B)^T,
// and when read back row-major you get exactly A*B.
auto launch = [&]{
cublasSgemm(handle,
CUBLAS_OP_N, CUBLAS_OP_N,
N, N, N, // m, n, k
&alpha,
dB, N, // cuBLAS's "A" = our B
dA, N, // cuBLAS's "B" = our A
&beta,
dC, N);
};
float ms = timeKernel(launch);
report("cuBLAS", N, ms);
cublasDestroy(handle);
מקמפלים עם קישור לספרייה:
קריטריון קבלה לשלב 3: checkResult על פלט cuBLAS עובר (הוכחה שהטריק ה-column-major נכון), ואתם מדווחים את הפער: זמן_tiled / זמן_cuBLAS ואת אחוז ה-peak של cuBLAS.
cuBLAS מגיע ל-86% מה-peak - פי 3.3 מה-kernel ה-tiled שלנו, ופי 28 מה-naive. זו התקרה שאליה נטפס בפרויקט 9.3 עם register tiling ו-Tensor Cores. שימו לב: cuBLAS מקבל את התוצאה הזו ב-FP32 טהור; אם היינו מרשים TF32 על Tensor Cores הוא היה עובר אף את מספר ה-66.9 TFLOPS הזה, כי אז מודדים מול peak אחר לגמרי - נושא לפרק 6.
טבלת התוצאות - the results table¶
הדליוורבל המרכזי של הפרויקט הוא הטבלה הבאה, ממולאת מהמדידות שלכם (עם N = 4096). הערכים כאן הם דוגמה מ-H100 - שלכם ישתנו:
| N | kernel | ms | GFLOP/s | % peak | speedup מול naive |
|---|---|---|---|---|---|
| 4096 | naive | 68.000 | 2021 | 3.0% | 1.0x |
| 4096 | coalesced | 20.000 | 6872 | 10.3% | 3.4x |
| 4096 | tiled (TILE=32) | 8.000 | 17180 | 25.7% | 8.5x |
| 4096 | cuBLAS | 2.400 | 57267 | 85.6% | 28.3x |
בנוסף לטבלה, כתבו פסקה קצרה לכל מעבר בין שלבים המסבירה את השיפור במושג שלמדנו:
- naive -> coalesced: לא שינינו את כמות העבודה ולא את הarithmetic intensity - רק תיקנו את מיפוי ה-thread לנתון כך שגישות ה-warp ל-global memory מאוחדות (coalesced). כל transaction של הזיכרון מנוצל במלואו במקום שמינית ממנו.
- coalesced -> tiled: כאן שינינו את הarithmetic intensity. staging דרך shared memory נותן ניצול-חוזר (reuse): כל איבר נטען מ-global memory פעם אחת ומשמש ב-
TILEפעולות FMA. חתכנו את תעבורת ה-global memory פיTILE, והזזנו את צוואר הבקבוק ממנה הלאה. - tiled -> cuBLAS: הפער שנשאר נובע מאופטימיזציות שעוד לא עשינו - register tiling (כל thread מחשב tile פלט מה-registers), double-buffering של הטעינות, וניצול Tensor Cores. אלה בדיוק שלבי 9.3.
קריטריוני קבלה - acceptance criteria¶
הפרויקט נחשב מושלם כאשר כל השורות הבאות מתקיימות:
[ ] stage 0: mmNaive passes checkResult (rtol=1e-3) vs matmulCPU for N=512
[ ] stage 0: GFLOP/s and %peak measurement for N=4096 recorded in the table
[ ] stage 1: mmCoalesced passes checkResult, and is significantly faster than stage 0
[ ] stage 2: mmTiled passes checkResult, and is several times faster than stage 1
[ ] stage 2: comparison of TILE=16 vs TILE=32 documented
[ ] stage 3: cuBLAS passes checkResult (proof the column-major trick is correct)
[ ] the table (N, kernel, ms, GFLOP/s, %peak) is complete for all four stages
[ ] explanation paragraph for each transition, in terms of coalescing / arithmetic intensity / reuse
[ ] all CUDA calls wrapped in CUDA_CHECK
הערה על נכונות מול ביצועים: אף פעם אל תדווחו GFLOP/s של kernel שלא עבר את בודק הנכונות. kernel מהיר ושגוי הוא חסר ערך, וטעות נפוצה היא "לשפר" ביצועים תוך שבירת הנכונות בשקט (למשל __syncthreads חסר בשלב 2 יאיץ מעט ויחזיר זבל לא-דטרמיניסטי). הריצו את checkResult אחרי כל שינוי.
לאן ממשיכים - what comes next¶
הגענו מ-3% ל-26% מה-peak בשלושה שלבים, וראינו שהתקרה (cuBLAS) עומדת על 86%. הפער שנשאר איננו קסם - הוא סדרה של אופטימיזציות מוגדרות היטב שנבנה בפרויקט 9.3, ה-GEMM המלא:
- register tiling: כל thread מחשב tile פלט של
8 x 8(למשל) מתוך ה-registers שלו, מה שמכפיל שוב את הarithmetic intensity ומקטין את הoverhead של המחסומים לכל פעולה. - double-buffering: טעינת הtile הבא ל-shared memory במקביל לחישוב על הtile הנוכחי, כדי להסתיר את ה-latency של global memory.
- Tensor Cores: המעבר מ-CUDA Cores ל-Tensor Cores (ליבות טנזור), שרוחב הפס האריתמטי שלהן גבוה בסדרי גודל. כפי שהעיר הגלוסארי, ה-matmul המהירים ביותר רצים על Tensor Cores, לא על CUDA Cores - וזו הסיבה שאפילו ה-kernel ה-tiled שלנו רחוק מהתקרה. את התכנות שלהן, ואת מודל ה-tile-based שנבנה סביבן, נפגוש בפרק 6 ובפרויקט 9.3.
הפרויקט הזה הוא ההוכחה שהמושגים של פרקים 0 עד 3 - kernels, אינדוקס, coalescing, ו-shared memory - הם לא תיאוריה: כל אחד מהם שווה קפיצת ביצועים מדידה, וכולם יחד לוקחים אותנו כמעט פי תשעה מהנקודה שממנה התחלנו.