2.3 כתיבת kernel אינדוקס ו launch configuration פתרון
פתרון - כתיבת kernel - אינדוקס ו-launch configuration¶
כאן נעבור על הפתרון המלא של כל תרגיל, עם הקוד והפלט המדויקים. שימו לב: הזמנים וה-GFLOP/s הם דוגמאות מהרצה על H100 SXM, והם ישתנו בין כרטיסים ובין הרצות (בעיקר במדידה). המבנה, ההיגיון והמסקנות זהים בכל מכונה. אם אתם על T4 ב-Colab, החליפו את -arch=sm_90a ב--arch=sm_75, ותקבלו GFLOP/s נמוכים יותר אך אותה תמונה איכותית של kernel memory-bound.
פתרון תרגיל 1 - מימוש ה-kernel הנאיבי¶
גוף ה-kernel המושלם:
__global__ void mm(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;
}
}
ההרצה עם N = 4 (הדפסת כל איברי C):
כל איבר הוא 2N = 8, כי כל איבר של C הוא סכום של N=4 מכפלות 1*2.
למה זה עבד: כל thread מזוהה עם איבר יחיד של C לפי (row, col), כאשר col בא מציר ה-x (הממד הרציף בזיכרון row-major) ו-row מציר ה-y. שומר הגבולות if (row < N && col < N) מגן על ה-threads העודפים שנוצרים כי ה-grid מעוגל כלפי מעלה. השיטוח row*N + col מתרגם את הקואורדינטה הדו-ממדית לאינדקס במערך החד-ממדי שבזיכרון.
איך להכליל: כל kernel שעובד על מטריצה או תמונה משתמש באותה תבנית 2D: row,col מהצירים, שומר גבולות, שיטוח. הבחירה ש-col הולך על ציר x אינה שרירותית - היא זו שתבטיח בהמשך שגישות של threads סמוכים ב-warp ייפלו על כתובות סמוכות (coalescing).
פתרון תרגיל 2 - אימות מול CPU וחישוב GFLOP/s¶
בלוק ה-TODO השני:
matmulCPU(h_A, h_B, h_ref, N);
double maxErr = 0.0;
for (int i = 0; i < N * N; i++)
maxErr = fmax(maxErr, fabs((double)h_C[i] - (double)h_ref[i]));
double flops = 2.0 * (double)N * N * N;
double gflops = flops / (ms * 1e-3) / 1e9;
printf("N = %d, kernel = %.3f ms, max error = %.6f, %.1f GFLOP/s\n",
N, ms, maxErr, gflops);
הפלט עבור N = 1024:
החשבון: 2 * 1024^3 = 2.147e9 FLOPs, חלקי 7.0e-3 שניות, חלקי 1e9, נותן כ-307 GFLOP/s. מול שיא ה-FP32 של ה-H100 (כ-66,900 GFLOP/s), זה כ-0.46%. השגיאה היא 0 כי A=1,B=2 ו-1024 * 2 = 2048 מיוצג בדיוק ב-float.
למה זה עבד: הרפרנס ב-CPU מבצע בדיוק אותו חישוב בסדר סדרתי, ולכן הוא אמת המידה לנכונות. חישוב ה-FLOPs נגזר מכך שכל אחד מ-N^2 איברי הפלט דורש N כפלים ו-N חיבורים, ומכאן 2*N^3. חילוק במספר השניות נותן קצב. החימום שהרצנו לפני התזמון מונע מדידה של קומפילציית ה-PTX ל-SASS הראשונית.
איך להכליל: 2*N^3 FLOPs הוא המדד הסטנדרטי להשוואת כל מימוש matmul (נאיבי, מרוצף, cuBLAS) - כולם עושים אותה כמות עבודה, אז GFLOP/s מודד ישירות יעילות. אחוז נמוך כל כך מהשיא הוא הרמז הראשון שמשהו יסודי חסר; תרגיל 3 יסביר בדיוק מה.
פתרון תרגיל 3 - arithmetic intensity ונקודת המפנה¶
שלב 1 - arithmetic intensity. בלולאה הפנימית: 2 קריאות global (A[...] ו-B[...]), כל אחת float = 4 בתים, סך 8 בתים. מולן 2 FLOPs (כפל אחד, חיבור אחד).
שלב 2 - נקודת המפנה של ה-H100:
שלב 3 - השוואה: 20 / 0.25 = 80. הarithmetic intensity של ה-kernel נמוכה פי 80 מנקודת המפנה. מכיוון שהוא הרבה משמאל לנקודת המפנה, הוא memory-bound באופן מובהק.
שלב 4 - התקרה memory-bound:
ceiling = 0.25 FLOP/byte * 3.35e12 bytes/s = 0.84e12 FLOP/s = 0.84 TFLOPS
= about 1.25% of the FP32 peak
שלב 5 - מול המדידה: מדדנו כ-0.307 TFLOPS (307 GFLOP/s), שאכן מתחת לתקרה של 0.84 TFLOPS. מסודר:
naive intensity = 2 FLOP / 8 bytes = 0.25 FLOP/byte
ridge point = 66.9e12 / 3.35e12 ≈ 20 FLOP/byte
distance from ridge = 20 / 0.25 = 80x (memory-bound)
ceiling = 0.25 * 3.35e12 ≈ 0.84 TFLOPS (~1.25% of peak)
measured = 2*1024^3 / 7.0e-3 ≈ 0.31 TFLOPS (below the ceiling)
למה זה עבד: מודל ה-Roofline אומר שבאזור הmemory-bound הביצועים חסומים ב-arithmetic intensity * רוחב פס. עם arithmetic intensity 0.25 בלבד, אפילו ניצול מושלם של ה-HBM3 נותן 0.84 TFLOPS - שבריר מהשיא. המדידה שלנו נמוכה עוד יותר כי הגישות ל-B הן לאורך עמודה (B[k*N+col]), פחות מאוחדות (coalesced), ולא כל הקריאות נתפסות בcache.
איך להכליל: זו השיטה לנתח כל kernel לפני שמייעלים אותו: ספרו FLOPs וספרו בתים, חלקו לקבלת arithmetic intensity, והשוו לנקודת המפנה של הכרטיס (שיא חישוב / רוחב פס). אם אתם משמאל לנקודה - אתם memory-bound, וזה בזבוז לייעל חישוב; צריך להעלות arithmetic intensity (tiling, פרק 3.3). אם מימין - compute-bound, וכדאי לפנות ל-Tensor Cores או להוריד עומס אריתמטי.
פתרון תרגיל 4 - המרה ללולאת grid-stride¶
__global__ void mmStride(const float* A, const float* B, float* C, int N) {
int col0 = blockIdx.x * blockDim.x + threadIdx.x;
int row0 = blockIdx.y * blockDim.y + threadIdx.y;
int strideX = blockDim.x * gridDim.x;
int strideY = blockDim.y * gridDim.y;
for (int row = row0; row < N; row += strideY)
for (int col = col0; col < N; col += strideX) {
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 עם grid קטן במכוון:
dim3 block(16, 16);
dim3 grid(4, 4); // 16 blocks, only 4096 threads
mmStride<<<grid, block>>>(d_A, d_B, d_C, N);
הפלט עבור שלושת הגדלים (grid של 4x4, block יחיד וכו'), כולם נכונים:
N = 1024 (grid 4x4) : max error = 0.000000 OK
N = 1000 (grid 4x4) : max error = 0.000000 OK
N = 1500 (grid 4x4) : max error = 0.000000 OK
N = 1024 (grid 1x1) : max error = 0.000000 OK (very slow)
4096 threads מטפלים במיליון איברי פלט: strideX = 16*4 = 64, strideY = 16*4 = 64, וכל thread מטפל בערך (1024/64) * (1024/64) = 16 * 16 = 256 איברים.
למה זה עבד: במקום איבר בודד לכל thread, כל thread מדלג על מטריצת התוצאה בצעדים בגודל ה-grid כולו בשני הצירים ומכסה את כל האיברים שנופלים על סריגו. התנאים row < N ו-col < N בלולאות משמשים גם כשומרי גבולות, ולכן ה-kernel נכון לכל N, כולל כאלה שאינם כפולה של גודל ה-block. הרצה עם block יחיד עדיין נכונה כי אותו thread פשוט מכסה יותר איברים - רק לאט.
איך להכליל: grid-stride היא התבנית המקצועית המומלצת. אותו kernel נכון לכל N ולכל תצורת launch, ואפשר לבחור gridDim שממקסם occupancy (כפולה של 132 ה-SMs של ה-H100) במקום להיגרר אחרי גודל הבעיה. היכולת להריץ block יחיד ולקבל תוצאה נכונה היא כלי בדיקה מצוין: אם ה-kernel נכון עם block יחיד אך שגוי עם grid מלא, הבאג הוא בתלות בין-בלוקית אסורה.
פתרון תרגיל 5 - השפעת גודל ה-block¶
טבלת מדידות ל-N = 2048, kernel נאיבי mm (H100, לדוגמה):
התשובות:
- הכי מהיר:
16x16ו-32x8(סביב 300-305 GFLOP/s) - שניהם 256 threads, גודל block מאוזן. - הכי איטי:
8x8(64 threads). הבלוקים קטנים מדי; מגבלת ה-32 בלוקים תושבים ל-SM חוסמת את התפוסה ל-32*64 = 2048threads רק אם 32 בלוקים נכנסים, ובפועל הoverhead גבוה. בנוסף,8בציר x אינו כפולה של 32, אז ה-warps נחתכים על פני שורות ופוגעים ב-coalescing. 32x32(block מקסימלי) מעט איטי יותר מ-256 כי רק 2 בלוקים נכנסים ל-SM, מה שמקשה על הסתרת latency.
למה זה עבד: ציר x בגודל כפולה של 32 (כמו ב-32x8) מיישר את ה-warps לאורך שורות רציפות בזיכרון, כך ש-32 threads warp קוראים 32 כתובות סמוכות של C ו-B - גישה מאוחדת. block בטווח 128-256 מאזן בין מספר בלוקים מספיק למלא את ה-SM לבין threads מספיקים בכל אחד.
איך להכליל: התחילו תמיד מ-16x16 או 256 threads כברירת מחדל, ודאו שציר x כפולה של 32, וכוונו משם. אבל שימו לב: כל ההבדלים כאן קטנים (פי 1.5 לכל היותר) כי כולנו memory-bound. שום בחירת block לא תציל kernel נאיבי - זה מחזק את מסקנת תרגיל 3, שהתיקון האמיתי הוא tiling, לא כיוונון תצורה.
פתרון תרגיל 6 (בונוס) - כמה רחוקים מהתקרה האמיתית¶
הרצת הפרופיילר על ה-kernel הנאיבי:
קטע אופייני מהפלט (מקוצר):
mm(const float *, const float *, float *, int)
Section: GPU Speed Of Light Throughput
----------------------------------------------------------------------
Compute (SM) Throughput [%] 12.4
Memory Throughput [%] 78.9
DRAM Throughput [%] 61.2
----------------------------------------------------------------------
Memory Throughput (כ-79%) גבוה בהרבה מ-Compute Throughput (כ-12%) - חתימה מובהקת של kernel memory-bound, בדיוק כפי שחזינו בתרגיל 3.
השוואה מול cuBLAS על אותה בעיה N=2048:
למה זה עבד: הפרופיילר מודד ישירות את שני צירי ה-Roofline. Memory Throughput גבוה ו-Compute נמוך פירושו שיחידות ה-SM מחכות לזיכרון - האבחנה הכמותית של memory-bound. cuBLAS קופץ פי 100 ומעלה כי הוא מעלה את הarithmetic intensity עם tiling רב-שכבתי (shared memory + register tiling) ומריץ את הליבה על Tensor Cores.
איך להכליל: ncu עם --set basic הוא הצעד הראשון בכל אופטימיזציה: הוא אומר לכם מיד אם אתם memory-bound או compute-bound, בלי לנחש. את ה-kernel המרוצף שיסגור חלק גדול מהפער נכתוב בעצמנו בפרק 3.3, ואת הדרך עד Tensor Cores ו-cuBLAS נשלים בפרויקט 3.5. זכרו ש-cuBLAS הוא column-major - כדי לקבל תוצאה נכונה בפריסת row-major מחשבים למעשה את המכפלה בסדר הפוך.