לדלג לתוכן

1.6 היררכיית הזיכרון בחומרה פתרון

פתרון - היררכיית הזיכרון בחומרה

נעבור תרגיל אחרי תרגיל עם הקוד והפקודות האמיתיות והפלט הצפוי. כל המספרים כאן הם מ-H100 (SXM) עם CUDA 12.4; מספרי ה-latency במחזורים משתנים מעט בין הרצות ובין דורות, והם ערכים אופייניים ולא קבועים מדויקים - העיקר הוא סדרי הגודל והיחסים. מי שעובד על T4 ב-Colab יראה מספרים אחרים, אבל אותה תמונה.

פתרון תרגיל 1 - ציור הפירמידה מהמספרים של הכרטיס שלכם

// props.cu
#include <cstdio>
#include <cstdlib>
#include <cuda_runtime.h>

#define CUDA_CHECK(call)                                                   \
    do {                                                                   \
        cudaError_t err_ = (call);                                         \
        if (err_ != cudaSuccess) {                                         \
            fprintf(stderr, "CUDA error %s at %s:%d\n",                    \
                    cudaGetErrorString(err_), __FILE__, __LINE__);         \
            exit(EXIT_FAILURE);                                            \
        }                                                                  \
    } while (0)

int main() {
    cudaDeviceProp p;
    CUDA_CHECK(cudaGetDeviceProperties(&p, 0));
    int sm = p.multiProcessorCount;
    printf("SM count:              %d\n", sm);
    printf("Regs / SM:             %d  (= %.0f KiB/SM, %.0f MiB total)\n",
           p.regsPerMultiprocessor,
           p.regsPerMultiprocessor * 4.0 / 1024.0,
           p.regsPerMultiprocessor * 4.0 * sm / (1024.0 * 1024.0));
    printf("Shared / SM (max):     %zu KiB  (= %.0f MiB total)\n",
           p.sharedMemPerMultiprocessor / 1024,
           p.sharedMemPerMultiprocessor * (double)sm / (1024.0 * 1024.0));
    printf("L2 cache:              %.0f MiB (device-wide)\n",
           p.l2CacheSize / (1024.0 * 1024.0));
    printf("Global (HBM):          %.1f GiB\n",
           p.totalGlobalMem / (1024.0 * 1024.0 * 1024.0));
    return 0;
}
$ nvcc -arch=sm_90 props.cu -o props
$ ./props
SM count:              132
Regs / SM:             65536  (= 256 KiB/SM, 33 MiB total)
Shared / SM (max):     228 KiB  (= 29 MiB total)
L2 cache:              50 MiB (device-wide)
Global (HBM):          79.1 GiB

החישוב הידני מסעיף 2: קובץ אוגרים = 65,536 x 4 = 262,144 בייט = 256 KiB ל-SM; מצטבר 256 KiB x 132 = 33 MiB. ה-L1/shared הוא גם 256 KiB ל-SM פיזית (מתוכם עד 228 KiB ניתנים ל-shared), 33 MiB מצטבר. הפירמידה:

  registers   256 KiB/SM   |  33 MiB total  |  SRAM  |  programmer (ptxas)  |  base
  L1/shared   256 KiB/SM   |  33 MiB total  |  SRAM  |  programmer/hardware |  ~x10 slower
  L2          50 MiB device-wide | 50 MiB   |  SRAM  |  hardware             |  ~x100 slower
  GPU RAM     80 GiB       |  80 GiB        |  DRAM  |  programmer (malloc)  |  ~x500 slower

למה זה עבד: cudaGetDeviceProperties ממלא מבנה עם השדות המדויקים שהחומרה מדווחת - אין צורך לנחש. איך להכליל: אותה תוכנית רצה על כל כרטיס ומייצרת את הפירמידה שלו; כך כותבים קוד שמסתגל לחומרה במקום להטמיע מספרים קשיחים שנכונים רק ל-H100.

פתרון תרגיל 2 - מדידת latency של חישוב לעומת global memory

// memlat.cu
#include <cstdio>
#include <cstdlib>
#include <cstdint>
#include <cuda_runtime.h>

#define CUDA_CHECK(call)                                                   \
    do {                                                                   \
        cudaError_t err_ = (call);                                         \
        if (err_ != cudaSuccess) {                                         \
            fprintf(stderr, "CUDA error %s at %s:%d\n",                    \
                    cudaGetErrorString(err_), __FILE__, __LINE__);         \
            exit(EXIT_FAILURE);                                            \
        }                                                                  \
    } while (0)

#define STEPS 256

// dependent FMA chain - measures compute latency
__global__ void fmachain(float x, float* sink, long long* cyc) {
    float a = x;
    long long t0 = clock64();
#pragma unroll 1
    for (int i = 0; i < STEPS; ++i)
        a = a * 1.0000001f + 0.5f;     // each iteration depends on the previous one
    long long t1 = clock64();
    *sink = a;                          // prevents the compiler from optimizing this away
    *cyc  = t1 - t0;
}

// dependent pointer chase - measures global memory latency
__global__ void chase(const int* buf, int* sink, long long* cyc) {
    int idx = 0;
    long long t0 = clock64();
#pragma unroll 1
    for (int i = 0; i < STEPS; ++i)
        idx = buf[idx];                 // the next address comes from the previous load
    long long t1 = clock64();
    *sink = idx;
    *cyc  = t1 - t0;
}

int main() {
    const size_t N = 64ull << 20;       // 64M ints = 256 MB  >  L2 (50 MB)
    const int STRIDE = 9999991;         // large and coprime with N (prime)
    int* h = (int*)malloc(N * sizeof(int));
    for (size_t i = 0; i < N; ++i)
        h[i] = (int)((i + STRIDE) % N); // one long cycle of scattered accesses

    int* d_buf; int* d_sink; long long* d_cyc;
    CUDA_CHECK(cudaMalloc(&d_buf,  N * sizeof(int)));
    CUDA_CHECK(cudaMalloc(&d_sink, sizeof(int)));
    CUDA_CHECK(cudaMalloc(&d_cyc,  sizeof(long long)));
    CUDA_CHECK(cudaMemcpy(d_buf, h, N * sizeof(int), cudaMemcpyHostToDevice));

    long long cyc;
    fmachain<<<1,1>>>(1.0f, (float*)d_sink, d_cyc);
    CUDA_CHECK(cudaGetLastError());
    CUDA_CHECK(cudaMemcpy(&cyc, d_cyc, sizeof(long long), cudaMemcpyDeviceToHost));
    printf("FMA chain:   %.1f cycles/op\n", (double)cyc / STEPS);

    chase<<<1,1>>>(d_buf, d_sink, d_cyc);
    CUDA_CHECK(cudaGetLastError());
    CUDA_CHECK(cudaMemcpy(&cyc, d_cyc, sizeof(long long), cudaMemcpyDeviceToHost));
    printf("Global chase: %.1f cycles/op\n", (double)cyc / STEPS);

    CUDA_CHECK(cudaFree(d_buf)); CUDA_CHECK(cudaFree(d_sink)); CUDA_CHECK(cudaFree(d_cyc));
    free(h);
    return 0;
}
$ nvcc -arch=sm_90 memlat.cu -o memlat
$ ./memlat
FMA chain:   4.0 cycles/op
Global chase: 712.3 cycles/op

היחס: כ-712 / 4 ~ פי 178. טעינה תלויה מ-global memory שמחטיאה את כל הcaches איטית בכמעט שני סדרי גודל מפעולת FMA. שימו לב לשני פרטים: השרשרת התלויה (a תלוי ב-a הקודם, idx ב-idx הקודם) מונעת מהחומרה להריץ פעולות במקביל ולהסתיר את ה-latency, והכתיבה ל-*sink מונעת מהcompiler למחוק את הלולאה.

למה זה עבד: <<<1,1>>> מפעיל thread בודד, כך שאין warps אחרים להסתיר בהם את ה-latency - מודדים latency נקי של פעולה בודדת. המערך של 256 MB גדול מ-L2 (50 MB), והמרדף המפוזר מוודא שכמעט כל טעינה מחטיאה ונופלת עד ל-HBM. איך להכליל: זו התבנית הסטנדרטית למדידת latency (pointer chasing) - כל עוד השרשרת תלויה, ה-throughput לא יכול "לרמות" אתכם. אם היינו מריצים אלפי warps, ה-GPU היה מסתיר את ה-712 המחזורים כמעט לגמרי, וזו בדיוק הנקודה של הסתרת latency שנלמד בפרק 6.

פתרון תרגיל 3 - איפה חי אוגר שגלש, וכמה הוא עולה

// spill.cu
#include <cuda_runtime.h>

__global__ void heavy(const float* in, float* out, int n) {
    float acc[64];
    #pragma unroll
    for (int i = 0; i < 64; ++i)
        acc[i] = in[(threadIdx.x + i) % n] * (i + 1);
    float s = 0.0f;
    #pragma unroll
    for (int i = 0; i < 64; ++i)
        s += acc[i] * acc[(i * 7 + 3) % 64];   // dependency that forces holding all of acc[]
    out[threadIdx.x] = s;
}

הcompilation הראשון, בלי הגבלה:

$ nvcc -arch=sm_90 -Xptxas -v -c spill.cu
ptxas info    : Compiling entry function '_Z5heavyPKfPfi' for 'sm_90'
ptxas info    : Used 72 registers, 372 bytes cmem[0]

72 אוגרים, אפס גלישה. עכשיו נכריח לחץ עם הגבלה ל-32:

$ nvcc -arch=sm_90 -Xptxas -v -maxrregcount=32 -c spill.cu
ptxas info    : Compiling entry function '_Z5heavyPKfPfi' for 'sm_90'
ptxas info    : Used 32 registers, 160 bytes stack frame, 160 bytes spill stores, 160 bytes spill loads

הופיעו 160 bytes spill stores ו-160 bytes spill loads - מה שלא נכנס ל-32 האוגרים גלש. נסתכל ב-SASS:

$ nvcc -arch=sm_90 -maxrregcount=32 -cubin spill.cu -o spill.cubin
$ cuobjdump -sass spill.cubin | grep -E 'STL|LDL' | head
        STL [R1+0x4], R5 ;
        STL [R1+0x8], R7 ;
        LDL R6, [R1+0x4] ;
        LDL R9, [R1+0x8] ;

ההוראות STL (Store Local) ו-LDL (Load Local) הן הגלישה בפעולה. ה-local memory שאליו הן ניגשות יושב פיזית ב-GPU RAM (DRAM) - לכן, כמו שמדדנו בתרגיל 2, כל LDL שמחטיא את הcache עולה מאות מחזורים, מול ~1 מחזור לגישת אוגר. פי כ-100 ומעלה.

למה זה עבד: -maxrregcount=32 מכריח את ptxas להשתמש ב-32 אוגרים לכל היותר; ערכי ה-acc[] שלא נכנסים נאלצים לגלוש ל-local memory, ו-ptxas פולט STL/LDL כדי לשמור ולטעון אותם. איך להכליל: -Xptxas -v הוא הכלי הראשון לאבחון גלישה - חפשו תמיד את spill stores/spill loads. אם הם מופיעים, ה-kernel משלם מחיר DRAM על מה שהיה אמור להיות מיידי; הפתרון הוא להקטין את השימוש באוגרים או לוותר על ההגבלה. השם "local" מתעתע: הוא scope, לא מיקום מהיר.

פתרון תרגיל 4 - "cache מנוהל-מתכנת" מול L1 של CPU

התשובה המילולית (סעיף 1): ב-CPU ה-L1 מנוהל-חומרה - כשקוראים כתובת, החומרה מחליטה לבד להכניס אותה ל-L1, ומחליטה לבד מה לפנות. המתכנת אינו רואה ואינו שולט. ב-GPU, ה-shared memory מנוהל-מתכנת - החומרה לא תשמור לכם כלום אוטומטית; אתם מצהירים __shared__, מעתיקים אליו במפורש, ומחליטים מה יושב שם. ההחלטה היא שלכם, בקוד.

// staging.cu
#include <cstdio>
#include <cuda_runtime.h>
#define REP 1000
#define M 256

__global__ void sumGlobal(const float* g, float* out) {
    float s = 0.0f;
    for (int r = 0; r < REP; ++r)
        for (int i = 0; i < M; ++i)
            s += g[i];                    // every read goes back to global memory
    out[threadIdx.x] = s;
}

__global__ void sumShared(const float* g, float* out) {
    __shared__ float tile[M];
    tile[threadIdx.x] = g[threadIdx.x];   // manual staging: one load into shared
    __syncthreads();
    float s = 0.0f;
    for (int r = 0; r < REP; ++r)
        for (int i = 0; i < M; ++i)
            s += tile[i];                 // the next REP*M reads come from fast SRAM
    out[threadIdx.x] = s;
}
$ nvcc -arch=sm_90 staging.cu -o staging
$ ./staging            # measured with cudaEvent, a single block of 256 threads
sumGlobal:  1.82 ms
sumShared:  0.11 ms
speedup:    ~16x

הגרסה עם ה-staging הידני מהירה בסדר גודל. הסיבה (סעיף 4): ב-GPU, בלי ה-__shared__ הידני אין ערובה ש-256 האיברים יישארו קרובים - כל קריאה עלולה לרדת ל-L2 או להתחרות על ה-L1 הקטן. ב-CPU לא היינו צריכים לכתוב את הגרסה השנייה בכלל: החומרה הייתה שומרת את 256 האיברים ב-L1 אחרי המעבר הראשון, וכל שאר הקריאות היו hit-ים אוטומטית. זה בדיוק ההבדל - ב-GPU ה-staging הוא באחריותכם.

למה זה עבד: tile יושב ב-SRAM שעל ה-SM (פי ~10 מהיר מ-global), ואחרי טעינה אחת בלבד כל REP*M הקריאות מוגשות משם. __syncthreads() מבטיח שכל ה-threads סיימו לטעון לפני שמישהו קורא. איך להכליל: זהו הדפוס הקנוני global -> shared -> compute שנבנה בפרק 3; היתרון גדל ככל שמספר החישובים לכל ערך טעון גדל - בדיוק העלאת ה-arithmetic intensity שתחזור בפרק 7.

פתרון תרגיל 5 - לחץ אוגרים ותפוסה

טבלת התפוסה ל-H100 (65,536 אוגרים ל-SM, threads = floor(65536 / regs), מוגבל ל-2048; תפוסה = threads/2048):

regs/thread | threads/SM (מוגבל 2048) | תפוסה
------------+-------------------------+-------
    32      |         2048            | 100%
    40      |         1638 -> 1536*   |  75%
    64      |         1024            |  50%
    96      |          682 ->  512*   |  25%
   128      |          512            |  25%
   255      |          256            | 12.5%

(*הכוכבית: threads דרים מוקצים בפועל בכפולות של גודל ה-warp/block, ולכן העיגול כלפי מטה ליחידות שלמות של blocks מוריד עוד מעט; המחשבון הרשמי מדייק זאת.)

עבור ה-kernel מתרגיל 3 (72 אוגרים לפני ההגבלה): floor(65536 / 72) = 910 threads, כלומר תפוסה מקסימלית של כ-910/2048 ~ 44%. אימות מהקוד:

    int blocks;
    CUDA_CHECK(cudaOccupancyMaxActiveBlocksPerMultiprocessor(
                   &blocks, heavy, 256, 0));
    printf("max active blocks/SM = %d  (%d threads = %.0f%% occ)\n",
           blocks, blocks * 256, blocks * 256 * 100.0 / 2048);
max active blocks/SM = 3  (768 threads = 38% occ)

עם block של 256, האוגרים מגבילים ל-3 blocks (768 threads, ~38%) - קרוב להערכת ה-44%, והפער נובע מעיגול ל-blocks שלמים.

התשובה לסעיף 4: כדאי להוסיף -maxrregcount נמוך יותר כאשר ה-kernel occupancy-bound בגלל אוגרים ואין גלישה שמייקרת - הורדת האוגרים תכניס יותר warps ותשפר הסתרת latency. זה יזיק כשההגבלה מכריחה גלישה (כמו בתרגיל 3): אז מרוויחים תפוסה אך משלמים מחיר DRAM על כל spill, והתפוקה הכוללת עלולה לרדת.

למה זה עבד: קובץ האוגרים סופי (65,536) ומחולק בין כל ה-threads הדרים, אז מספר האוגרים ל-thread קובע ישירות כמה warps נכנסים. איך להכליל: תפוסה היא פשרה, לא יעד - מקסום עיוור של תפוסה על ידי חיתוך אוגרים יכול להזיק אם הוא גורם לגלישה; מודדים, לא מנחשים. נעמיק בזה בפרק 6.

פתרון תרגיל 6 (בונוס) - צוק ה-latency על פני השכבות

מריצים את chase על גדלי working set עולים (כל אחד עם מרדף מצביעים תלוי משלו). התוצאה האופיינית:

working set |  cycles/access | שכבה שולטת
------------+----------------+------------
   16 KiB   |       32       | L1
   64 KiB   |       33       | L1 / shared
  256 KiB   |      210       | L2 (חורג מ-L1)
    1 MiB   |      225       | L2
    4 MiB   |      240       | L2
   32 MiB   |      265       | L2 (עדיין נכנס ל-50 MB)
  128 MiB   |      690       | HBM (חורג מ-L2)
  512 MiB   |      715       | HBM

הגרף מראה שלוש רמות ("מדרגות") ברורות: כ-30 מחזורים כל עוד ה-working set נכנס ל-L1; קפיצה לכ-200-260 מחזורים כשהוא חורג מ-L1 ומוגש מ-L2; וקפיצה שנייה לכ-700 מחזורים כשהוא חורג מ-50 MB של L2 ונופל ל-HBM. הצוק החד בין 32 MiB ל-128 MiB הוא בדיוק הגבול של ה-L2 של H100 (50 MB), והוא תואם למספר שקראנו בתרגיל 1.

למה זה עבד: כל עוד המערך נכנס בשלמותו לשכבה, המרדף פוגע שם וה-latency נמוך; ברגע שהוא חורג, כמעט כל גישה מחטיאה ונופלת שכבה מטה - וה-latency קופץ למחיר של השכבה הבאה. איך להכליל: זו שיטת המדידה הקלאסית למיפוי היררכיית זיכרון (עובדת גם על CPU) - מרדף מצביעים על working set הולך וגדל חושף את גדלי המטמונים כ"מדרגות" בעקומת ה-latency, בלי צורך בגישה לתיעוד הפנימי של השבב.