code/day57-graph-update/graph_update.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 57: graph updates and conditional nodes.
//
// What it does: solves one tridiagonal system with weighted Jacobi sweeps,
// three times, changing only who decides when to stop. Strategy 1 is the
// host loop with a blocking residual copy per iteration. Strategy 2 launches
// one four-kernel graph per iteration and still copies the residual back.
// Strategy 3 is a conditional WHILE graph (CUDA 12.4+): one launch, and the
// convergence test runs on the device until the residual drops below the
// bound. The program then changes the relaxation weight baked into the
// captured graph with cudaGraphExecUpdate instead of re-instantiating.
//
// The three strategies must agree: same iteration count, bit-identical final
// state, and a true |Ax - b| residual checked in double on the host. Any
// disagreement is a failure, not a footnote.
//
// On a toolkit older than 12.4 the conditional path compiles out and the
// program says so; strategies 1 and 2 still run, which is the fallback the
// lesson names for that toolkit.
//
// Build: nvcc -std=c++17 -O3 -arch=sm_75 -lineinfo -o graph_update \
// graph_update.cu
// Run: ./graph_update
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <vector>
#include <cuda_runtime.h>
#include <nvtx3/nvToolsExt.h>
// The one error macro. This file is standalone, so it carries its own
// verbatim copy. `err_` has a trailing underscore so it cannot collide with a
// variable at the call site, and the do/while makes the macro one statement
// so it survives a braceless `if`.
#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)
// The system is A x = b with A = tridiag(-1, 2.5, -1). The Jacobi iteration
// matrix for it has infinity-norm 2/2.5 = 0.8, so the sweep-to-sweep change
// shrinks by at least 0.8 per sweep and the loop below cannot run forever.
// kMaxIterations is the cap that makes that a guarantee rather than a hope:
// a device-side loop with no cap has no ctrl-C.
//
// 2^22 floats is 16 MiB per buffer, past this card's 4 MiB L2, so each sweep
// is a real trip through DRAM and one iteration is hundreds of microseconds,
// wide enough to read on a timeline.
constexpr size_t kElems = size_t{1} << 22;
constexpr int kThreadsPerBlock = 256; // 8 warps
constexpr int kReduceBlocks = 256;
constexpr int kWarmupRuns = 3;
constexpr int kTimedRuns = 10;
constexpr float kDiag = 2.5f;
constexpr float kInvDiag = 1.0f / kDiag;
constexpr float kTolerance = 1e-5f;
constexpr int kMaxIterations = 500;
constexpr float kOmegaFirst = 1.0f;
constexpr float kOmegaSecond = 0.7f;
// Measured convergence is |x_k - x_{k-1}|, not |A x - b|. For this system
// the two differ by at most kDiag times one more contraction step, so a
// converged state must satisfy |A x - b| <= kDiag * 0.8 * kTolerance in
// exact arithmetic. The factor 10 is slack for float rounding and for the
// omega = 0.7 rerun, whose map contracts at 0.86 instead of 0.8.
constexpr double kResidualBound = 10.0 * kDiag * kTolerance;
static_assert(kThreadsPerBlock % 32 == 0,
"block size must be a whole number of warps");
static_assert((kReduceBlocks & (kReduceBlocks - 1)) == 0,
"the final reduce halves its stride, so the partial count must "
"be a power of two");
static_assert(kReduceBlocks <= 1024,
"the final reduce runs as one block of kReduceBlocks threads");
// One weighted Jacobi sweep: xNext[i] = (1-omega)*x[i] + omega*(b[i] +
// x[i-1] + x[i+1]) / 2.5, with a zero boundary on both ends. One thread owns
// one element.
//
// Memory: consecutive lanes read consecutive elements of x three times at
// offsets -1, 0 and +1, three overlapping 128-byte spans, all coalesced; b
// and xNext are one contiguous span each.
//
// Launch assumption: gridDim.x * blockDim.x >= n. omega arrives by value,
// which matters on this day: a value captured into a graph is frozen in the
// instantiated executable until cudaGraphExecUpdate replaces it.
__global__ void relaxJacobi(const float* __restrict__ x,
const float* __restrict__ b,
float* __restrict__ xNext, float omega, size_t n) {
const size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
if (i < n) {
const float left = (i > 0) ? x[i - 1] : 0.0f;
const float right = (i + 1 < n) ? x[i + 1] : 0.0f;
const float jacobi = (b[i] + left + right) * kInvDiag;
xNext[i] = (1.0f - omega) * x[i] + omega * jacobi;
}
}
// Per-block maximum of |xNew[i] - xOld[i]| into partials[blockIdx.x]. Each
// thread walks a grid-stride loop, so the launch is always kReduceBlocks
// blocks whatever n is.
//
// Memory: consecutive lanes read consecutive elements of both inputs,
// coalesced; each block writes one float.
//
// Launch assumption: exactly kThreadsPerBlock threads per block, because the
// shared array and the halving loop are sized from it. Every thread reaches
// the barrier; the guard covers the loads, not the barrier.
__global__ void residualPartial(const float* __restrict__ xNew,
const float* __restrict__ xOld,
float* __restrict__ partials, size_t n) {
__shared__ float tile[kThreadsPerBlock];
const unsigned int tid = threadIdx.x;
const size_t step = gridDim.x * static_cast<size_t>(blockDim.x);
float best = 0.0f;
for (size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + tid; i < n;
i += step) {
best = fmaxf(best, fabsf(xNew[i] - xOld[i]));
}
tile[tid] = best;
__syncthreads();
for (unsigned int half = kThreadsPerBlock / 2; half > 0; half /= 2) {
if (tid < half) {
tile[tid] = fmaxf(tile[tid], tile[tid + half]);
}
__syncthreads();
}
if (tid == 0) {
partials[blockIdx.x] = tile[0];
}
}
// Reduces the kReduceBlocks partial maxima to one residual. This is the
// finish the host-decided strategies use; the host then copies residual[0]
// back and compares, which is the round trip this whole day is about.
//
// Memory: one coalesced read of kReduceBlocks floats, one write.
//
// Launch assumption: one block of exactly kReduceBlocks threads.
__global__ void residualFinish(const float* __restrict__ partials,
float* __restrict__ residual) {
__shared__ float tile[kReduceBlocks];
const unsigned int tid = threadIdx.x;
tile[tid] = partials[tid];
__syncthreads();
for (unsigned int half = kReduceBlocks / 2; half > 0; half /= 2) {
if (tid < half) {
tile[tid] = fmaxf(tile[tid], tile[tid + half]);
}
__syncthreads();
}
if (tid == 0) {
residual[0] = tile[0];
}
}
#if CUDART_VERSION >= 12040
// The same finish, plus the decision. Thread 0 bumps the iteration counter
// and sets the conditional handle: 1 keeps the WHILE body running, 0 stops
// it. The reduce is line for line the one in residualFinish, so the value
// the device compares against kTolerance is bit-identical to the value the
// host-decided strategies copied back, and the three iteration counts can be
// required to match exactly.
//
// Memory: one coalesced read of kReduceBlocks floats, two scalar writes.
//
// Launch assumption: one block of exactly kReduceBlocks threads.
__global__ void decideAndCount(const float* __restrict__ partials,
float* __restrict__ residual,
int* __restrict__ iterations, float tol,
int maxIterations,
cudaGraphConditionalHandle handle) {
__shared__ float tile[kReduceBlocks];
const unsigned int tid = threadIdx.x;
tile[tid] = partials[tid];
__syncthreads();
for (unsigned int half = kReduceBlocks / 2; half > 0; half /= 2) {
if (tid < half) {
tile[tid] = fmaxf(tile[tid], tile[tid + half]);
}
__syncthreads();
}
// snippet: device-decision
if (tid == 0) {
residual[0] = tile[0];
const int done = *iterations + 1;
*iterations = done;
const unsigned int keepGoing =
(tile[0] > tol && done < maxIterations) ? 1u : 0u;
cudaGraphSetConditional(handle, keepGoing);
}
// end snippet
}
#endif // CUDART_VERSION >= 12040
// Times a launch with CUDA events and returns the mean milliseconds per run.
//
// This is the course's one timing helper, copied verbatim. Here each
// "launch" is a whole solve, including its resets and, for the host-decided
// strategies, every blocking copy inside it. The gaps the host causes sit
// between the two events, which is exactly what this day wants measured.
template <typename LaunchFn>
static float timeKernel(LaunchFn launch) {
cudaEvent_t start, stop;
CUDA_CHECK(cudaEventCreate(&start));
CUDA_CHECK(cudaEventCreate(&stop));
// Warm up this kernel, not just the first kernel in the program. Lazy
// module loading has been the default since CUDA 12.2 on Linux, so the
// first launch of each kernel pays its own load.
for (int i = 0; i < kWarmupRuns; ++i) {
launch();
}
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaEventRecord(start));
for (int i = 0; i < kTimedRuns; ++i) {
launch();
}
CUDA_CHECK(cudaEventRecord(stop));
CUDA_CHECK(cudaEventSynchronize(stop));
CUDA_CHECK(cudaGetLastError());
float ms = 0.0f;
CUDA_CHECK(cudaEventElapsedTime(&ms, start, stop));
CUDA_CHECK(cudaEventDestroy(start));
CUDA_CHECK(cudaEventDestroy(stop));
return ms / kTimedRuns;
}
// The state every solve starts from: x = 0 and the iteration counter at
// zero. Both memsets are stream-ordered on the default stream, ahead of
// whatever the solve enqueues next.
static void resetState(float* d_x, int* d_iterations) {
CUDA_CHECK(cudaMemset(d_x, 0, kElems * sizeof(float)));
CUDA_CHECK(cudaMemset(d_iterations, 0, sizeof(int)));
}
// One body iteration: two sweeps (x -> y -> x, so the buffer roles never
// change, which is what lets a fixed graph re-run it), then the residual of
// the last two states. Shared by strategy 1 (stream 0) and the captures.
static void enqueueSweeps(float* d_x, float* d_y, const float* d_b,
float* d_partials, float omega, cudaStream_t stream) {
const int blocks =
static_cast<int>((kElems + kThreadsPerBlock - 1) / kThreadsPerBlock);
relaxJacobi<<<blocks, kThreadsPerBlock, 0, stream>>>(d_x, d_b, d_y, omega,
kElems);
relaxJacobi<<<blocks, kThreadsPerBlock, 0, stream>>>(d_y, d_b, d_x, omega,
kElems);
residualPartial<<<kReduceBlocks, kThreadsPerBlock, 0, stream>>>(
d_x, d_y, d_partials, kElems);
}
// Strategy 1: the host owns the loop. Four launches, then a blocking 4-byte
// copy that cannot return until everything queued before it has finished.
// That copy is the per-iteration synchronisation this day prices; day 41
// measured its shape on a timeline. The copy is also the execution-error
// check for the iteration, so the only added check per pass is
// cudaGetLastError for a refused launch.
static int solveHostLoop(float* d_x, float* d_y, const float* d_b,
float* d_partials, float* d_residual, float omega,
float* residualOut) {
int iterations = 0;
float r = 0.0f;
do {
enqueueSweeps(d_x, d_y, d_b, d_partials, omega, nullptr);
residualFinish<<<1, kReduceBlocks>>>(d_partials, d_residual);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(
cudaMemcpy(&r, d_residual, sizeof(float), cudaMemcpyDeviceToHost));
++iterations;
} while (r > kTolerance && iterations < kMaxIterations);
*residualOut = r;
return iterations;
}
// Strategy 2: the same loop with its four launches folded into one graph
// launch. The round trip per iteration is still there; only the launch
// overhead went away. This is the fallback path on a toolkit older than
// 12.4, and it is also the executable that cudaGraphExecUpdate patches
// later.
static int solveGraphLoop(cudaGraphExec_t exec, float* d_residual,
float* residualOut) {
int iterations = 0;
float r = 0.0f;
do {
CUDA_CHECK(cudaGraphLaunch(exec, nullptr));
CUDA_CHECK(
cudaMemcpy(&r, d_residual, sizeof(float), cudaMemcpyDeviceToHost));
++iterations;
} while (r > kTolerance && iterations < kMaxIterations);
*residualOut = r;
return iterations;
}
// Captures one body iteration as a graph: two sweeps, the partial residual,
// and the finish that leaves the residual where the host can copy it from.
// The capture happens on a created stream because the legacy default stream
// cannot be captured; the executable still launches on stream 0 later.
//
// No cudaGetLastError inside the capture: nothing runs during capture, and
// a bad launch surfaces as an error on cudaStreamEndCapture instead.
static cudaGraph_t captureIterationGraph(float* d_x, float* d_y,
const float* d_b, float* d_partials,
float* d_residual, float omega) {
cudaStream_t stream;
CUDA_CHECK(cudaStreamCreate(&stream));
CUDA_CHECK(cudaStreamBeginCapture(stream, cudaStreamCaptureModeGlobal));
enqueueSweeps(d_x, d_y, d_b, d_partials, omega, stream);
residualFinish<<<1, kReduceBlocks, 0, stream>>>(d_partials, d_residual);
cudaGraph_t graph = nullptr;
CUDA_CHECK(cudaStreamEndCapture(stream, &graph));
CUDA_CHECK(cudaStreamDestroy(stream));
return graph;
}
#if CUDART_VERSION >= 12040
// Strategy 3's graph: a conditional WHILE node whose body is one iteration,
// with decideAndCount in place of residualFinish. The handle defaults to 1
// on every graph launch (cudaGraphCondAssignDefault), so the body always
// runs at least once and the shape is a do-while, matching the two host
// loops above. Conditional nodes and cudaGraphSetConditional are CUDA 12.4+.
static cudaGraph_t buildConditionalGraph(float* d_x, float* d_y,
const float* d_b, float* d_partials,
float* d_residual, int* d_iterations,
float omega) {
cudaGraph_t graph = nullptr;
CUDA_CHECK(cudaGraphCreate(&graph, 0));
// snippet: while-node
cudaGraphConditionalHandle handle;
CUDA_CHECK(cudaGraphConditionalHandleCreate(&handle, graph, 1,
cudaGraphCondAssignDefault));
cudaGraphNodeParams params = {};
params.type = cudaGraphNodeTypeConditional;
params.conditional.handle = handle;
params.conditional.type = cudaGraphCondTypeWhile;
params.conditional.size = 1;
cudaGraphNode_t node = nullptr;
#if CUDART_VERSION >= 13000
CUDA_CHECK(cudaGraphAddNode(&node, graph, nullptr, nullptr, 0, ¶ms));
#else
CUDA_CHECK(cudaGraphAddNode(&node, graph, nullptr, 0, ¶ms));
#endif
// The body graph is owned by the conditional node; capturing into it is
// how kernels get inside. It is never instantiated or destroyed here.
cudaGraph_t body = params.conditional.phGraph_out[0];
cudaStream_t stream;
CUDA_CHECK(cudaStreamCreate(&stream));
CUDA_CHECK(cudaStreamBeginCaptureToGraph(stream, body, nullptr, nullptr, 0,
cudaStreamCaptureModeRelaxed));
enqueueSweeps(d_x, d_y, d_b, d_partials, omega, stream);
decideAndCount<<<1, kReduceBlocks, 0, stream>>>(d_partials, d_residual,
d_iterations, kTolerance,
kMaxIterations, handle);
CUDA_CHECK(cudaStreamEndCapture(stream, nullptr));
CUDA_CHECK(cudaStreamDestroy(stream));
// end snippet
return graph;
}
#endif // CUDART_VERSION >= 12040
// The true residual |A x - b| in double, on the host. This is the gate that
// says the strategies did not merely agree with each other, they solved the
// system: agreeing on a wrong answer three times would pass every other
// check in this file.
static double trueResidual(const float* x, const float* b, size_t n) {
double worst = 0.0;
for (size_t i = 0; i < n; ++i) {
const double left = (i > 0) ? static_cast<double>(x[i - 1]) : 0.0;
const double right = (i + 1 < n) ? static_cast<double>(x[i + 1]) : 0.0;
const double r = static_cast<double>(b[i]) + left + right -
static_cast<double>(kDiag) * x[i];
worst = std::fmax(worst, std::fabs(r));
}
return worst;
}
int main() {
const int device = 0;
CUDA_CHECK(cudaSetDevice(device));
cudaDeviceProp prop;
CUDA_CHECK(cudaGetDeviceProperties(&prop, device));
std::printf("GPU: %s (compute capability %d.%d)\n", prop.name, prop.major,
prop.minor);
std::printf(
"n = %zu (%.1f MiB per buffer), diag %.2f, tol %.0e, cap %d "
"iterations\n",
kElems, kElems * sizeof(float) / (1024.0 * 1024.0),
static_cast<double>(kDiag), static_cast<double>(kTolerance),
kMaxIterations);
#if CUDART_VERSION >= 12040
std::printf("conditional nodes: available (CUDART %d)\n\n", CUDART_VERSION);
#else
std::printf(
"conditional nodes: need CUDA 12.4+, this is CUDART %d. Strategy 3 "
"is compiled out; strategies 1 and 2 still run.\n\n",
CUDART_VERSION);
#endif
nvtxRangePushA("setup");
const size_t bytes = kElems * sizeof(float);
std::vector<float> h_b(kElems);
for (size_t i = 0; i < kElems; ++i) {
// Deterministic, in [1.0, 1.5]. The exact values are irrelevant; the
// solve must reach the same x from them by all three routes.
h_b[i] = 1.0f + 0.5f * static_cast<float>(i % 97) / 96.0f;
}
float* d_x = nullptr;
float* d_y = nullptr;
float* d_b = nullptr;
float* d_partials = nullptr;
float* d_residual = nullptr;
int* d_iterations = nullptr;
CUDA_CHECK(cudaMalloc(&d_x, bytes));
CUDA_CHECK(cudaMalloc(&d_y, bytes));
CUDA_CHECK(cudaMalloc(&d_b, bytes));
CUDA_CHECK(cudaMalloc(&d_partials, kReduceBlocks * sizeof(float)));
CUDA_CHECK(cudaMalloc(&d_residual, sizeof(float)));
CUDA_CHECK(cudaMalloc(&d_iterations, sizeof(int)));
CUDA_CHECK(cudaMemcpy(d_b, h_b.data(), bytes, cudaMemcpyHostToDevice));
nvtxRangePop();
bool ok = true;
// Correctness pass, one run per strategy, before anything is timed.
std::printf("Correctness pass at omega = %.2f\n",
static_cast<double>(kOmegaFirst));
nvtxRangePushA("host-loop-sync");
resetState(d_x, d_iterations);
float r1 = 0.0f;
const int iter1 =
solveHostLoop(d_x, d_y, d_b, d_partials, d_residual, kOmegaFirst, &r1);
nvtxRangePop();
std::vector<float> h_x1(kElems);
CUDA_CHECK(cudaMemcpy(h_x1.data(), d_x, bytes, cudaMemcpyDeviceToHost));
std::printf(" host-loop-sync %4d iterations, residual %.3e\n", iter1,
static_cast<double>(r1));
nvtxRangePushA("capture");
cudaGraph_t iterationGraph = captureIterationGraph(
d_x, d_y, d_b, d_partials, d_residual, kOmegaFirst);
nvtxRangePop();
nvtxRangePushA("instantiate");
cudaGraphExec_t iterationExec = nullptr;
CUDA_CHECK(cudaGraphInstantiate(&iterationExec, iterationGraph, 0));
nvtxRangePop();
nvtxRangePushA("graph-per-iteration");
resetState(d_x, d_iterations);
float r2 = 0.0f;
const int iter2 = solveGraphLoop(iterationExec, d_residual, &r2);
nvtxRangePop();
std::vector<float> h_x2(kElems);
CUDA_CHECK(cudaMemcpy(h_x2.data(), d_x, bytes, cudaMemcpyDeviceToHost));
std::printf(" graph-per-iteration %4d iterations, residual %.3e\n", iter2,
static_cast<double>(r2));
#if CUDART_VERSION >= 12040
nvtxRangePushA("build-conditional");
cudaGraph_t condGraph = buildConditionalGraph(
d_x, d_y, d_b, d_partials, d_residual, d_iterations, kOmegaFirst);
nvtxRangePop();
nvtxRangePushA("instantiate-conditional");
cudaGraphExec_t condExec = nullptr;
CUDA_CHECK(cudaGraphInstantiate(&condExec, condGraph, 0));
nvtxRangePop();
nvtxRangePushA("conditional-graph");
resetState(d_x, d_iterations);
CUDA_CHECK(cudaGraphLaunch(condExec, nullptr));
CUDA_CHECK(cudaDeviceSynchronize());
nvtxRangePop();
float r3 = 0.0f;
int iter3 = 0;
CUDA_CHECK(
cudaMemcpy(&r3, d_residual, sizeof(float), cudaMemcpyDeviceToHost));
CUDA_CHECK(
cudaMemcpy(&iter3, d_iterations, sizeof(int), cudaMemcpyDeviceToHost));
std::vector<float> h_x3(kElems);
CUDA_CHECK(cudaMemcpy(h_x3.data(), d_x, bytes, cudaMemcpyDeviceToHost));
std::printf(" conditional-graph %4d iterations, residual %.3e\n", iter3,
static_cast<double>(r3));
if (iter1 != iter2 || iter1 != iter3) {
std::fprintf(stderr,
"iteration counts disagree: host %d, graph %d, "
"conditional %d\n",
iter1, iter2, iter3);
ok = false;
}
if (std::memcmp(h_x1.data(), h_x3.data(), bytes) != 0) {
std::fprintf(stderr,
"conditional-graph state differs from host-loop-sync\n");
ok = false;
}
#else
if (iter1 != iter2) {
std::fprintf(stderr, "iteration counts disagree: host %d, graph %d\n",
iter1, iter2);
ok = false;
}
#endif
if (std::memcmp(h_x1.data(), h_x2.data(), bytes) != 0) {
std::fprintf(stderr,
"graph-per-iteration state differs from "
"host-loop-sync\n");
ok = false;
}
if (iter1 >= kMaxIterations || r1 > kTolerance) {
std::fprintf(stderr,
"host-loop-sync did not converge: %d "
"iterations, residual %.3e\n",
iter1, static_cast<double>(r1));
ok = false;
}
const double axb = trueResidual(h_x1.data(), h_b.data(), kElems);
std::printf(" |A x - b|_inf = %.3e (bound %.3e): %s\n", axb,
kResidualBound, axb <= kResidualBound ? "pass" : "FAIL");
if (axb > kResidualBound) {
ok = false;
}
// Timed pass. Each run is a whole solve from x = 0, resets included, so
// the three rows differ only in who makes the per-iteration decision.
std::printf("\nTimed, mean of %d whole solves after %d warm-ups\n",
kTimedRuns, kWarmupRuns);
std::printf(" strategy total (ms) per iteration (ms)\n");
nvtxRangePushA("timed-host-loop");
const float msHost = timeKernel([&] {
resetState(d_x, d_iterations);
float r = 0.0f;
solveHostLoop(d_x, d_y, d_b, d_partials, d_residual, kOmegaFirst, &r);
});
nvtxRangePop();
std::printf(" host-loop-sync %10.3f %18.3f\n",
static_cast<double>(msHost),
static_cast<double>(msHost) / iter1);
nvtxRangePushA("timed-graph-loop");
const float msGraph = timeKernel([&] {
resetState(d_x, d_iterations);
float r = 0.0f;
solveGraphLoop(iterationExec, d_residual, &r);
});
nvtxRangePop();
std::printf(" graph-per-iteration %10.3f %18.3f\n",
static_cast<double>(msGraph),
static_cast<double>(msGraph) / iter2);
#if CUDART_VERSION >= 12040
nvtxRangePushA("timed-conditional");
const float msCond = timeKernel([&] {
resetState(d_x, d_iterations);
CUDA_CHECK(cudaGraphLaunch(condExec, nullptr));
});
nvtxRangePop();
std::printf(" conditional-graph %10.3f %18.3f\n",
static_cast<double>(msCond),
static_cast<double>(msCond) / iter3);
#endif // CUDART_VERSION >= 12040
// Update pass: change the omega baked into the iteration graph without
// re-instantiating. The recipe is capture again with the new value, then
// hand the topologically identical graph to cudaGraphExecUpdate. The
// result gate is real: an update that silently failed would leave the
// old omega running and the iteration count below would not move.
std::printf("\nUpdate pass: omega %.2f -> %.2f via cudaGraphExecUpdate\n",
static_cast<double>(kOmegaFirst),
static_cast<double>(kOmegaSecond));
// snippet: exec-update
nvtxRangePushA("recapture");
cudaGraph_t updatedGraph = captureIterationGraph(d_x, d_y, d_b, d_partials,
d_residual, kOmegaSecond);
nvtxRangePop();
nvtxRangePushA("exec-update");
cudaGraphExecUpdateResultInfo info = {};
CUDA_CHECK(cudaGraphExecUpdate(iterationExec, updatedGraph, &info));
nvtxRangePop();
CUDA_CHECK(cudaGraphDestroy(updatedGraph));
if (info.result != cudaGraphExecUpdateSuccess) {
std::fprintf(stderr, "cudaGraphExecUpdate result %d, not success\n",
static_cast<int>(info.result));
ok = false;
}
// end snippet
nvtxRangePushA("graph-after-update");
resetState(d_x, d_iterations);
float rSlow = 0.0f;
const int iterSlow = solveGraphLoop(iterationExec, d_residual, &rSlow);
nvtxRangePop();
CUDA_CHECK(cudaMemcpy(h_x2.data(), d_x, bytes, cudaMemcpyDeviceToHost));
const double axbSlow = trueResidual(h_x2.data(), h_b.data(), kElems);
std::printf(
" omega %.2f: %d iterations (was %d), residual %.3e, "
"|A x - b|_inf %.3e: %s\n",
static_cast<double>(kOmegaSecond), iterSlow, iter2,
static_cast<double>(rSlow), axbSlow,
axbSlow <= kResidualBound ? "pass" : "FAIL");
if (iterSlow >= kMaxIterations || rSlow > kTolerance ||
axbSlow > kResidualBound) {
ok = false;
}
if (iterSlow == iter2) {
std::fprintf(stderr,
"omega update had no visible effect: still %d "
"iterations\n",
iterSlow);
ok = false;
}
#if CUDART_VERSION >= 12040
CUDA_CHECK(cudaGraphExecDestroy(condExec));
CUDA_CHECK(cudaGraphDestroy(condGraph));
#endif // CUDART_VERSION >= 12040
CUDA_CHECK(cudaGraphExecDestroy(iterationExec));
CUDA_CHECK(cudaGraphDestroy(iterationGraph));
CUDA_CHECK(cudaFree(d_x));
CUDA_CHECK(cudaFree(d_y));
CUDA_CHECK(cudaFree(d_b));
CUDA_CHECK(cudaFree(d_partials));
CUDA_CHECK(cudaFree(d_residual));
CUDA_CHECK(cudaFree(d_iterations));
if (!ok) {
return EXIT_FAILURE;
}
std::printf("\nall gates passed\n");
return EXIT_SUCCESS;
}