לדלג לתוכן

0.3 התוכנית הראשונה חיבור וקטורים פתרון

פתרון - התוכנית הראשונה - חיבור וקטורים

כאן נעבור על הפתרון המלא של כל תרגיל, עם הקוד והפלט המדויקים. שימו לב: הזמנים והמספרים בהמשך הם דוגמאות מהרצה על H100, והם ישתנו בין כרטיסים - בעיקר בתרגיל המדידה. המבנה, ההודעות וההיגיון זהים בכל מכונה. אם אתם על T4 ב-Colab, החליפו את -arch=sm_90a ב--arch=sm_75.

פתרון תרגיל 1 - הקלדה, קומפילציה והרצה

הקוד הוא בדיוק vecadd.cu מההרצאה. הקומפילציה וההרצה:

nvcc -O2 -arch=sm_90a -o vecadd vecadd.cu
./vecadd

הפלט עבור n = 1 << 20:

n = 1048576, max error = 0.000000
c[0] = 3.0, c[1048575] = 3.0

לאחר שינוי n ל-1 << 25:

n = 33554432, max error = 0.000000
c[0] = 3.0, c[33554431] = 3.0

למה זה עבד: כל איבר במערך הוקצה ל-thread אחד לפי i = blockIdx.x * blockDim.x + threadIdx.x, ושומר הגבולות if (i < n) הגן על הthreads העודפים. מכיוון ש-1.0f + 2.0f שווה בדיוק ל-3.0f בייצוג נקודה צפה, השגיאה המקסימלית היא אפס מדויק. הגדלת n ל-1 << 25 רק הגדילה את מספר הבלוקים ((33554432 + 255)/256 = 131072), והנוסחה עם העיגול כלפי מעלה טיפלה בזה מבלי שנגענו ב-kernel.

איך להכליל: אותו שלד - הקצה, העתק H2D, launch, העתק D2H, שחרר - הוא המבנה של כמעט כל תוכנית CUDA. שינוי גודל הבעיה לא דורש שינוי ב-kernel כל עוד תצורת ה-grid נגזרת מ-n. הגבול היחיד הוא זיכרון ה-GPU: שלושה מערכים של 4 * n בתים חייבים להיכנס ל-global memory, ו-CUDA_CHECK על cudaMalloc הוא מה שיתפוס out of memory בקצה.

פתרון תרגיל 2 - שגיאה מכוונת וטיפול בה

שלב 2 - שגיאת launch. עם blockSize = 2048 הפלט הוא:

CUDA error at vecadd.cu:52
  'cudaGetLastError()'
  -> invalid configuration argument

השגיאה נתפסה כבר ב-cudaGetLastError(), לפני שה-kernel בכלל התחיל, כי 2048 threads ל-block חורגים מהמקסימום 1024.

שלב 3 - שגיאת ריצה. אחרי הסרת if (i < n) והוספת numBlocks += 100;, הthreads העודפים ניגשים מעבר לסוף המערך:

CUDA error at vecadd.cu:53
  'cudaDeviceSynchronize()'
  -> an illegal memory access was encountered

הפעם cudaGetLastError() אחרי הlaunch עבר בשלום (התצורה חוקית), אבל הגישה מחוץ לתחום התגלתה רק כשה-kernel רץ בפועל, ולכן נתפסה ב-cudaDeviceSynchronize().

שלב 4 - ההסבר. שני סוגי שגיאות, שתי נקודות תפיסה:

סוג השגיאה מתי מתגלה נתפס על ידי
תצורת הרצה לא חוקית (blockSize > 1024) ברגע הlaunch cudaGetLastError()
גישה לא חוקית לזיכרון בתוך ה-kernel בזמן הריצה על ה-GPU cudaDeviceSynchronize()

למה זה עבד: הlaunch kernel היא אסינכרונית - השורה חוזרת לhost לפני שה-GPU סיים. שגיאה שהחומרה יכולה לזהות מיד (תצורה פסולה) זמינה מ-cudaGetLastError מיד; שגיאה שמתרחשת רק במהלך החישוב (גישה מחוץ לתחום) דורשת סנכרון כדי לחכות לה ולתפוס אותה. שתי הבדיקות משלימות זו את זו.

איך להכליל: תמיד עטפו launch kernel בזוג cudaGetLastError() ואז cudaDeviceSynchronize() בזמן פיתוח וניפוי. שגיאות CUDA "דביקות": ברגע שקרתה גישה לא חוקית, ה-context שבור וכל קריאה שאחריה תיכשל עם אותה הודעה - אז אם אתם רואים שגיאה מוזרה, חפשו את ה-kernel הראשון שנכשל, לא את הקריאה שבה הבחנתם בה. בפרודקשן אפשר להסיר את הסנכרון היקר, אבל את בדיקת cudaGetLastError משאירים.

פתרון תרגיל 3 - מימוש SAXPY

__global__ void saxpy(int n, float a, const float* x, float* y) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < n) {
        y[i] = a * x[i] + y[i];
    }
}

השינויים ב-main (על גבי שלד vecadd.cu, עם שני מערכים בלבד d_x, d_y):

float a = 3.0f;
for (int i = 0; i < n; i++) { h_x[i] = 1.0f; h_y[i] = 2.0f; }
// ... cudaMalloc d_x, d_y ; cudaMemcpy H2D of both ...
int blockSize = 256;
int numBlocks = (n + blockSize - 1) / blockSize;
saxpy<<<numBlocks, blockSize>>>(n, a, d_x, d_y);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaMemcpy(h_y, d_y, bytes, cudaMemcpyDeviceToHost));  // the result is in y

הפלט (כל איבר צריך להיות 3*1 + 2 = 5):

n = 1048576, a = 3.0, max error = 0.000000
y[0] = 5.0, y[1048575] = 5.0

ה-max error הצפוי הוא 0, כי 3.0f * 1.0f + 2.0f == 5.0f בדיוק.

למה זה עבד: SAXPY זהה במבנה ל-vecAdd - אותו חישוב אינדקס, אותו שומר גבולות - ההבדל היחיד הוא שהמערך y משמש גם כקלט וגם כפלט (קריאה-שינוי-כתיבה על אותה כתובת), וש-a הוא סקלר שעובר by-value. סקלר לא דורש cudaMalloc: רק מצביעים מפנים לזיכרון device; int ו-float פשוט נדחפים לתצורת הארגומנטים של הlaunch.

איך להכליל: רוב ה-kernels ה"אלמנטריים" (חיבור, כפל, ReLU, scale) הם אותו תבנית map: חשב i, שמור גבולות, בצע פעולה נקודתית. SAXPY במיוחד היא לבנת יסוד של BLAS Level 1 וממומשת ב-cuBLAS כ-cublasSaxpy. כשמערך משמש קלט-ופלט, שימו לב שאין תלות בין איברים שונים (כל i עצמאי) - אחרת היינו זקוקים לסנכרון או לחיץ נפרד.

פתרון תרגיל 4 - המרה ללולאת grid-stride

__global__ void vecAddStride(const float* a, const float* b,
                             float* c, int n) {
    int idx    = blockIdx.x * blockDim.x + threadIdx.x;
    int stride = blockDim.x * gridDim.x;
    for (int i = idx; i < n; i += stride) {
        c[i] = a[i] + b[i];
    }
}

הlaunch עם grid קטן קבוע:

int blockSize = 256;
int numBlocks = 128;                       // constant, not derived from n
vecAddStride<<<numBlocks, blockSize>>>(d_a, d_b, d_c, n);

הפלט:

n = 1048576, max error = 0.000000
c[0] = 3.0, c[1048575] = 3.0

שלב 4 - כמה איברים לכל thread. stride = 256 * 128 = 32768 threads. כל thread מטפל ב-n / stride = 1048576 / 32768 = 32 איברים בדיוק.

למה זה עבד: למרות שיש רק 32,768 threads לעומת 1,048,576 איברים, כל thread מדלג על פני המערך בצעדים בגודל stride ומכסה 32 איברים לאורך הלולאה. התנאי i < n משמש גם כשומר הגבולות, כך שאין צורך ב-if נפרד. כל האיברים כוסו, בלי חפיפה, כי הthreads מחלקים ביניהם את המרחב לפי שאריות (thread מספר k מטפל ב-k, k+stride, k+2·stride...).

איך להכליל: לולאת grid-stride היא התבנית המקצועית המומלצת: אותו kernel נכון לכל n, ואפשר לבחור numBlocks שממקסם תפוסה (occupancy) על החומרה - נהוג לכוון לכפולה של מספר ה-SMs (ב-H100 יש 132) - במקום להיגרר אחרי גודל הבעיה. בונוס: הרצה עם block יחיד (numBlocks = 1) היא דרך מצוינת לבדוק נכונות של kernel לפני שמגדילים, כי היא מסדרת (serializes) את העבודה.

פתרון תרגיל 5 - מדידת זמנים: CPU מול GPU

שלד המדידה (עם חימום לפני התזמון):

#include <chrono>

// --- CPU ---
auto t0 = std::chrono::high_resolution_clock::now();
for (int i = 0; i < n; i++) h_ref[i] = h_a[i] + h_b[i];
auto t1 = std::chrono::high_resolution_clock::now();
double cpu_ms = std::chrono::duration<double, std::milli>(t1 - t0).count();

// --- warmup ---
vecAdd<<<numBlocks, blockSize>>>(d_a, d_b, d_c, n);
CUDA_CHECK(cudaDeviceSynchronize());

// --- GPU: kernel only ---
cudaEvent_t s, e;
CUDA_CHECK(cudaEventCreate(&s));
CUDA_CHECK(cudaEventCreate(&e));
CUDA_CHECK(cudaEventRecord(s));
vecAdd<<<numBlocks, blockSize>>>(d_a, d_b, d_c, n);
CUDA_CHECK(cudaEventRecord(e));
CUDA_CHECK(cudaEventSynchronize(e));
float kernel_ms = 0.0f;
CUDA_CHECK(cudaEventElapsedTime(&kernel_ms, s, e));

טבלת תוצאות לדוגמה (H100 SXM; CPU מודרני יחיד-ליבה; זמנים במילישניות, ישתנו אצלכם):

        n     CPU (ms)   kernel (ms)   kernel+H2D+D2H (ms)
     1,024      0.002       0.006             0.09
    65,536      0.09        0.007             0.22
 1,048,576      1.4         0.021             2.6
16,777,216     22.0         0.28             39.5
134,217,728    178.0        2.20            310.0

התשובות:

  • kernel לבדו מנצח את ה-CPU כבר סביב n = 1 << 20 ומעלה: ב-n = 134M, ה-kernel לוקח כ-2.2 ms מול 178 ms ב-CPU - פי 80 מהיר.
  • כשמוסיפים את ההעברות המסקנה מתהפכת: kernel + H2D + D2H לוקח כ-310 ms - איטי יותר מה-CPU. ההעברות על ה-PCIe (בסדר גודל של עשרות GB/s) בולעות לחלוטין את הרווח, כי החישוב עצמו כה קל.

למה זה עבד: חיבור וקטורים הוא memory-bound: פעולת חיבור אחת מול 12 בתים של תעבורה לאיבר. זמן ה-kernel נשלט על ידי רוחב הפס של ה-HBM3 (ב-H100 בסדר גודל 3.35 TB/s), שגבוה בהרבה מרוחב הפס של זיכרון ה-CPU - ומכאן הניצחון בזמן ה-kernel לבד. אבל ההעברות ל-GPU וממנו עוברות ב-PCIe האיטי, וכשמשלמים עליהן בכל קריאה, הן צוואר הבקבוק. החימום חשוב כי ההרצה הראשונה כוללת JIT של PTX ל-SASS והקצאות עצלות שמזייפות את המדידה.

איך להכליל: ה-GPU משתלם כשעצימות החישוב לנתון גבוהה מספיק (compute-bound, לא חסום-העברה), או - וזה הנפוץ בפועל - כשהנתונים נשארים על ה-GPU לאורך שרשרת שלמה של פעולות (matmul, קונבולוציה, שכבות רשת), כך שההעברה משתלמת פעם אחת ומתחלקת על פני הרבה חישוב. את המסגרת הכמותית לזה - מודל ה-Roofline, arithmetic intensity ונקודת המפנה בין memory-bound לcompute-bound - נפתח בפרק על ביצועים.

פתרון תרגיל 6 (בונוס) - בדיקת רוחב פס

עבור n = 1 << 27 = 134,217,728 וזמן kernel לדוגמה kernel_ms = 2.20:

  • תעבורה כוללת: 12 בתים * n = 12 * 134217728 = 1,610,612,736 בתים (כ-1.5 GiB).
  • רוחב פס אפקטיבי: 1.610e9 בתים / (2.20e-3 שניות) / 1e9 ≈ 732 GB/s.

בהמרה מדויקת יותר לפי הזמן הזה מתקבל כ-0.73 TB/s. מול רוחב פס תיאורטי של כ-3.35 TB/s ב-H100 SXM, זה בערך 22% ניצול. (הערכים אצלכם ישתנו לפי הכרטיס והזמן שמדדתם.)

traffic      = 12 * 134217728        = 1.61e9  bytes
effective BW = 1.61e9 / 2.20e-3 / 1e9 ≈ 0.73   TB/s
utilization  = 0.73 / 3.35            ≈ 22%

למה זה עבד: ספרנו 12 בתים לאיבר (שני float נקראים, אחד נכתב) והכפלנו במספר האיברים, ואז חילקנו בזמן. התוצאה - מאות GB/s עד יחידות של TB/s - מאשרת שהחישוב אכן memory-bound ולא compute-bound: אנחנו רחוקים מ-66.9 TFLOPS השיא של ה-CUDA Cores, וקרובים לתקרת רוחב הפס. אחוז ניצול נמוך יחסית (עשרות אחוזים) נובע מכך שהמערך היחיד גדול מדי מכדי להיות תושב בcache, ומהoverhead של launch grid בודד.

איך להכליל: כשמדדתם רוחב פס אפקטיבי גבוה, זה סימן שהגישות מאוחדות (coalesced) - כלומר threads סמוכים ניגשים לכתובות סמוכות, וה-warp קורא בלוק רציף ביעילות. ה-kernel שלנו עושה זאת מטבעו כי thread i ניגש ל-a[i]. את המפה המלאה - מדוע coalescing חשוב, מהי נקודת המפנה בין memory-bound ל-compute-bound, ואיך מודל ה-Roofline קושר את רוחב הפס לarithmetic intensity - נלמד בפרק על ביצועים.