2.3 כתיבת kernel אינדוקס ו launch configuration הרצאה

בשיעורים 2.1 ו-2.2 בנינו את המודל התיאורטי של CUDA - היררכיית ה-threads (thread, block, grid), הכלל שבלוקים חייבים להיות עצמאיים, והמיפוי של block שלם ל-SM (מעבד רב-זרמי) יחיד. בשיעור 0.3 אפילו כתבנו kernel רץ ראשון, חיבור וקטורים חד-ממדי. עכשיו נסגור את הפער בין התיאוריה לקוד: נלמד לכתוב ולהשיק kernel נכון - את תצורת ההרצה המלאה בת ארבעת הפרמטרים, את תבניות האינדוקס לבעיות חד-ממדיות, דו-ממדיות ותלת-ממדיות, ואת לולאת ה-grid-stride שהופכת kernel לנייד לכל גודל בעיה. ואז נכתוב את ה-kernel שילווה אותנו לאורך כל הקורס: כפל המטריצות (matmul) הנאיבי. הוא יעבוד, ייתן תוצאה נכונה, ובכל זאת ינצל שבריר זעום מכוח ה-GPU. הוא ילמד אותנו את השיעור הכי חשוב בביצועי GPU - קיר רוחב הפס - שבו נבין למה עצם הרצת חישוב על ה-GPU אינה מבטיחה שום דבר, ומאיזה כיוון מגיע התיקון. זה השיעור שמכין את הקרקע לכל פרק 3.

אנטומיה של ה-kernel - פונקציית `global`¶

נחזור על ההגדרה המדויקת של kernel, כי כל טעות עליה מתגלגלת הלאה. ה-kernel הוא פונקציית C++ רגילה, המסומנת במילת המפתח __global__, עם שלושה מאפיינים מחייבים:

טיפוס ההחזרה הוא תמיד void. kernel לא מחזיר ערך. אם כתבתם __global__ int foo(...), ה-nvcc יסרב לקמפל. התוצר היחיד של kernel הוא תופעת הלוואי שלו - כתיבות לזיכרון.
התקשורת היא דרך מצביעים ל-global memory (זיכרון גלובלי). הhost (ה-CPU) מעביר ל-kernel כתובות של זיכרון שיושב על ה-GPU, וה-kernel קורא וכותב דרכן. אין ערך מוחזר; יש רק מצביעים ש"מְשַנים זיכרון".
מושק פעם אחת, מבוצע פעמים רבות. זו האמירה המגדירה של CUDA. קריאה אחת ל-kernel (launch) יוצרת grid שלם של threads, וכל אחד מהם מריץ את אותו גוף קוד בדיוק. בעולם ה-CPU קריאה אחת = הרצה אחת; ב-CUDA launch אחת = אלפי הרצות מקבילות.

הנה השלד:

__global__ void myKernel(const float* in, float* out, int n) {
    // body executed by each of many threads; distinguished from one another by index
}

מלכוד שמות שחשוב לחזור עליו: המילה __global__ אין לה שום קשר ל-global memory. היא מסמנת "פונקציה שמושקת מהhost ורצה על הdevice". אחד ממעצבי CUDA אמר בבדיחות הדעת שהשם נבחר כדי למקסם בלבול. זכרו: __global__ = kernel; global memory = ה-RAM של הכרטיס. שני מושגים שונים לגמרי שחולקים מילה במקרה. (למען השלמות: __device__ מסמן פונקציה שרצה על הdevice ונקראת רק מתוך device code, ו-__host__ היא ברירת המחדל של host code. __global__ הוא הגשר היחיד ביניהם - host code שמפעיל device code.)

תצורת ההרצה המלאה - the full launch configuration¶

בשיעור 0.3 השקנו kernel עם שני פרמטרים בסוגריים המשולשים: <<<numBlocks, blockSize>>>. זו הצורה המקוצרת. הצורה המלאה כוללת ארבעה פרמטרים:

myKernel<<<gridDim, blockDim, sharedBytes, stream>>>(args...);

נפרק אותם אחד-אחד:

פרמטר	טיפוס	משמעות	ברירת מחדל
`gridDim`	`dim3`	מבנה ה-grid: כמה בלוקים, ובאילו ממדים	חובה
`blockDim`	`dim3`	מבנה כל block: כמה threads, ובאילו ממדים	חובה
`sharedBytes`	`size_t`	כמות ה-shared memory הדינמי לכל block, בבתים	0
`stream`	`cudaStream_t`	הזרם שאליו הlaunch נכנסת (סדר הביצוע)	0 (הזרם הראשי)

שני הפרמטרים הראשונים הם מטיפוס dim3 - מבנה בן שלושה שדות x, y, z. אפשר לאתחל אותו בממד אחד, בשניים או בשלושה; שדות שלא צוינו מקבלים 1:

dim3 a(256);            // a.x=256, a.y=1, a.z=1   (one-dimensional)
dim3 b(16, 16);         // b.x=16,  b.y=16, b.z=1  (two-dimensional, 256 threads)
dim3 c(8, 8, 8);        // 512 threads in three dimensions

כשמעבירים מספר שלם פשוט כמו ב-<<<numBlocks, blockSize>>>, הקומפיילר עוטף אותו אוטומטית ל-dim3 חד-ממדי. לכן הצורה מ-0.3 היא בסך הכל מקרה פרטי של הצורה המלאה, עם sharedBytes = 0 ו-stream = 0.

sharedBytes קובע כמה shared memory (זיכרון משותף) דינמי יוקצה לכל block. זה רלוונטי רק ל-kernel שמשתמש ב-extern __shared__; נעמיק בזה בפרק 3. עד אז נשאיר אותו 0 (או נשמיט).
stream קובע לאיזה זרם (stream) הlaunch שייכת. פעולות באותו זרם מבוצעות בסדר; פעולות בזרמים שונים יכולות לחפוף. הזרם 0 הוא ברירת המחדל. את הזרמים והחפיפה בין העברות זיכרון לחישוב נלמד בפרק על ריצה אסינכרונית.

מגבלות חומרה שכדאי לזכור: המכפלה blockDim.x * blockDim.y * blockDim.z חייבת להיות עד 1024 threads ל-block. ממדי ה-gridDim מוגבלים ל-2^31 - 1 בציר x, ו-65535 בצירים y ו-z. חריגה מכל אלה נותנת cudaErrorInvalidConfiguration כבר בlaunch.

חישוב אינדקס - 1D, 2D and 3D indexing¶

השקנו grid של אלפי threads שמריצים אותו קוד; איך כל thread יודע על איזה נתון לעבוד? התשובה, כפי שראינו ב-2.1, היא אינדוקס מבוסס מזהים. כל thread רואה ארבעה משתני מבנה מובנים - threadIdx, blockIdx, blockDim, gridDim - כל אחד עם שדות .x, .y, .z, ומהם הוא גוזר את מיקומו.

המקרה החד-ממדי (1D) - למערכים, כפי שראינו בחיבור וקטורים:

int i = blockIdx.x * blockDim.x + threadIdx.x;

המקרה הדו-ממדי (2D) - למטריצות ותמונות. שני צירים בלתי תלויים, כל אחד עם אותה נוסחה:

int col = blockIdx.x * blockDim.x + threadIdx.x;   // column index
int row = blockIdx.y * blockDim.y + threadIdx.y;   // row index
int idx = row * width + col;                        // flatten to a 1D array in memory

שימו לב לנקודה שמבלבלת מתחילים: הזיכרון תמיד חד-ממדי. מטריצה N x N שמורה כמערך רציף של N*N איברים בסדר שורות (row-major); האיבר (row, col) יושב באינדקס row * N + col. הצירים הדו-ממדיים של ה-grid הם רק דרך נוחה למספֵּר threads; הגישה לזיכרון עדיין דורשת שיטוח ידני.

   2D grid of 16x16 blocks, covering an NxN matrix
   ┌──────────┬──────────┬──────────┐
   │ block    │ block    │ block    │   blockIdx.y = 0
   │ (0,0)    │ (1,0)    │ (2,0)    │
   ├──────────┼──────────┼──────────┤
   │ block    │ block    │ block    │   blockIdx.y = 1
   │ (0,1)    │ (1,1)    │ (2,1)    │
   └──────────┴──────────┴──────────┘
     each block contains 16x16 = 256 threads;
     thread (tx,ty) within block (bx,by) handles
     col = bx*16 + tx ,  row = by*16 + ty

המקרה התלת-ממדי (3D) - לנפחים (סימולציות פיזיקליות, וידאו). אותו רעיון, ציר שלישי:

int x = blockIdx.x * blockDim.x + threadIdx.x;
int y = blockIdx.y * blockDim.y + threadIdx.y;
int z = blockIdx.z * blockDim.z + threadIdx.z;
int idx = (z * height + y) * width + x;

בכל שלושת המקרים, כמו בחיבור הווקטורים, נדרש שומר גבולות (bounds guard) כי ה-grid כמעט אף פעם לא מחלק בדיוק את מידות הבעיה: if (row < N && col < N). הthreads העודפים בשוליים חייבים לצאת בלי לגעת בזיכרון.

בחירת גודל ה-block - choosing the block size¶

כמה threads לשים ב-block? יש שני שיקולים, אחד מחייב ואחד של כיוונון:

כלל מחייב: כפולה של 32. בחומרה, הthreads מתבצעים בקבוצות של 32 הנקראות warp (קבוצת 32 threads המבצעים אותה פקודה יחד, כפי שראינו ב-2.2). block שאינו כפולה של 32 יגרור warp אחרון חלקי, שבו חלק מנתיבי החישוב פשוט מושבתים ומתבזבזים. לכן blockDim (המכפלה הכוללת) צריך תמיד להיות כפולה של 32.

כיוונון: 128 עד 256 בדרך כלל. ערכים נפוצים הם 128, 256 ו-512. 256 הוא ברירת מחדל טובה שכמעט תמיד סבירה. בבעיות דו-ממדיות בלוק של 16 x 16 = 256 הוא הבחירה הקלאסית, וגם 32 x 8 = 256 או 32 x 32 = 1024 נפוצים.

הסיבה לטווח הזה קשורה ל-occupancy (תפוסה) - כמה warps יכולים לשבת (resident) בו-זמנית על SM. ב-H100 כל SM מחזיק עד 2048 threads תושבים (64 warps) ועד 32 בלוקים תושבים. אם ה-block קטן מדי (נניח 32 threads), מגבלת ה-32 בלוקים תחסום אתכם ל-32 * 32 = 1024 threads בלבד - חצי תפוסה. אם ה-block גדול (1024), מכסת ה-threads מאפשרת רק 2 בלוקים לכל SM, מה שמקשה על החומרה להסתיר latency (השהיה) כשבלוק תקוע בהמתנה לזיכרון. block בטווח 128-256 מאזן בין השניים: מספיק בלוקים כדי למלא את ה-SM, מספיק threads בכל אחד כדי לארוז warps שלמים. את החישוב המדויק של occupancy - איך registers ו-shared memory מגבילים אותה - נעשה בפרק 8.1; כאן די בכלל האצבע.

כפל מטריצות נאיבי - the naive matmul kernel¶

עכשיו נכתוב את ה-kernel המרכזי של הקורס. הבעיה: בהינתן שתי מטריצות מרובעות A ו-B בגודל N x N, לחשב את מכפלתן C = A * B, כאשר

C[row][col] = sum over k of  A[row][k] * B[k][col]

ההקבלה הטבעית ל-CUDA: thread אחד לכל איבר של C. כל thread מזהה את השורה והעמודה שלו לפי האינדוקס הדו-ממדי, ומריץ לולאה על k שסוכמת את המכפלה. זהו ה-kernel כפי שהוא מופיע ב-glossary:

__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];   // multiply-add on two global elements
        }
        C[row * N + col] = sum;
    }
}

והlaunch, עם block דו-ממדי ו-grid דו-ממדי שמכסה את C בעיגול כלפי מעלה:

dim3 block(16, 16);                              // 256 threads
dim3 grid((N + block.x - 1) / block.x,
          (N + block.y - 1) / block.y);
mm<<<grid, block>>>(d_A, d_B, d_C, N);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());

נזכיר את מקרו הבדיקה CUDA_CHECK מ-0.3, כי נשתמש בו לאורך כל הקורס. הוא לוכד את קוד השגיאה שכל קריאת Runtime API מחזירה, מדפיס קובץ, שורה, הקריאה וההודעה, ויוצא:

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

ה-kernel נכון. אם תריצו אותו מול רפרנס CPU תקבלו התאמה מלאה (עד רעש נקודה צפה זעום). קמפלו עם nvcc -O2 -arch=sm_90a -o matmul matmul.cu והריצו - הוא יעבוד. ובכל זאת, אנחנו עומדים לגלות שהוא אסון ביצועים.

קיר רוחב הפס - the memory bandwidth wall¶

בואו נספור את התעבורה מול הזיכרון בלולאה הפנימית, איטרציה אחת:

   sum += A[row*N + k] * B[k*N + col];
          └── global read ──┘   └── global read ──┘
          2 reads (8 bytes)  ,  2 FLOPs (multiply + add)

בכל איטרציה אנחנו קוראים איבר אחד מ-A ואיבר אחד מ-B - שתי קריאות מ-global memory, 8 בתים - ומבצעים כפל וחיבור, כלומר שתי FLOPs. היחס:

   arithmetic intensity = 2 FLOPs / 8 bytes = 0.25 FLOP per byte

זו arithmetic intensity נמוכה באופן קיצוני. גרוע מזה: אין כאן שום ניצול-חוזר (reuse). כל N ה-threads שמחשבים שורה נתונה של C קוראים את אותה שורה של A - אבל כל אחד קורא אותה מ-global memory מחדש, בנפרד. אותו נתון נשלף מה-HBM3 שוב ושוב. מטריצה N x N מבצעת 2*N^3 FLOPs ומחוללת בסדר גודל של 2*N^3 קריאות global.

עכשיו נשווה לחומרה. ב-H100 SXM:

שיא ה-FP32 (ללא Tensor Cores) הוא כ-66.9 TFLOPS = 66.9e12 FLOP בשנייה.
רוחב הפס של ה-HBM3 הוא כ-3.35 TB/s = 3.35e12 בתים בשנייה.

הכלי לחבר בין השניים הוא מודל ה-Roofline. הוא מגדיר את נקודת המפנה (ridge point) - הarithmetic intensity שבה החישוב עובר מלהיות memory-bound לcompute-bound:

   ridge point = peak compute / bandwidth
              = 66.9e12 FLOP/s / 3.35e12 bytes/s
              ≈ 20 FLOP per byte

כלומר: כדי לא להיות memory-bound על H100, חישוב חייב לבצע לפחות כ-20 FLOPs על כל בית שהוא קורא. וה-matmul הנאיבי מספק 0.25. הוא רחוק פי 80 מנקודת המפנה:

        memory-bound          │      compute-bound
   ─────────────────────────┼──────────────────────►  arithmetic intensity (FLOP/byte)
   0.25          4        20 (ridge point, H100)
    ▲            ▲         ▲
   naive      tiled(16)   from here on the CUDA Cores are the bottleneck

מה זה אומר על המהירות המקסימלית האפשרית? באזור הmemory-bound, המהירות המושגת חסומה על ידי arithmetic intensity * רוחב פס:

   naive ceiling = 0.25 FLOP/byte * 3.35e12 bytes/s
              = 0.84e12 FLOP/s = 0.84 TFLOPS

גם אם ננצל את מלוא רוחב הפס של ה-HBM3, ה-kernel הנאיבי חסום ל-0.84 TFLOPS - כ-1.25% בלבד מ-66.9 TFLOPS השיא. זו בדיוק האזהרה שה-glossary מנסח בחדות: "לעולם לא תנצלו כך את כל ה-GPU, כי רוחב הפס האריתמטי של ה-CUDA Cores ב-FLOPs/s גבוה בהרבה מרוחב הפס בין ה-GPU RAM ל-SMs." יחידות החישוב מורעבות; הן מחכות לזיכרון. בפועל, בגלל L1/L2 caches שקולטים חלק מהקריאות החוזרות, נמדוד סביב 300 GFLOP/s - בערך שליש מהתקרה הmemory-bound, ועדיין פי 200 מתחת לשיא ה-FP32.

כיוון הפתרון - raising arithmetic intensity¶

אם הבעיה היא arithmetic intensity נמוכה, הפתרון חייב להיות להעלות אותה: לקרוא כל נתון מ-global memory פעם אחת, להניח אותו במחסן מהיר על ה-SM, ולסחוט ממנו כמה שיותר FLOPs לפני שזורקים אותו. המחסן הזה הוא ה-shared memory, והשיטה נקראת tiling (חלוקה-לאריחים):

מחלקים את C לtiles בגודל TILE_WIDTH x TILE_WIDTH, block אחד לכל tile.
ה-block טוען tile של A וtile של B ל-shared memory פעם אחת, ואז כל ה-threads שלו מחשבים מתוך המחסן המהיר.
כל איבר טעון מנוצל TILE_WIDTH פעמים במקום פעם אחת.

התוצאה: הarithmetic intensity מזנקת מ-1 FLOP לקריאה ל-TILE_WIDTH FLOPs לקריאה. ל-TILE_WIDTH = 16 מקבלים 16 FLOPs לקריאה, כלומר 16 FLOP / 4 בתים = 4 FLOP/בית - פי 16 מהנאיבי, וכבר קרוב יותר לנקודת המפנה (עדיין מתחתיה, ולכן עדיין קצת memory-bound, אבל הרבה פחות). ל-TILE_WIDTH = 32 מגיעים ל-8 FLOP/בית. את המימוש המלא של ה-kernel המרוצף, על שני מחסומי ה-__syncthreads() שלו, נכתוב בפרק 3 (שיעור 3.3), והוא ייתן קפיצת מהירות של פי 10 עד פי 15.

ואפילו זה אינו התקרה. גם ה-kernel המרוצף רץ על CUDA Cores. ה-matmul המהיר ביותר רץ על Tensor Cores (ליבות טנזור), שרוחב הפס האריתמטי שלהן גבוה בסדר גודל, ומגיע לעשרות אלפי GFLOP/s. את זה נראה בפרויקט 3.5 ובפרק על Tensor Cores. הרעיון המרכזי שכדאי לקחת מכאן: עצם הרצת חישוב על ה-GPU לא אומרת כלום. מה שקובע הוא היחס בין חישוב לתעבורת זיכרון, ואת היחס הזה מעלים בעבודה מודעת עם היררכיית הזיכרון.

לולאת grid-stride - portability for any N¶

ל-kernel הנאיבי שכתבנו יש הנחה סמויה: ה-grid חייב להיות גדול מספיק כדי לתת thread אחד לכל איבר של C. עבור N ענק, או כשרוצים לקבע את מספר הבלוקים למספר ה-SMs במקום ל-N, ההנחה הזו מגבילה. הפתרון הנקי, כפי שראינו במקרה החד-ממדי ב-0.3, הוא לולאת 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;   // total threads on the x axis
    int strideY = blockDim.y * gridDim.y;   // total threads on the y axis

    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;
        }
    }
}

הרעיון: כל thread לא מטפל באיבר בודד אלא מדלג על פני מטריצת התוצאה בצעדים בגודל ה-grid כולו, בשני הצירים, ומטפל בכל האיברים שנופלים על סריגו. שני התנאים row < N ו-col < N משמשים גם כשומרי גבולות, כך שאין צורך ב-if נפרד. היתרונות:

גמישות מלאה בגודל. אותו kernel נכון לכל N - קטן, ענק, לא כפולה של גודל ה-block - בלי לגעת בתצורת הlaunch.
ניתוק ה-grid מהבעיה. אפשר לבחור gridDim שממקסם occupancy (למשל כפולה של 132 ה-SMs של ה-H100) במקום להיגרר אחרי N.
בדיקת נכונות קלה. אפשר להריץ את אותו kernel עם grid זעום במכוון - אפילו block יחיד <<<1, dim3(16,16)>>> - ולוודא שהוא נכון לפני שמגדילים.

בתרגול תממשו את שתי הגרסאות, תמדדו את הנאיבית ותחשבו GFLOP/s, תאשרו במספרים שהיא memory-bound מול נקודת המפנה של ה-H100, ותהפכו אותה ל-grid-stride נייד.

סיכום¶

ה-kernel מסומן ב-__global__, מחזיר תמיד void, מתקשר רק דרך מצביעים ל-global memory, ומושק פעם אחת אך מבוצע בידי grid שלם של threads; __global__ אין לו קשר ל-global memory.
תצורת ההרצה המלאה היא <<<gridDim, blockDim, sharedBytes, stream>>>: שני הראשונים מטיפוס dim3 (ממדי ה-grid וה-block), השלישי הוא בתי shared memory דינמי (ברירת מחדל 0), הרביעי הוא הזרם (ברירת מחדל 0).
dim3 הוא מבנה בן שלושה שדות שממלא 1 בממדים שלא צוינו; מספר שלם פשוט בlaunch נעטף אוטומטית ל-dim3 חד-ממדי.
האינדקס נגזר מ-blockIdx.* * blockDim.* + threadIdx.* בכל ציר; בעיות 1D/2D/3D משתמשות באותה נוסחה לכל ציר, אך הזיכרון תמיד חד-ממדי ודורש שיטוח ידני (row * N + col) ושומר גבולות.
בוחרים blockDim ככפולה של 32 (בגלל ה-warp), בטווח טיפוסי של 128-256, כדי לאזן בין מספר הבלוקים למספר ה-threads ולמקסם occupancy; block קטן מדי או גדול מדי פוגע בתפוסה.
ה-matmul הנאיבי מקצה thread אחד לכל איבר של C ומריץ לולאת k על global memory; הוא נכון אך מבצע כ-2 FLOPs לכל 8 בתים - arithmetic intensity של 0.25 FLOP/בית, בלי שום ניצול-חוזר.
נקודת המפנה (ridge point) של ה-H100 היא 66.9 TFLOPS / 3.35 TB/s ≈ 20 FLOP/בית; הarithmetic intensity הנאיבית 0.25 רחוקה פי 80 מתחתיה, ולכן ה-kernel memory-bound וחסום לתקרה של כ-0.84 TFLOPS - כ-1.25% מהשיא.
הכיוון לתיקון הוא להעלות את הarithmetic intensity עם tiling דרך shared memory (מגיע ל-TILE_WIDTH FLOPs לקריאה, פרק 3.3), והתקרה המוחלטת היא Tensor Cores; עצם ההרצה על GPU אינה מבטיחה ביצועים.
לולאת grid-stride דו-ממדית (for row += blockDim.y*gridDim.y, for col += blockDim.x*gridDim.x) מנתקת את ה-grid מגודל הבעיה, נכונה לכל N ולכל תצורת launch, וניתנת להרצה עם block יחיד לבדיקת נכונות.