3.3 Shared memory ו tiling פתרון
פתרון - Shared memory ו-tiling¶
הפתרונות כאן הם דוגמה מייצגת. המספרים המדויקים (GFLOP/s, זמנים, אחוזי occupancy) משתנים בין דורות חומרה, גרסאות דרייבר ו-nvcc, ואפילו בין הרצות; כל המדידות למטה נלקחו על H100 SXM כסדר גודל. מה שחשוב הוא היחסים - הקפיצה מהנאיבי למרוצף, והתלות ב-TILE_WIDTH - לא הספרה האחרונה.
פתרון תרגיל 1 - מימוש ה-matmul המרוצף ומדידת המהירות¶
ה-kernel המרוצף המלא:
__global__ void mmTiled(const float* A, const float* B, float* C, int N) {
__shared__ float As[TILE_WIDTH][TILE_WIDTH];
__shared__ float Bs[TILE_WIDTH][TILE_WIDTH];
int tx = threadIdx.x, ty = threadIdx.y;
int row = blockIdx.y * TILE_WIDTH + ty;
int col = blockIdx.x * TILE_WIDTH + tx;
float sum = 0.0f;
int numTiles = (N + TILE_WIDTH - 1) / TILE_WIDTH;
for (int t = 0; t < numTiles; t++) {
int aCol = t * TILE_WIDTH + tx;
int bRow = t * TILE_WIDTH + ty;
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 1
for (int k = 0; k < TILE_WIDTH; k++)
sum += As[ty][k] * Bs[k][tx];
__syncthreads(); // barrier 2
}
if (row < N && col < N)
C[row * N + col] = sum;
}
לאחר החלפת mmNaive ב-mmTiled בשתי הקריאות שב-main, הפלט להרצות N = 2048 ו-N = 4096:
# naive, N = 4096
N = 4096, time = 458.210 ms, 300.0 GFLOP/s
C[0] = 8192.0 (expected 8192.0)
# tiled (TILE_WIDTH = 16), N = 4096
N = 4096, time = 45.900 ms, 2994.4 GFLOP/s
C[0] = 8192.0 (expected 8192.0)
מקדם ההאצה כאן הוא כ-458 / 46 ≈ 10. הבדיקה C[0] = 8192 (שהוא 2*N = 2*4096) מאשרת נכונות.
למה זה עבד: ה-tiling הפך את דפוס הגישה. במקום שכל thread יקרא את שורת A ועמודת B שלו מ-global memory, ה-block טוען כל tile פעם אחת ל-shared memory, ואז כל איבר טעון משמש TILE_WIDTH threads. הarithmetic intensity עלתה מ-1 ל-16 FLOPs לכל קריאת global, וה-kernel חדל להיות מורעב-זיכרון. שני מחסומי ה-__syncthreads() הם שמבטיחים שהtile מלא לפני הקריאה ושלא נדרס לפני שסיימו איתו.
איך להכליל: אותה תבנית staging (טעינה שיתופית -> מחסום -> חישוב מ-shared -> מחסום) עובדת בכל kernel memory-bound עם ניצול-חוזר של נתונים: קונבולוציה (stencil), היסטוגרמות, reductions, וסריקות. כל אימת שאותו נתון global נקרא בידי כמה threads של אותו block, שקלו לרצף אותו דרך shared memory.
פתרון תרגיל 2 - אימות נכונות מול ה-CPU¶
הוספת אתחול אקראי ואימות:
srand(1234);
for (int i = 0; i < N * N; i++) {
h_A[i] = (float)rand() / RAND_MAX;
h_B[i] = (float)rand() / RAND_MAX;
}
// ... after copying back:
mmHost(h_A, h_B, h_ref, N); // for small N only (256)
double maxErr = 0.0;
for (int i = 0; i < N * N; i++)
maxErr = fmax(maxErr, fabs((double)h_C[i] - (double)h_ref[i]));
printf("max abs error vs CPU = %.3e\n", maxErr);
הפלט האופייני עבור N = 256:
השגיאה אינה אפס, אך היא זעירה - בסדר גודל של 1e-4, כצפוי מ-float.
למה זה עבד: ה-float מדויק ל-24 סיביות מנטיסה בלבד. כשצוברים N = 256 מחצבים, כל חיבור מעגל, וסדר החיבור קובע את שגיאת העיגול המצטברת. ה-GPU צובר בסדר שונה מלולאת ה-CPU (ובtiles), ולכן התוצאות מתפצלות בספרות הנמוכות. שגיאה יחסית קטנה, ולא אפס מדויק, היא הסימן לנכונות בחישוב נקודה-צפה.
איך להכליל: לעולם אל תשוו תוצאות נקודה-צפה ב-==. הגדירו סף (tolerance) יחסי או מוחלט המתאים לטיפוס ולגודל הבעיה. לדיוק גבוה יותר עברו ל-FP64 (במחיר תפוקה) או השתמשו בצבירה בעלת דיוק מעורב (mixed precision) עם צובר FP32 מעל מחצבי FP16 - בדיוק מה ש-Tensor Cores עושים.
פתרון תרגיל 3 - סריקת TILE_WIDTH¶
הטבלה האופיינית עבור N = 4096 על H100:
TILE_WIDTH | threads/block | shared/block | GFLOP/s | FLOPs/read
-----------+---------------+--------------+---------+-----------
8 | 64 | 1024 B | 1650 | 8
16 | 256 | 2048 B | 2994 | 16
32 | 1024 | 8192 B | 4480 | 32
למה זה עבד: הarithmetic intensity שווה בדיוק ל-TILE_WIDTH (8, 16, 32 FLOPs לכל קריאת global), ולכן tile גדול יותר קורא פחות מ-global memory לכל FLOP. לכן 32 מנצח כאן את 16. אבל זה לא חוק ברזל: TILE_WIDTH = 8 נותן רק שני warps ל-block וarithmetic intensity נמוכה, ולכן הוא האיטי ביותר; 32 מגיע לגבול של 1024 threads ל-block ו-8 KiB shared, מה שמוריד את מספר ה-blocks התושבים ל-SM (רק 2 לפי מכסת ה-threads) ועלול, ב-kernel אחר, לפגוע ב-occupancy עד כדי האטה.
איך להכליל: אין TILE_WIDTH אופטימלי יחיד - הוא טרייד-אוף בין ניצול-חוזר (רוצים tile גדול) לבין occupancy וצריכת משאבים (שמעדיפים tile קטן). תמיד סורקים את הפרמטר על החומרה היעד ומודדים. אותו עיקרון חל על כל פרמטר כוונון (blockDim, register tiles, מקדם unrolling): קבעו אותו אמפירית, לא לפי אינטואיציה.
פתרון תרגיל 4 - המרה ל-shared memory דינמי¶
הגרסה הדינמית:
__global__ void mmTiledDyn(const float* A, const float* B, float* C,
int N, int tile) {
extern __shared__ float smem[];
float* As = smem;
float* Bs = smem + tile * tile;
int tx = threadIdx.x, ty = threadIdx.y;
int row = blockIdx.y * tile + ty;
int col = blockIdx.x * tile + tx;
float sum = 0.0f;
int numTiles = (N + tile - 1) / tile;
for (int t = 0; t < numTiles; t++) {
int aCol = t * tile + tx;
int bRow = t * tile + ty;
As[ty * tile + tx] = (row < N && aCol < N) ? A[row * N + aCol] : 0.0f;
Bs[ty * tile + tx] = (bRow < N && col < N) ? B[bRow * N + col] : 0.0f;
__syncthreads();
for (int k = 0; k < tile; k++)
sum += As[ty * tile + k] * Bs[k * tile + tx];
__syncthreads();
}
if (row < N && col < N)
C[row * N + col] = sum;
}
והlaunch, עם הפרמטר השלישי במספר בתים:
int tile = 16;
dim3 block(tile, tile);
dim3 grid((N + tile - 1) / tile, (N + tile - 1) / tile);
size_t smemBytes = (size_t)2 * tile * tile * sizeof(float);
mmTiledDyn<<<grid, block, smemBytes>>>(d_A, d_B, d_C, N, tile);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
הפלט עבור N = 4096, tile = 16:
הביצועים כמעט זהים לגרסה הסטטית (2994 מול 2968 GFLOP/s) - בתוך רעש המדידה.
למה זה עבד: ה-shared memory הדינמי מוקצה כבלוק אחד רציף בגודל שביקשנו בפרמטר השלישי (2*tile*tile*sizeof(float) בתים), ואנחנו חילקנו אותו ידנית לשני חצאים. האינדוקס המשוטח As[ty*tile + tx] מחשב בדיוק את אותה כתובת שהקומפיילר חישב עבורנו ב-As[ty][tx] הסטטי, ולכן אין הבדל ביצועים.
איך להכליל: shared memory דינמי חיוני כשגודל הtile נקבע בזמן ריצה (למשל kernel ספרייתי שמקבל tile כפרמטר), או כשצריך לחרוג ממגבלת ה-48 KiB הסטטית - אז חובה דינמי יחד עם cudaFuncSetAttribute(..., cudaFuncAttributeMaxDynamicSharedMemorySize, bytes). המחיר היחיד הוא ויתור על אינדוקס דו-ממדי נוח והצורך לחלק את הבלוק בעצמכם.
פתרון תרגיל 5 - הזרקת הבאג הקלאסי ותיקונו¶
עם המחסום השני מוסר (mmTiledBug בלי ה-__syncthreads() שאחרי הלולאה הפנימית), הפלט של הרצות חוזרות עבור N = 1024:
run 1: max abs error vs CPU = 3.7e+01
run 2: max abs error vs CPU = 0.0e+00
run 3: max abs error vs CPU = 1.1e+02
run 4: max abs error vs CPU = 5.2e+01
השגיאה משתנה בין הרצות - לפעמים אפילו אפס במקרה. זהו החתם של תנאי מרוץ (race condition). כשמוסרים את המחסום הראשון במקום זאת, מקבלים תסמין דומה (תוצאות שגויות ולא יציבות), הפעם מסכנת קריאה-לפני-כתיבה. עם שני המחסומים במקומם, השגיאה חוזרת להיות ~1e-4 יציבה בכל הרצה.
למה זה עבד (כלומר, למה הבאג מתנהג כך): בלי המחסום השני, warp שסיים את הלולאה הפנימית ממשיך לשלב t+1 ומתחיל לכתוב את הtile החדש אל אותם As/Bs, בעוד warp אחר באותו block עדיין קורא את הtile הישן. זו סכנת כתיבה-לפני-קריאה. איזה warp "מנצח" את המרוץ תלוי בתזמון הרגעי של המתזמן, ולכן התוצאה לא דטרמיניסטית - וזו בדיוק הסיבה שבאגים כאלה כה קשים לאיתור: הם נעלמים ב-printf, ב-debugger, וב-N קטן שנכנס בtile יחיד.
איך להכליל: כל אימת ש-threads של block קוראים נתון ש-threads אחרים כתבו ל-shared memory (או להפך), נדרש __syncthreads() ביניהם. הכלל: מחסום אחרי כל שלב כתיבה-משותפת ולפני השלב הקורא, ומחסום נוסף לפני שימוש חוזר בזיכרון. תוצאה לא יציבה בין הרצות היא כמעט תמיד תנאי מרוץ; חפשו מחסום חסר או גישה מרובת-threads לא מוגנת. כלים כמו compute-sanitizer --tool racecheck ./matmul מאתרים בדיוק את הסכנות האלה אוטומטית.
פתרון תרגיל 6 (בונוס) - פרופיילינג ואחוז מהשיא¶
השוואת מדדי ncu בין הנאיבי למרוצף (מקוצר):
mmNaive mmTiled(16)
Compute (SM) Throughput 6.1 % 38.5 %
DRAM Throughput 82.4 % 21.0 %
Achieved Occupancy 78.0 % 62.3 %
Shared Bank Conflicts 0 0 (compute loop)
הנאיבי צבוע בבירור memory-bound: DRAM Throughput קרוב לרוויה (82%) בעוד יחידות החישוב כמעט בטלות (6%). המרוצף מזיז את המשקל - ה-DRAM יורד ל-21% וה-SM עולה ל-38%.
חישוב אחוז השיא: הגרסה הטובה ביותר (TILE=32, ~4480 GFLOP/s) נותנת 4480 / 66900 ≈ 6.7% משיא ה-FP32 של H100.
למה זה עבד: ה-ncu מודד ישירות את שני צווארי הבקבוק המתחרים - תפוקת ה-DRAM ותפוקת ה-SM. הנתונים מאשרים כמותית את מה שחזינו: ה-tiling הוריד את הלחץ על ה-DRAM (פחות קריאות global לכל FLOP) והעלה את ניצול יחידות החישוב.
איך להכליל: ncu הוא הבורר האמפירי לשאלה "האם אני memory-bound או compute-bound". תמיד מפרופלים לפני שמכווננים - המדד שקרוב לרוויה מצביע על צוואר הבקבוק, והוא זה שיש לתקוף. הפער בין 6.7% ל-cuBLAS נסגר בשכבות הבאות: ריצוף ברמת ה-register (כל thread מחשב מיקרו-tile של C), coalescing מיטבי, ולבסוף מעבר ל-Tensor Cores. את כולן נרכיב בפרויקט 3.5, ואת התנגשויות הבנקים והריפוד [TILE_WIDTH + 1] נמצה בפרק 8.4.