3.4 זיכרון קבוע וטקסטורה פתרון
פתרון - זיכרון קבוע וטקסטורה¶
הפתרונות כאן מלאים ורצים כפי שהם. הזמנים, יחסי המהירות והפלטים המדויקים תלויים בכרטיס, בדרייבר ובגרסת CUDA - המספרים שמובאים הם לדוגמה מהרצה על H100 SXM, וישתנו אצלכם. מה שלא ישתנה הוא היחסים האיכותיים: broadcast מנצח כתובות מתפצלות, const __restrict__ מפיק LDG.E.CI, והאינטרפולציה הבי-לינארית מחזירה ממוצעים צפויים.
פתרון תרגיל 1 - מקדמי סטנסיל בזיכרון קבוע מול global¶
הקובץ stencil.cu המלא:
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <cuda_runtime.h>
#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)
#define RADIUS 4
#define DIAM (2 * RADIUS + 1) // 9
__constant__ float c_coef[DIAM];
__global__ void stencil_const(const float* in, float* out, int n) {
int i = blockIdx.x * blockDim.x + threadIdx.x;
if (i >= RADIUS && i < n - RADIUS) {
float acc = 0.0f;
#pragma unroll
for (int k = -RADIUS; k <= RADIUS; ++k)
acc += c_coef[k + RADIUS] * in[i + k];
out[i] = acc;
}
}
__global__ void stencil_global(const float* in, const float* coef,
float* out, int n) {
int i = blockIdx.x * blockDim.x + threadIdx.x;
if (i >= RADIUS && i < n - RADIUS) {
float acc = 0.0f;
#pragma unroll
for (int k = -RADIUS; k <= RADIUS; ++k)
acc += coef[k + RADIUS] * in[i + k];
out[i] = acc;
}
}
int main(void) {
int n = 1 << 26; // ~67M
size_t bytes = (size_t)n * sizeof(float);
float h_coef[DIAM];
for (int k = 0; k < DIAM; ++k) h_coef[k] = 1.0f / DIAM;
float* h_in = (float*)malloc(bytes);
for (int i = 0; i < n; ++i) h_in[i] = 1.0f;
float *d_in, *d_out, *d_coef;
CUDA_CHECK(cudaMalloc(&d_in, bytes));
CUDA_CHECK(cudaMalloc(&d_out, bytes));
CUDA_CHECK(cudaMalloc(&d_coef, sizeof(h_coef)));
CUDA_CHECK(cudaMemcpy(d_in, h_in, bytes, cudaMemcpyHostToDevice));
CUDA_CHECK(cudaMemcpyToSymbol(c_coef, h_coef, sizeof(h_coef))); // -> constant memory
CUDA_CHECK(cudaMemcpy(d_coef, h_coef, sizeof(h_coef), // -> global
cudaMemcpyHostToDevice));
int bs = 256, gs = (n + bs - 1) / bs;
cudaEvent_t s, e; float ms;
CUDA_CHECK(cudaEventCreate(&s)); CUDA_CHECK(cudaEventCreate(&e));
stencil_const<<<gs, bs>>>(d_in, d_out, n); // warm-up
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaEventRecord(s));
stencil_const<<<gs, bs>>>(d_in, d_out, n);
CUDA_CHECK(cudaEventRecord(e)); CUDA_CHECK(cudaEventSynchronize(e));
CUDA_CHECK(cudaEventElapsedTime(&ms, s, e));
printf("constant : %.3f ms\n", ms);
CUDA_CHECK(cudaEventRecord(s));
stencil_global<<<gs, bs>>>(d_in, d_coef, d_out, n);
CUDA_CHECK(cudaEventRecord(e)); CUDA_CHECK(cudaEventSynchronize(e));
CUDA_CHECK(cudaEventElapsedTime(&ms, s, e));
printf("global : %.3f ms\n", ms);
// verification: uniform input 1.0, coefficient sum 1.0 -> output 1.0
float* h_out = (float*)malloc(bytes);
CUDA_CHECK(cudaMemcpy(h_out, d_out, bytes, cudaMemcpyDeviceToHost));
printf("out[100] = %.4f (expected 1.0000)\n", h_out[100]);
cudaFree(d_in); cudaFree(d_out); cudaFree(d_coef);
free(h_in); free(h_out);
return 0;
}
הרצה:
פלט לדוגמה (H100):
למה זה עבד: בלולאה הפנימית האינדקס k + RADIUS אחיד על פני כל ה-warp, ולכן c_coef[k+RADIUS] הוא אותה כתובת ל-32 הthreads - תבנית ה-broadcast של זיכרון קבוע, שמשרתת את כל ה-warp בקריאה אחת. גרסת ה-global קוראת את אותם מקדמים דרך נתיב ה-L1 הרגיל; מכיוון שרק 9 מקדמים נקראים שוב ושוב הם נשארים חמים ב-L1, ולכן ההפרש קטן - הסטנסיל memory-bound על קריאת מערך הקלט in, לא על המקדמים. איך להכליל: זיכרון קבוע נותן את היתרון הגדול ביותר כשהערך המשודר נקרא הרבה יחסית לשאר עבודת הזיכרון, וכשהcache הרגיל היה נדחק. הרווח האמיתי של broadcast מודגם בתרגיל 4, שם הכתובות מתפצלות ואין נתיב אחר.
פתרון תרגיל 2 - קלט לקריאה-בלבד ובדיקת SASS¶
מוסיפים גרסה עם __restrict__ ומוודאים את הנתיב ב-SASS. ה-kernel:
__global__ void stencil_ro(const float* __restrict__ in,
float* __restrict__ out, int n) {
int i = blockIdx.x * blockDim.x + threadIdx.x;
if (i >= RADIUS && i < n - RADIUS) {
float acc = 0.0f;
#pragma unroll
for (int k = -RADIUS; k <= RADIUS; ++k)
acc += c_coef[k + RADIUS] * __ldg(&in[i + k]); // explicit, to confirm the path
out[i] = acc;
}
}
חילוץ ה-SASS:
nvcc -O2 -arch=sm_90a -cubin -o stencil.cubin stencil.cu
cuobjdump -sass stencil.cubin | grep -E "LDG|LDC"
פלט לדוגמה (מקוצר):
/*0090*/ LDG.E.CI R6, [R4.64] ; # stencil_ro: input via the read-only cache
/*00a0*/ LDG.E.CI R8, [R4.64+0x4] ;
...
/*0110*/ LDG.E R6, [R4.64] ; # stencil_const: input on the regular L1 path
/*0120*/ LDG.E R8, [R4.64+0x4] ;
המקדמים מזיכרון הקבוע אינם LDG כלל - הם נטענים כאופרנד קבוע. אפשר לראות זאת עם:
פלט לדוגמה:
/*0098*/ FFMA R7, R6, c[0x3][0x0], R7 ; # coefficient c_coef[0] as an operand from constant bank 0x3
למה זה עבד: const __restrict__ (או __ldg מפורש) מבטיח לcompiler שהמצביע in נקרא-בלבד ואינו aliased, ולכן ה-ptxas חופשי להנפיק LDG.E.CI - טעינה דרך cache הקריאה-בלבד הלא-קוהרנטי (הסיומת .CI = Cache in Invalidate). המקדמים מזיכרון הקבוע נטמעים כאופרנד c[0x3][...] בתוך פקודת ה-FFMA עצמה (כי הם משודרים אחיד), ולכן אין להם LDG נפרד. איך להכליל: כל מצביע קלט שנקרא-בלבד ב-kernel ראוי ל-const __restrict__ - זה גם מרמז על cache הקריאה-בלבד וגם מתיר לcompiler לתזמן טעינות אגרסיבי יותר. קריאת SASS עם cuobjdump -sass היא הדרך המוסמכת לוודא איזה נתיב זיכרון קיבלתם בפועל, בכל אופטימיזציה - נשתמש בה שוב בפרק 4 ובפרק 8.
פתרון תרגיל 3 - דגימת תמונה עם אינטרפולציה בי-לינארית בחומרה¶
הקובץ resample.cu:
#include <cstdio>
#include <cstdlib>
#include <cuda_runtime.h>
#define CUDA_CHECK(call) \
do { cudaError_t e_=(call); if(e_!=cudaSuccess){ \
fprintf(stderr,"CUDA %s:%d %s -> %s\n",__FILE__,__LINE__,#call, \
cudaGetErrorString(e_)); exit(1);} } while(0)
__global__ void sample(cudaTextureObject_t tex, float* out,
const float2* pts, int m) {
int i = blockIdx.x * blockDim.x + threadIdx.x;
if (i < m) out[i] = tex2D<float>(tex, pts[i].x, pts[i].y);
}
int main(void) {
const int W = 4, H = 4;
float h_img[W * H];
for (int y = 0; y < H; ++y)
for (int x = 0; x < W; ++x)
h_img[y * W + x] = (float)(x + y); // known gradient
cudaChannelFormatDesc ch = cudaCreateChannelDesc<float>();
cudaArray_t arr;
CUDA_CHECK(cudaMallocArray(&arr, &ch, W, H));
CUDA_CHECK(cudaMemcpy2DToArray(arr, 0, 0, h_img,
W * sizeof(float),
W * sizeof(float), H,
cudaMemcpyHostToDevice));
cudaResourceDesc resDesc = {};
resDesc.resType = cudaResourceTypeArray;
resDesc.res.array.array = arr;
cudaTextureDesc texDesc = {};
texDesc.addressMode[0] = cudaAddressModeClamp;
texDesc.addressMode[1] = cudaAddressModeClamp;
texDesc.filterMode = cudaFilterModeLinear;
texDesc.readMode = cudaReadModeElementType;
texDesc.normalizedCoords = 0;
cudaTextureObject_t tex = 0;
CUDA_CHECK(cudaCreateTextureObject(&tex, &resDesc, &texDesc, NULL));
// texel center (px,py) = (px+0.5, py+0.5)
float2 h_pts[3] = { {1.5f, 1.5f}, // center of texel (1,1) -> 2.0
{1.0f, 1.0f}, // corner of 4 texels (0,0),(1,0),(0,1),(1,1) -> average 1.0
{2.0f, 1.5f} }; // midpoint between texels (1,1)=2 and (2,1)=3 -> 2.5
int m = 3;
float2* d_pts; float* d_out;
CUDA_CHECK(cudaMalloc(&d_pts, sizeof(h_pts)));
CUDA_CHECK(cudaMalloc(&d_out, m * sizeof(float)));
CUDA_CHECK(cudaMemcpy(d_pts, h_pts, sizeof(h_pts), cudaMemcpyHostToDevice));
sample<<<1, 32>>>(tex, d_out, d_pts, m);
CUDA_CHECK(cudaGetLastError());
float h_out[3];
CUDA_CHECK(cudaMemcpy(h_out, d_out, sizeof(h_out), cudaMemcpyDeviceToHost));
printf("sample(1.5,1.5) = %.4f (expect 2.0)\n", h_out[0]);
printf("sample(1.0,1.0) = %.4f (expect 1.0)\n", h_out[1]);
printf("sample(2.0,1.5) = %.4f (expect 2.5)\n", h_out[2]);
cudaDestroyTextureObject(tex);
cudaFreeArray(arr);
cudaFree(d_pts); cudaFree(d_out);
return 0;
}
הרצה:
פלט לדוגמה:
sample(1.5,1.5) = 2.0000 (expect 2.0)
sample(1.0,1.0) = 1.0000 (expect 1.0)
sample(2.0,1.5) = 2.5000 (expect 2.5)
למה זה עבד: עם cudaFilterModeLinear וקואורדינטות לא-מנורמלות, מרכזי הטקסלים יושבים ב-integer + 0.5. דגימה במרכז טקסל (1.5,1.5) נופלת בדיוק על טקסל אחד ומחזירה את ערכו (2.0). דגימה ב-(1.0,1.0) נופלת בדיוק בפינה בין ארבעה טקסלים והחומרה מחזירה את השקלול הבי-לינארי שלהם, שהוא הממוצע (0+1+1+2)/4 = 1.0. דגימה ב-(2.0,1.5) נופלת באמצע בין שני טקסלים אופקיים ומחזירה את ממוצעם (2+3)/2 = 2.5. כל השקלולים נעשו ביחידת החומרה בפעולה אחת לכל דגימה, בלי חשבון ידני. איך להכליל: אותה תבנית - cudaArray + cudaTextureObject_t + cudaFilterModeLinear - היא הבסיס לכל דגימה-מחדש, סיבוב, עיוות (warping) והגדלה/הקטנה של תמונות ב-GPU. זכרו שהמשקלים מחושבים בדיוק שבור של 8 ביטים, כך שלשימוש מדעי הדורש דיוק מלא עדיף לחשב את האינטרפולציה ידנית מ-global/shared.
פתרון תרגיל 4 - מתי זיכרון קבוע מזיק¶
ה-kernels וה-main המרכזי:
__constant__ float c_table[256];
__global__ void read_broadcast(float* out, int n) {
int i = blockIdx.x * blockDim.x + threadIdx.x;
if (i < n) out[i] = c_table[0]; // the whole warp: same address
}
__global__ void read_divergent(float* out, int n) {
int i = blockIdx.x * blockDim.x + threadIdx.x;
if (i < n) out[i] = c_table[threadIdx.x & 255]; // each thread: a different address
}
int main(void) {
int n = 1 << 26;
float h_table[256];
for (int k = 0; k < 256; ++k) h_table[k] = (float)k;
CUDA_CHECK(cudaMemcpyToSymbol(c_table, h_table, sizeof(h_table)));
float* d_out;
CUDA_CHECK(cudaMalloc(&d_out, (size_t)n * sizeof(float)));
int bs = 256, gs = (n + bs - 1) / bs;
cudaEvent_t s, e; float ms;
CUDA_CHECK(cudaEventCreate(&s)); CUDA_CHECK(cudaEventCreate(&e));
read_broadcast<<<gs, bs>>>(d_out, n); // warm-up
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaEventRecord(s));
read_broadcast<<<gs, bs>>>(d_out, n);
CUDA_CHECK(cudaEventRecord(e)); CUDA_CHECK(cudaEventSynchronize(e));
CUDA_CHECK(cudaEventElapsedTime(&ms, s, e));
printf("broadcast : %.3f ms\n", ms);
CUDA_CHECK(cudaEventRecord(s));
read_divergent<<<gs, bs>>>(d_out, n);
CUDA_CHECK(cudaEventRecord(e)); CUDA_CHECK(cudaEventSynchronize(e));
CUDA_CHECK(cudaEventElapsedTime(&ms, s, e));
printf("divergent : %.3f ms\n", ms);
cudaFree(d_out);
return 0;
}
פלט לדוגמה (H100):
למה זה עבד: ב-read_broadcast כל 32 הthreads של warp מבקשים את c_table[0] - כתובת אחת, טרנזקציה אחת מזיכרון הקבוע, מהיר כמו register. ב-read_divergent threads ה-warp מבקשים 32 כתובות שונות (threadIdx.x & 255), וחומרת הזיכרון הקבוע מסדרת אותן בטור - עד 32 טרנזקציות ל-warp. הפער בזמן (כאן פי ~3, ולא פי 32 מלא כי הקצה השני חסום-כתיבה על out) מדגים שהזמן פרופורציוני למספר הכתובות השונות. אילו c_table היה ב-global memory, דווקא read_divergent היה מהיר יותר, כי 32 הכתובות הרציפות של warp מתאחדות (coalesce) לטרנזקציה אחת. איך להכליל: בחרו זיכרון קבוע רק כשהאינדקס אליו אחיד על פני ה-warp (מקדמים, פרמטרים). טבלת חיפוש שנקראת באינדקסים תלויי-threadIdx שייכת ל-global (מאוחד) או ל-texture - שם גישה מתפצלת אך מרחבית-מקומית נהנית מהcache במקום להיענש.
פתרון תרגיל 5 (בונוס) - סטנסיל 2D: טקסטורה מול global¶
ה-kernels המרכזיים (ה-main מקים texture כמו בתרגיל 3, על תמונה 4096 x 4096):
__global__ void jacobi_tex(cudaTextureObject_t tex, float* out, int W, int H) {
int x = blockIdx.x * blockDim.x + threadIdx.x;
int y = blockIdx.y * blockDim.y + threadIdx.y;
if (x < W && y < H) {
float xc = x + 0.5f, yc = y + 0.5f;
float v = tex2D<float>(tex, xc, yc)
+ tex2D<float>(tex, xc - 1.f, yc) // clamp handles the boundary
+ tex2D<float>(tex, xc + 1.f, yc)
+ tex2D<float>(tex, xc, yc - 1.f)
+ tex2D<float>(tex, xc, yc + 1.f);
out[y * W + x] = 0.2f * v;
}
}
__global__ void jacobi_glob(const float* __restrict__ in, float* out,
int W, int H) {
int x = blockIdx.x * blockDim.x + threadIdx.x;
int y = blockIdx.y * blockDim.y + threadIdx.y;
if (x < W && y < H) {
int i = y * W + x;
float c = in[i];
float l = (x > 0) ? in[i - 1] : c; // manual boundary checks
float r = (x < W - 1) ? in[i + 1] : c;
float u = (y > 0) ? in[i - W] : c;
float d = (y < H - 1) ? in[i + W] : c;
out[i] = 0.2f * (c + l + r + u + d);
}
}
הlaunch עם grid דו-ממדי:
dim3 block(16, 16);
dim3 grid((W + 15) / 16, (H + 15) / 16);
jacobi_tex<<<grid, block>>>(tex, d_out, W, H);
פלט לדוגמה (H100, 4096 x 4096):
texture : 0.283 ms (eff. BW ~1423 GB/s)
global : 0.361 ms (eff. BW ~1115 GB/s)
max diff (interior) = 0.000000
למה זה עבד: בסטנסיל הדו-ממדי, השכנים האנכיים (in[i - W], in[i + W]) רחוקים W איברים בזיכרון, ולכן בגרסת ה-global הם נופלים בשורות cache שונות וגורמים ל-cache misses. cache הטקסטורה, שממטמן tiles דו-ממדיים, שומר את השכנים האנכיים והאופקיים קרוב, ולכן משיג רוחב פס אפקטיבי גבוה יותר; בנוסף cudaAddressModeClamp מבטל את ענפי בדיקות-הגבול (x > 0 וכו') שגורמים לפיצול קל. חישוב רוחב הפס: 24 בתים * 4096^2 / (0.283e-3 s) / 1e9 = ~1423 GB/s. איך להכליל: ל-texture יתרון בסטנסילים ובעיבוד תמונה בזכות הלוקליות הדו-ממדית והטיפול-הגבולות החינמי, אבל בחומרה מודרנית גרסת shared-memory tiling (טעינת tile (16+2) x (16+2) פעם אחת ל-block, כפי שראינו ב-3.3) בדרך כלל מנצחת את שתיהן, כי היא ממחזרת כל טעינה בין כל threads ה-block ונמנעת מ-cache misses לגמרי. הבחירה הסופית - texture, global מאוחד, או shared-tiled - היא בדיוק סוג ההחלטה שתקבלו בפרויקט 3.5 ותכמתו בפרק 8.