2.2 היררכיית ה threads blocks ו grids פתרון
פתרון - היררכיית ה-threads - blocks ו-grids¶
כאן נעבור על הפתרון המלא של כל תרגיל, עם הקוד והפלט המדויקים. שימו לב: סדר הבלוקים בתרגיל 3 והזמנים בתרגילים 4 ו-6 הם דוגמאות מהרצה על H100 - הם ישתנו בין כרטיסים ובין הרצות. המבנה, ההודעות וההיגיון זהים בכל מכונה. אם אתם על T4 ב-Colab, החליפו את -arch=sm_90a ב--arch=sm_75. בכל הפתרונות מונח שמקרו CUDA_CHECK מוגדר בראש הקובץ, כמו בשיעורים 0.3 ו-1.2.
פתרון תרגיל 1 - היפוך תמונה עם אינדוקס דו-ממדי¶
ה-kernel וlaunch שלו:
__global__ void invert(unsigned char* img, int width, int height) {
int col = blockIdx.x * blockDim.x + threadIdx.x; // .x -> column
int row = blockIdx.y * blockDim.y + threadIdx.y; // .y -> row
if (row < height && col < width) {
int idx = row * width + col;
img[idx] = 255 - img[idx];
}
}
int main(void) {
int width = 1920, height = 1080;
size_t n = (size_t)width * height;
unsigned char* h_img = (unsigned char*)malloc(n);
for (size_t i = 0; i < n; i++) h_img[i] = (unsigned char)(i % 256);
unsigned char* d_img;
CUDA_CHECK(cudaMalloc(&d_img, n));
CUDA_CHECK(cudaMemcpy(d_img, h_img, n, cudaMemcpyHostToDevice));
dim3 block(16, 16);
dim3 grid((width + block.x - 1) / block.x,
(height + block.y - 1) / block.y);
invert<<<grid, block>>>(d_img, width, height);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaMemcpy(h_img, d_img, n, cudaMemcpyDeviceToHost));
// verification: pixel (row,col) should be 255 - ((row*width+col) % 256)
int errors = 0;
for (size_t i = 0; i < n; i++)
if (h_img[i] != (unsigned char)(255 - (i % 256))) errors++;
printf("grid = %d x %d blocks, block = %d x %d\n",
grid.x, grid.y, block.x, block.y);
printf("errors = %d\n", errors);
CUDA_CHECK(cudaFree(d_img));
free(h_img);
return 0;
}
הפלט:
למה זה עבד: ה-grid הוא 120 x 68 בלוקים כי 1920/16 = 120 בדיוק, אבל 1080/16 = 67.5, וכלפי מעלה (1080+15)/16 = 68. השורה ה-68 של הבלוקים (threads בשורות 1072..1087) גולשת מעבר לגובה 1080, ושומר הגבולות row < height מוודא שאותם threads עודפים לא נוגעים בזיכרון. ה-.x מופה לעמודה וה-.y לשורה, והגישה row*width+col משטחת את הקואורדינטות הדו-ממדיות למערך החד-ממדי בפריסת row-major.
איך להכליל: זו התבנית הקנונית לכל עיבוד דו-ממדי - תמונות, מטריצות, שדות. הבחירה block(16,16) = 256 threads היא ברירת מחדל טובה (כפולה של 32, מתחלקת יפה). כשעוברים ל-3D מוסיפים .z ומחשבים depth באותה נוסחה. חשוב לשמור על עקביות: אם בחרתם .x = עמודה, שמרו על זה בכל ה-kernel כדי שthreads סמוכים ייגשו לכתובות סמוכות (coalescing) - נושא שנפתח בפרק הביצועים.
פתרון תרגיל 2 - כיסוי מטריצת NxN עם שומר גבולות¶
חישוב ידני עבור N = 1000, block(16,16):
blocks per dimension = (1000 + 15) / 16 = 63
threads per dimension = 63 * 16 = 1008
excess threads = 1008 - 1000 = 8 (per dimension)
total threads = 1008 * 1008 = 1,016,064
total elements = 1000 * 1000 = 1,000,000
excess = 16,064 threads
ה-kernel:
__global__ void scale(float* m, int N, float s) {
int col = blockIdx.x * blockDim.x + threadIdx.x;
int row = blockIdx.y * blockDim.y + threadIdx.y;
if (row == 0 && col == 0)
printf("gridDim = %d x %d, blockDim = %d x %d\n",
gridDim.x, gridDim.y, blockDim.x, blockDim.y);
if (row < N && col < N)
m[row * N + col] *= s;
}
הlaunch עם N = 1000, s = 2.0f:
dim3 block(16, 16);
dim3 grid((N + block.x - 1) / block.x, (N + block.y - 1) / block.y);
scale<<<grid, block>>>(d_m, N, 2.0f);
הפלט (אם אתחלנו כל איבר ל-1.0f, כל איבר צריך להיות 2.0f):
למה זה עבד: העיגול כלפי מעלה מבטיח grid של 63 x 63 בלוקים שמכסה את כל 1000 x 1000 האיברים, במחיר של 8 threads עודפים בכל ממד. שומר הגבולות הכפול if (row < N && col < N) הוא מה שמנטרל את 16,064 הthreads העודפים - בלעדיו הם היו ניגשים ל-m[row*N+col] עם row או col שגדולים מ-999 וכותבים מחוץ למערך, מה שהיה גורם ל-illegal memory access שנתפס ב-cudaDeviceSynchronize.
איך להכליל: בכל פעם שממדי ה-grid מעוגלים כלפי מעלה (וזה כמעט תמיד), נוצרים threads עודפים, ושומר הגבולות אינו קישוט אלא הכרח. בבעיה דו-ממדית שני התנאים חייבים להתקיים - thread יכול להיות תקף בממד אחד ועודף בשני. הכלל: מספר הthreads שהושקו הוא תמיד gridDim * blockDim >= N, אף פעם לא בדיוק N, אלא אם N מתחלק בדיוק.
פתרון תרגיל 3 - סדר הרצת הבלוקים¶
ה-kernel וlaunch שלו:
__global__ void order(int* ticket) {
if (threadIdx.x == 0) {
int t = atomicAdd(ticket, 1);
printf("block %d -> order %d\n", blockIdx.x, t);
}
}
int main(void) {
int* d_ticket;
CUDA_CHECK(cudaMalloc(&d_ticket, sizeof(int)));
CUDA_CHECK(cudaMemset(d_ticket, 0, sizeof(int)));
order<<<16, 1>>>(d_ticket);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaFree(d_ticket));
return 0;
}
פלט לדוגמה משתי הרצות עוקבות (אצלכם יהיה שונה):
# first run
block 0 -> order 0
block 4 -> order 1
block 1 -> order 2
block 8 -> order 3
block 5 -> order 4
block 2 -> order 5
...
# second run (different order!)
block 0 -> order 0
block 1 -> order 1
block 8 -> order 2
block 4 -> order 3
...
התשובות: סדר ההרצה אינו 0,1,...,15, והוא משתנה בין הרצות. זה מוכיח בפועל את חוק העצמאות: הנהג מתזמן את הבלוקים בכל סדר שנוח לו, ואסור לכתוב קוד שמניח סדר כלשהו.
למה זה עבד: ה-atomicAdd(ticket, 1) מחזיר לכל בלוק מספר תור לפי הרגע שבו הוא באמת הגיע ל-SM והתחיל לרוץ. מכיוון שהנהג מפזר בלוקים על ה-SMs לפי זמינות משאבים ולא לפי מספר הבלוק, המספרים המוחזרים אינם עוקבים. השתמשנו ב-atomic ולא בעדכון רגיל כי כמה בלוקים עשויים להגיע כמעט בו-זמנית; int t = (*ticket)++ לא-אטומי היה יוצר data race ומאבד ספירות.
איך להכליל: זו הדגמה חיה של מדוע ההרחבה השקופה עובדת - אם הקוד נכון תחת כל סדר בלוקים, הוא ירוץ נכון על SM אחד או על 132. אם אי פעם תמצאו את עצמכם מניחים ש"בלוק i רץ לפני בלוק j", עצרו - זה באג. תלות סדר אמיתית בין שלבים פותרים בשני kernel launches נפרדים, כי סיום kernel הוא נקודת הסנכרון החוצה-בלוקים היחידה שאפשר לסמוך עליה.
פתרון תרגיל 4 - סכום כלל-grid עם atomicAdd¶
ה-kernel וlaunch שלו:
__global__ void sumAtomic(const float* x, int n, float* result) {
int i = blockIdx.x * blockDim.x + threadIdx.x;
if (i < n) atomicAdd(result, x[i]);
}
int main(void) {
int n = 1 << 20;
size_t bytes = (size_t)n * sizeof(float);
float* h_x = (float*)malloc(bytes);
for (int i = 0; i < n; i++) h_x[i] = 1.0f;
float *d_x, *d_result;
CUDA_CHECK(cudaMalloc(&d_x, bytes));
CUDA_CHECK(cudaMalloc(&d_result, sizeof(float)));
CUDA_CHECK(cudaMemcpy(d_x, h_x, bytes, cudaMemcpyHostToDevice));
CUDA_CHECK(cudaMemset(d_result, 0, sizeof(float)));
int blockSize = 256;
int numBlocks = (n + blockSize - 1) / blockSize;
sumAtomic<<<numBlocks, blockSize>>>(d_x, n, d_result);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
float h_result = 0.0f;
CUDA_CHECK(cudaMemcpy(&h_result, d_result, sizeof(float),
cudaMemcpyDeviceToHost));
printf("sum = %.1f, expected = %d\n", h_result, n);
CUDA_CHECK(cudaFree(d_x)); CUDA_CHECK(cudaFree(d_result));
free(h_x);
return 0;
}
הפלט:
שאלת ההבנה - למה __syncthreads לא יכול: __syncthreads() מסנכרן רק את threads הבלוק שהם, לפי ערובת ה-SM המשותף, כולם על אותו SM. הוא אין לו שום ידיעה או השפעה על בלוקים אחרים, שיושבים על SMs אחרים ואולי אפילו עוד לא התחילו לרוץ. לכן אפשר לצרף עם __syncthreads את התוצאות של threads בתוך בלוק, אבל אי אפשר לצרף את הסכומים החלקיים בין בלוקים - אין מחסום ברמת ה-grid. הדרך היחידה לצרף בין בלוקים היא global memory עם atomics.
למה זה עבד: atomicAdd מבצע קריאה-שינוי-כתיבה על result כפעולה אחת בלתי-ניתנת-לחלוקה, כך שגם כשמיליון threads מכל 4096 הבלוקים מוסיפים בו-זמנית לאותה כתובת, אף עדכון לא נדרס. הסכום המדויק 1048576.0 מתקבל כי כל האיברים הם 1.0f והסכום עדיין מיוצג בדיוק ב-float.
איך להכליל: atomics הם הכלי הכללי לתקשורת חוצה-בלוקים (מונים, היסטוגרמות, רדוקציות, רשימות עבודה). החיסרון: כשכל הthreads מתנגשים על כתובת אחת, ה-atomics מסתרלים (מבוצעים כמעט בזה אחר זה) וההאצה נפגעת. לכן הגרסה הזו נכונה אך לא אופטימלית - ראו את תרגיל 6 לתבנית היעילה שמצמצמת את מספר ה-atomics פי blockSize.
פתרון תרגיל 5 - חישוב occupancy ביד ואימות¶
הקומפילציה עם דיווח משאבים:
פלט ה-ptxas עבור sumAtomic (דוגמה - הערכים ישתנו מעט לפי גרסת ה-toolkit):
ptxas info : Compiling entry function '_Z9sumAtomicPKfiPf' for 'sm_90a'
ptxas info : Used 16 registers, 372 bytes cmem[0]
חישוב ידני עבור blockSize = 256 (8 warps), 16 אוגרים/thread, בלי shared memory, על H100:
register limit : 65,536 / (16 * 256) = 65,536 / 4,096 = 16 blocks
warps limit : 64 / 8 = 8 blocks
shared limit : no shared -> unlimited
blocks limit : hard maximum = 32 blocks
------------------------------------------------
limiting factor = minimum = 8 blocks
האימות עם ה-Runtime:
int maxBlocks = 0;
CUDA_CHECK(cudaOccupancyMaxActiveBlocksPerMultiprocessor(
&maxBlocks, sumAtomic, 256, 0));
printf("max active blocks / SM = %d\n", maxBlocks);
printf("occupancy = %.0f%%\n", (maxBlocks * 256) / 2048.0 * 100.0);
הפלט:
למה זה עבד: ל-kernel קטן (16 אוגרים, בלי shared) הגורם המגביל אינו האוגרים (שהיו מאפשרים 16 בלוקים) אלא חריצי ה-warp: 64 warps ל-SM חלקי 8 warps לבלוק נותנים 8 בלוקים. שמונה בלוקים של 256 threads הם 2,048 threads - בדיוק המקסימום התושב, כלומר תפוסה מלאה. ה-API של cudaOccupancyMaxActiveBlocksPerMultiprocessor מבצע בדיוק את החישוב הזה בשבילנו ומאשר את המספר.
איך להכליל: תמיד חשבו את המשאב המגביל לפי רצפה(תקציב / צריכה-לבלוק) על פני כל המשאבים, וקחו את המינימום. כדי לראות את לחץ האוגרים בפעולה, הוסיפו __launch_bounds__(256, 4) או כתבו kernel שצורך הרבה אוגרים - תראו את מספר הבלוקים והתפוסה יורדים. ה-API הזה הוא הבסיס גם ל-cudaOccupancyMaxPotentialBlockSize, שבוחר עבורכם את גודל הבלוק שממקסם תפוסה. את התמונה המלאה - ומתי תפוסה גבוהה בכלל עוזרת - נפתח בפרק 8.
פתרון תרגיל 6 (בונוס) - רדוקציה בשני שלבים: shared memory ואז atomic¶
ה-kernel המלא:
#define BLOCK 256
__global__ void sumReduce(const float* x, int n, float* result) {
__shared__ float sdata[BLOCK];
int tid = threadIdx.x;
int i = blockIdx.x * blockDim.x + threadIdx.x;
sdata[tid] = (i < n) ? x[i] : 0.0f; // load into shared (0 for excess threads)
__syncthreads(); // barrier at the block level
for (int s = blockDim.x / 2; s > 0; s >>= 1) { // tree reduction
if (tid < s) sdata[tid] += sdata[tid + s];
__syncthreads(); // barrier between tree stages
}
if (tid == 0) atomicAdd(result, sdata[0]); // cross-block synchronization
}
הlaunch זהה לתרגיל 4 (numBlocks = (n + BLOCK - 1) / BLOCK, blockSize = BLOCK). הפלט וההשוואת זמנים לדוגמה (H100):
sumAtomic (atomic per element) : sum = 1048576.0, time = 0.62 ms
sumReduce (shared + atomic) : sum = 1048576.0, time = 0.021 ms
שלב 3 - היכן כל סנכרון: ה-__syncthreads() (שני מופעים - אחרי הטעינה ובתוך לולאת העץ) הוא סנכרון ברמת הבלוק: הוא מבטיח שכל threads הבלוק סיימו לכתוב ל-sdata לפני שמישהו קורא ממנה, אפשרי כי כולם על אותו SM וחולקים את אותה shared memory. ה-atomicAdd בסוף הוא סנכרון חוצה-בלוקים: הוא מצרף את 4096 הסכומים החלקיים (אחד לבלוק) לתוצאה גלובלית אחת דרך global memory, כי אין דרך אחרת לתאם בין בלוקים.
למה זה עבד: כל בלוק מצמצם את 256 האיברים שלו לערך יחיד sdata[0] בעזרת רדוקציה בעץ ב-shared memory המהיר, ואז מבצע atomicAdd אחד בלבד. במקום 1,048,576 פעולות atomic מתנגשות (גרסת תרגיל 4), יש רק 4,096 - פי 256 פחות התנגשות על ה-global - ומכאן ההאצה של פי כ-30. ה-__syncthreads() בתוך הלולאה הכרחי לנכונות: בלעדיו thread עלול לקרוא sdata[tid+s] לפני שהthread השני סיים לכתוב אליו, וזה data race שמחזיר תוצאה שגויה ולא דטרמיניסטית.
איך להכליל: התבנית "צמצם בתוך הבלוק ב-shared, ואז צרף בין בלוקים ב-atomic" היא הדרך הסטנדרטית לכל רדוקציה כלל-grid (סכום, מקסימום, נורמה, היסטוגרמה). היא ממחישה בול את שתי רמות שיתוף הפעולה של ההיררכיה: מהיר וברמת הבלוק בפנים (shared + barriers), איטי ומבוסס-atomics בין הבלוקים (global). את הגרסאות המשופרות עוד יותר - warp shuffles שמוותרים על חלק מה-__syncthreads, ורדוקציה ללא bank conflicts - נפתח בפרק על shared memory וב-פרק הביצועים.