100 DAYS OF GPU Course preview
Day 7Module 1
in-technical-review

Two-dimensional grids and image kernels

Someone reading NVIDIA's programming guide, on the MatAdd example every 2D CUDA tutorial descends from:

"Here, i is associated with x as if the matrix was transposed, since A[i][j] will skip one 'row' for i and increment for j."

https://forums.developer.nvidia.com/t/confusion-regarding-cuda-2d-indexing-found-in-the-official-programming-guide/328107 (checked 2026-08-30)

The kernel they were reading, as that thread quotes it:

__global__ void MatAdd(float A[N][N], float B[N][N], float C[N][N]) {
  int i = threadIdx.x;
  int j = threadIdx.y;
  C[i][j] = A[i][j] + B[i][j];
}

The reader is right. The reply is also right: either mapping can work if its memory access pattern is sound.

Both mappings give the same answer for a square matrix add. For a 611 by 397 image, one writes the wrong pixels while the other lets each warp touch consecutive bytes.

This page shows how to launch a 2D grid over an image, map each axis, and guard both rounded-up edges.

A picture is one flat array, and one formula makes it two-dimensional

unsigned char* gray points to a flat byte array in global memory. Row-major order stores all of row 0 first, then row 1, and so on. The pixel at (row, col) sits at

row * width + col

The pointer type does not store the image shape. A wrong formula may still produce valid addresses, so CUDA does not report an error.

A grid can have two dimensions. dim3 stores three unsigned values and defaults missing values to 1, so dim3 block(32, 8) means 32 by 8 by 1.

The launch configuration takes one dim3 for the grid and one for the block. The guide states: "The use of multi-dimensional thread blocks and grids is for convenience only and does not affect performance." (https://docs.nvidia.com/cuda/cuda-programming-guide/02-basics/writing-cuda-kernels.html , section 2.3.2 "Thread Hierarchy", checked 2026-08-30.)

Every built-in index variable grows a .y, and the global thread index from day 4 is written once per axis:

col = blockIdx.x * blockDim.x + threadIdx.x;
row = blockIdx.y * blockDim.y + threadIdx.y;

x carries the column and y carries the row. That is not a house style. The same guide section says why:

"The threads of a block are linearized predictably: the first index x moves the fastest, followed by y and then z. This means that in the linearization of a thread indices, consecutive values of threadIdx.x indicate consecutive threads, threadIdx.y has a stride of blockDim.x, and threadIdx.z has a stride of blockDim.x * blockDim.y. This affects how threads are assigned to warps."

CUDA forms a warp from consecutive threads. In a 32 by 8 block, lanes 0 to 31 of the first warp use threadIdx.x values 0 to 31 while threadIdx.y stays 0.

Map col to x and the warp covers 32 adjacent pixels in one row, which are 32 consecutive output bytes. Map col to y and the warp covers 32 pixels in one column, each width bytes apart. Day 11 measures the cost.

The intuition that survives a square matrix

In mathematics and nested loops, A[i][j] uses row i and column j. But threadIdx.x appearing first in threadIdx does not mean it must map to i.

The guide's MatAdd combines those two conventions, which caused the forum reader's confusion.

A square matrix hides the mistake. Both mappings visit each element once, so an elementwise operation gives the right result. A rectangular buffer exposes the difference.

Then it stops being a matter of taste:

const size_t i = col * height + row;  // not row * width + col

For a 611 by 397 image, that line stays within the allocation. It maps the same 242,567 pairs to the same number of addresses, but in a different order.

CUDA reports no fault. The program therefore compares every output pixel with a CPU reference.

The course rule in CUDA-CODE-STYLE.md is: row comes from y, and col comes from x. For another layout, transpose the data or change the index formula and explain it in a comment. Do not hide that choice in the launch setup.

One thread, one pixel, and forty lines of file format

The full program is in code/day07-2d-grids/grayscale.cu. It creates a 611 by 397 test image, writes and reads a binary PPM, converts it on the GPU, checks each pixel, and writes a binary PGM.

One thread owns one pixel, and two implementations flatten the index. The kernel uses a 2D grid, while the CPU reference uses nested loops. Their separate index code makes the comparison useful.

__global__ void grayscale(const unsigned char* rgb, unsigned char* gray,
                          size_t width, size_t height) {
    const size_t col =
        blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
    const size_t row =
        blockIdx.y * static_cast<size_t>(blockDim.y) + threadIdx.y;
    if (row < height && col < width) {
        const size_t i = row * width + col;
        gray[i] = luma(rgb[kChannels * i + 0], rgb[kChannels * i + 1],
                       rgb[kChannels * i + 2]);
    }
}

The guard tests both axes because the grid rounds up twice. ceil(611 / 32) is 20 and ceil(397 / 8) is 50, so the launch uses 640 by 400 threads for 242,567 pixels.

The final block is partial on both axes: 611 is 19 blocks of 32 plus 3, while 397 is 49 blocks of 8 plus 5. Either half of the guard alone leaves one edge unchecked.

    const dim3 block(kBlockDimX, kBlockDimY);
    const dim3 grid(static_cast<unsigned int>(
                        ceilDiv(kWidth, static_cast<size_t>(kBlockDimX))),
                    static_cast<unsigned int>(
                        ceilDiv(kHeight, static_cast<size_t>(kBlockDimY))));

    grayscale<<<grid, block>>>(d_rgb, d_gray, kWidth, kHeight);
    CUDA_CHECK(cudaGetLastError());
    CUDA_CHECK(cudaDeviceSynchronize());

No image library. A binary PGM starts P5, a binary PPM starts P6, then width, height and the maximum sample value as ASCII decimal, one whitespace character, then the bytes (https://netpbm.sourceforge.net/doc/ppm.html and https://netpbm.sourceforge.net/doc/pgm.html , both checked 2026-08-30). The writer is the whole of it. The reader is longer only because a # comment may sit between any two header fields.

static bool writeNetpbm(const char* path, const unsigned char* pixels,
                        size_t width, size_t height, size_t channels) {
    std::FILE* file = std::fopen(path, "wb");
    if (file == nullptr) {
        return false;
    }
    std::fprintf(file, "P%d\n%zu %zu\n255\n", (channels == 1) ? 5 : 6, width,
                 height);
    const size_t n = width * height * channels;
    const bool wrote = std::fwrite(pixels, 1, n, file) == n;
    return (std::fclose(file) == 0) && wrote;
}

Nothing is timed. Day 9 is where timing arrives, for the reasons day 5 gave.

The luma weights live in one __host__ __device__ function (execution space specifiers from day 1). The kernel and CPU reference run the same arithmetic, so this check tests only the indexing.

The weights are integers that sum to 256. This avoids a false mismatch if the GPU fuses a multiply and add but the host compiler does not.

Note. Real image buffers are usually padded. cudaMallocPitch hands back a row stride in bytes, the pitch, which is generally larger than width * sizeof(T) so that every row starts at a friendly address (https://docs.nvidia.com/cuda/cuda-runtime-api/group__CUDART__MEMORY.html , checked 2026-08-30). Indexing then reads row * pitch + col * sizeof(T), in bytes, and cudaMemcpy2D moves such a buffer without the padding. This program uses a plain cudaMalloc so that row * width + col is the only formula on the page.

Results

Re-verified without behavioral drift on a Tesla T4 with driver 580.173.02 and CUDA 13.0 (V13.0.88) on 2026-09-02. The original CUDA 12.6 transcript remains beside the new run in the page's evidence array.

GPU: Tesla T4 (compute capability 7.5)
image 611 x 397, 242567 pixels
block (32, 8), grid (20, 50), 256000 threads
13433 threads have no pixel
wrote day07_input.ppm and day07_gray.pgm
all 242567 pixels match the CPU reference

The 611 x 397 image has 242,567 pixels. A (32, 8) block needs a (20, 50) grid, which launches 256,000 threads. The bounds check rejects 13,433 threads, or 5.2 percent of the launch, because both axes round up.

The bounds check has to test both. A thread can be inside the image horizontally and past the bottom edge, or the reverse, and either one writes outside the buffer.

Run it yourself

Run the program on a CUDA GPU. The build line is in the repo's README:

nvcc -std=c++17 -O3 -arch=sm_75 -o grayscale grayscale.cu
./grayscale

The program needs a writable working directory and about 1 MB for its two output files. This page has no Compiler Explorer example because the project has not checked whether its CUDA sandbox permits file writes; FACT-SHEET.md section 4 tracks that limit.

Exercise

Write solve() for the grayscale kernel: a 2D grid over an RGB image, one thread per pixel, one byte of output each. The harness supplies the buffers and the shapes.

Time: 25 to 40 minutes. Submit: your grayscale.cu, plus one sentence saying which of the seven cases you expected to fail before you ran it.

Check: the harness owns main(), the buffers, and the stream. It runs your function on seven shapes: (0,0), (1,1), (31,33), (611,613), (1024,1024), (1025,1023), and (4096,4096).

It reports the smallest failing case, your launch configuration, and the first ten wrong indices with their block and thread. If (1024,1024) passes but (1025,1023) fails, half of the guard is missing. Full contract at /reference/harness.

Hint 1

Count the threads, not the pixels. How many threads does your launch start for a 611 by 613 image, and what is every one of them told to do? Now answer that for 1024 by 1024.

Hint 2

(1025,1023) rounds up on one axis and down to nothing on the other. If your kernel passes (1024,1024) and fails that one, which comparison was carrying the whole guard, and what was the other one going to do?

Solution

Two index lines and one conjunction:

const size_t col =
    blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
const size_t row =
    blockIdx.y * static_cast<size_t>(blockDim.y) + threadIdx.y;
if (row < height && col < width) {
    const size_t i = row * width + col;
    // ...
}

Not if (i < width * height), which is the one-dimensional guard in disguise. It admits every thread whose flattened index lands inside the buffer, including the ones whose col ran off the right-hand edge into the next row.

Not if (row >= height || col >= width) return; either. The early return behaves the same today and hangs on day 13, when a __syncthreads() appears below it and the returned threads never arrive.

A 2D launch rounds up each axis on its own, so the bounds check needs one comparison per axis. A single test on the flat index cannot show which axis exceeded its limit.

Pitfalls

Your output is transposed and nothing crashed. col * height + row covers the same address range as row * width + col, so no thread leaves the buffer and no check fires. Keep a CPU reference, so a wrong bijection cannot pass. Day 12 is where swapping the two on purpose becomes the point.

You guard row only, and every pixel still matches. A thread with col = width + c writes the address for pixel (row + 1, c), which another thread also writes. The invalid writes start on the last row, where no next row exists.

For a 611 by 397 image with a 32-wide block, 396 * 611 + col leaves the allocation at col = 611. The next 29 threads write past the allocation.

compute-sanitizer --tool memcheck ./grayscale can find them, but the pixel check cannot. Compute Sanitizer is day 61; this page quotes no sanitizer run.

The launch fails with invalid configuration argument. A 2D or 3D block can exceed 1,024 threads even when each axis looks small: dim3 block(32, 32, 2) has 2,048 threads.

The grid limits also differ by axis. gridDim.x reaches 2^31 - 1, while gridDim.y and gridDim.z stop at 65,535 (Table 30, https://docs.nvidia.com/cuda/cuda-programming-guide/05-appendices/compute-capabilities.html , checked 2026-08-30).

cudaGetLastError() reports either error at the launch. See /errors/invalid-configuration-argument.

You size the grid with width / blockDim.x. Integer division truncates: 611 over 32 is 19, so the right-hand three columns get no thread. The harness calls those never written rather than wrong, which is the faster diagnosis. The ceiling is (n + b - 1) / b, kept here in one constexpr function so it can size the grid and appear in a static_assert.

Go deeper

Next

Day 8 takes the guard you just wrote and asks what happens when one grid has to cover any size at all, which is where the rounded-up edge becomes a loop condition. Day 10 sweeps the block shape across this same image kernel, so the 32 in dim3 block(32, 8) stops being something you were told. Both sit in module 1.