code/day99-forward-pass/forward_pass.cuThis is the source used by the lesson and its recorded evidence. Compile commands and expected output live in the directory README.
// SPDX-License-Identifier: MIT
// Day 99: one forward and backward pass.
// H = ReLU(X W^T + b), followed by stable softmax cross-entropy.
// Backward computes dZ, dW = dZ^T X, db, and dX = dZ W.
// Build: nvcc -std=c++17 -O3 -arch=sm_75 -o forward_pass forward_pass.cu
// VERIFIED: Tesla T4, driver 580.173.02, CUDA 12.6, 2026-09-02.
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <vector>
#include <cuda_runtime.h>
#define CUDA_CHECK(call) \
do { \
cudaError_t err_ = (call); \
if (err_ != cudaSuccess) { \
std::fprintf(stderr, "CUDA error %s:%d: %s: %s\n", __FILE__, \
__LINE__, #call, cudaGetErrorString(err_)); \
std::exit(EXIT_FAILURE); \
} \
} while (0)
constexpr int kBatch = 32;
constexpr int kInputs = 16;
constexpr int kClasses = 8;
constexpr int kSoftmaxThreads = 32;
constexpr double kAtol = 2e-5;
constexpr double kRtol = 2e-4;
// X is [batch, inputs]. W is [classes, inputs], so forward reads W
// transposed. Z is retained because backward needs its sign for ReLU.
__global__ void linearReluForward(const float* __restrict__ x,
const float* __restrict__ w,
const float* __restrict__ bias,
float* __restrict__ z,
float* __restrict__ h) {
const int c = blockIdx.x * blockDim.x + threadIdx.x;
const int row = blockIdx.y * blockDim.y + threadIdx.y;
if (row >= kBatch || c >= kClasses) return;
float sum = bias[c];
for (int i = 0; i < kInputs; ++i) {
sum += x[row * kInputs + i] * w[c * kInputs + i];
}
z[row * kClasses + c] = sum;
h[row * kClasses + c] = fmaxf(sum, 0.0f);
}
// One block owns one row. Subtracting rowMax before expf is the stability
// step. The probabilities also give dH = (p - one_hot(label)) / batch.
__global__ void stableSoftmaxCrossEntropy(const float* __restrict__ h,
const int* __restrict__ labels,
float* __restrict__ rowLoss,
float* __restrict__ dH) {
__shared__ float scratch[kSoftmaxThreads];
const int row = blockIdx.x;
const int lane = threadIdx.x;
float value = lane < kClasses ? h[row * kClasses + lane] : -INFINITY;
scratch[lane] = value;
__syncthreads();
for (int stride = kSoftmaxThreads / 2; stride > 0; stride /= 2) {
if (lane < stride) scratch[lane] = fmaxf(scratch[lane], scratch[lane + stride]);
__syncthreads();
}
const float rowMax = scratch[0];
const float numerator = lane < kClasses ? expf(value - rowMax) : 0.0f;
scratch[lane] = numerator;
__syncthreads();
for (int stride = kSoftmaxThreads / 2; stride > 0; stride /= 2) {
if (lane < stride) scratch[lane] += scratch[lane + stride];
__syncthreads();
}
const float denominator = scratch[0];
if (lane < kClasses) {
const float target = lane == labels[row] ? 1.0f : 0.0f;
dH[row * kClasses + lane] =
(numerator / denominator - target) / static_cast<float>(kBatch);
}
if (lane == 0) {
const float labelled = h[row * kClasses + labels[row]];
rowLoss[row] = logf(denominator) - (labelled - rowMax);
}
}
// dZ = dH * 1[Z > 0].
__global__ void reluBackward(const float* __restrict__ z,
const float* __restrict__ dH,
float* __restrict__ dZ) {
const int idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < kBatch * kClasses) dZ[idx] = z[idx] > 0.0f ? dH[idx] : 0.0f;
}
// dW = dZ^T X. This is the transposed GEMM shape training adds.
// snippet: transposed-gradient
__global__ void weightGradient(const float* __restrict__ dZ,
const float* __restrict__ x,
float* __restrict__ dW) {
const int i = blockIdx.x * blockDim.x + threadIdx.x;
const int c = blockIdx.y * blockDim.y + threadIdx.y;
if (i >= kInputs || c >= kClasses) return;
float sum = 0.0f;
for (int row = 0; row < kBatch; ++row) {
sum += dZ[row * kClasses + c] * x[row * kInputs + i];
}
dW[c * kInputs + i] = sum;
}
// end snippet
// db[c] is a column reduction over dZ's batch dimension.
__global__ void biasGradient(const float* __restrict__ dZ,
float* __restrict__ db) {
const int c = blockIdx.x * blockDim.x + threadIdx.x;
if (c >= kClasses) return;
float sum = 0.0f;
for (int row = 0; row < kBatch; ++row) sum += dZ[row * kClasses + c];
db[c] = sum;
}
// dX = dZ W. Unlike forward, this reads W in its stored orientation.
__global__ void inputGradient(const float* __restrict__ dZ,
const float* __restrict__ w,
float* __restrict__ dX) {
const int i = blockIdx.x * blockDim.x + threadIdx.x;
const int row = blockIdx.y * blockDim.y + threadIdx.y;
if (row >= kBatch || i >= kInputs) return;
float sum = 0.0f;
for (int c = 0; c < kClasses; ++c) {
sum += dZ[row * kClasses + c] * w[c * kInputs + i];
}
dX[row * kInputs + i] = sum;
}
static void cpuReference(const std::vector<float>& x,
const std::vector<float>& w,
const std::vector<float>& bias,
const std::vector<int>& labels,
std::vector<double>* z, std::vector<double>* h,
std::vector<double>* rowLoss,
std::vector<double>* dH, std::vector<double>* dZ,
std::vector<double>* dW, std::vector<double>* db,
std::vector<double>* dX) {
for (int row = 0; row < kBatch; ++row) {
for (int c = 0; c < kClasses; ++c) {
double sum = bias[c];
for (int i = 0; i < kInputs; ++i) {
sum += static_cast<double>(x[row * kInputs + i]) *
w[c * kInputs + i];
}
(*z)[row * kClasses + c] = sum;
(*h)[row * kClasses + c] = std::fmax(sum, 0.0);
}
double rowMax = (*h)[row * kClasses];
for (int c = 1; c < kClasses; ++c) {
rowMax = std::fmax(rowMax, (*h)[row * kClasses + c]);
}
double denominator = 0.0;
for (int c = 0; c < kClasses; ++c) {
denominator += std::exp((*h)[row * kClasses + c] - rowMax);
}
(*rowLoss)[row] = std::log(denominator) -
((*h)[row * kClasses + labels[row]] - rowMax);
for (int c = 0; c < kClasses; ++c) {
const double p = std::exp((*h)[row * kClasses + c] - rowMax) /
denominator;
const int idx = row * kClasses + c;
(*dH)[idx] = (p - (c == labels[row] ? 1.0 : 0.0)) / kBatch;
(*dZ)[idx] = (*z)[idx] > 0.0 ? (*dH)[idx] : 0.0;
}
}
for (int c = 0; c < kClasses; ++c) {
for (int i = 0; i < kInputs; ++i) {
double sum = 0.0;
for (int row = 0; row < kBatch; ++row) {
sum += (*dZ)[row * kClasses + c] * x[row * kInputs + i];
}
(*dW)[c * kInputs + i] = sum;
}
double sum = 0.0;
for (int row = 0; row < kBatch; ++row) sum += (*dZ)[row * kClasses + c];
(*db)[c] = sum;
}
for (int row = 0; row < kBatch; ++row) {
for (int i = 0; i < kInputs; ++i) {
double sum = 0.0;
for (int c = 0; c < kClasses; ++c) {
sum += (*dZ)[row * kClasses + c] * w[c * kInputs + i];
}
(*dX)[row * kInputs + i] = sum;
}
}
}
static bool checkTensor(const char* name, const std::vector<float>& got,
const std::vector<double>& want) {
double maxError = 0.0;
size_t firstBad = got.size();
for (size_t i = 0; i < got.size(); ++i) {
const double error = std::fabs(static_cast<double>(got[i]) - want[i]);
maxError = std::fmax(maxError, error);
const double bound = kAtol + kRtol * std::fabs(want[i]);
if ((!std::isfinite(got[i]) || error > bound) && firstBad == got.size()) {
firstBad = i;
}
}
std::printf(" %-8s max abs error %.3e %s\n", name, maxError,
firstBad == got.size() ? "PASS" : "FAIL");
if (firstBad == got.size()) return true;
std::fprintf(stderr, "%s[%zu] got %.9g, want %.12g\n", name, firstBad,
got[firstBad], want[firstBad]);
return false;
}
int main() {
cudaDeviceProp prop;
int driverVersion = 0;
int runtimeVersion = 0;
CUDA_CHECK(cudaGetDeviceProperties(&prop, 0));
CUDA_CHECK(cudaDriverGetVersion(&driverVersion));
CUDA_CHECK(cudaRuntimeGetVersion(&runtimeVersion));
std::printf("GPU: %s (compute capability %d.%d)\n", prop.name, prop.major,
prop.minor);
std::printf("CUDA driver API %d, runtime %d\n", driverVersion,
runtimeVersion);
const size_t nAct = static_cast<size_t>(kBatch) * kClasses;
const size_t nX = static_cast<size_t>(kBatch) * kInputs;
const size_t nW = static_cast<size_t>(kClasses) * kInputs;
std::vector<float> x(nX), w(nW), bias(kClasses);
std::vector<int> labels(kBatch);
for (size_t i = 0; i < nX; ++i)
x[i] = static_cast<float>(static_cast<int>((i * 13 + 5) % 23) - 11) / 11.0f;
for (size_t i = 0; i < nW; ++i)
w[i] = static_cast<float>(static_cast<int>((i * 7 + 3) % 19) - 9) / 12.0f;
for (int c = 0; c < kClasses; ++c) bias[c] = (c - 3) * 0.075f;
for (int row = 0; row < kBatch; ++row) labels[row] = (row * 5 + 1) % kClasses;
std::vector<double> refZ(nAct), refH(nAct), refLoss(kBatch), refDH(nAct);
std::vector<double> refDZ(nAct), refDW(nW), refDb(kClasses), refDX(nX);
cpuReference(x, w, bias, labels, &refZ, &refH, &refLoss, &refDH, &refDZ,
&refDW, &refDb, &refDX);
float *dX, *dW, *dBias, *dZ, *dH, *dLoss, *dDH, *dDZ, *dDW, *dDb, *dDX;
int* dLabels;
CUDA_CHECK(cudaMalloc(&dX, nX * sizeof(float)));
CUDA_CHECK(cudaMalloc(&dW, nW * sizeof(float)));
CUDA_CHECK(cudaMalloc(&dBias, kClasses * sizeof(float)));
CUDA_CHECK(cudaMalloc(&dLabels, kBatch * sizeof(int)));
CUDA_CHECK(cudaMalloc(&dZ, nAct * sizeof(float)));
CUDA_CHECK(cudaMalloc(&dH, nAct * sizeof(float)));
CUDA_CHECK(cudaMalloc(&dLoss, kBatch * sizeof(float)));
CUDA_CHECK(cudaMalloc(&dDH, nAct * sizeof(float)));
CUDA_CHECK(cudaMalloc(&dDZ, nAct * sizeof(float)));
CUDA_CHECK(cudaMalloc(&dDW, nW * sizeof(float)));
CUDA_CHECK(cudaMalloc(&dDb, kClasses * sizeof(float)));
CUDA_CHECK(cudaMalloc(&dDX, nX * sizeof(float)));
CUDA_CHECK(cudaMemcpy(dX, x.data(), nX * sizeof(float), cudaMemcpyHostToDevice));
CUDA_CHECK(cudaMemcpy(dW, w.data(), nW * sizeof(float), cudaMemcpyHostToDevice));
CUDA_CHECK(cudaMemcpy(dBias, bias.data(), kClasses * sizeof(float), cudaMemcpyHostToDevice));
CUDA_CHECK(cudaMemcpy(dLabels, labels.data(), kBatch * sizeof(int), cudaMemcpyHostToDevice));
const dim3 tile(16, 16);
// snippet: pass-launches
linearReluForward<<<dim3(1, 2), tile>>>(dX, dW, dBias, dZ, dH);
stableSoftmaxCrossEntropy<<<kBatch, kSoftmaxThreads>>>(dH, dLabels, dLoss, dDH);
reluBackward<<<(nAct + 255) / 256, 256>>>(dZ, dDH, dDZ);
weightGradient<<<dim3(1, 1), tile>>>(dDZ, dX, dDW);
biasGradient<<<1, 32>>>(dDZ, dDb);
inputGradient<<<dim3(1, 2), tile>>>(dDZ, dW, dDX);
// end snippet
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
std::vector<float> gotZ(nAct), gotH(nAct), gotLoss(kBatch), gotDH(nAct);
std::vector<float> gotDZ(nAct), gotDW(nW), gotDb(kClasses), gotDX(nX);
#define COPY_BACK(host, device) \
CUDA_CHECK(cudaMemcpy((host).data(), device, (host).size() * sizeof(float), \
cudaMemcpyDeviceToHost))
COPY_BACK(gotZ, dZ); COPY_BACK(gotH, dH); COPY_BACK(gotLoss, dLoss);
COPY_BACK(gotDH, dDH); COPY_BACK(gotDZ, dDZ); COPY_BACK(gotDW, dDW);
COPY_BACK(gotDb, dDb); COPY_BACK(gotDX, dDX);
#undef COPY_BACK
std::printf("one-layer pass: batch %d, inputs %d, classes %d\n", kBatch,
kInputs, kClasses);
std::printf("forward gates\n");
bool ok = checkTensor("Z", gotZ, refZ);
ok = checkTensor("H", gotH, refH) && ok;
ok = checkTensor("loss", gotLoss, refLoss) && ok;
ok = checkTensor("dH", gotDH, refDH) && ok;
std::printf("backward gradient gates\n");
ok = checkTensor("dZ", gotDZ, refDZ) && ok;
ok = checkTensor("dW", gotDW, refDW) && ok;
ok = checkTensor("db", gotDb, refDb) && ok;
ok = checkTensor("dX", gotDX, refDX) && ok;
double meanLoss = 0.0;
for (float value : gotLoss) meanLoss += value;
std::printf("mean cross-entropy %.6f\n", meanLoss / kBatch);
CUDA_CHECK(cudaFree(dX)); CUDA_CHECK(cudaFree(dW));
CUDA_CHECK(cudaFree(dBias)); CUDA_CHECK(cudaFree(dLabels));
CUDA_CHECK(cudaFree(dZ)); CUDA_CHECK(cudaFree(dH));
CUDA_CHECK(cudaFree(dLoss)); CUDA_CHECK(cudaFree(dDH));
CUDA_CHECK(cudaFree(dDZ)); CUDA_CHECK(cudaFree(dDW));
CUDA_CHECK(cudaFree(dDb)); CUDA_CHECK(cudaFree(dDX));
std::printf("%s\n", ok ? "PASS" : "FAIL");
return ok ? EXIT_SUCCESS : EXIT_FAILURE;
}