# Watch the Residual Fall

*Task 2 of 5 · [Iterative Linear Solvers](https://gpu.rocks/learn/iterative-solvers-e73b8e1f.md) · GPU.js Learn*

Iterating is easy; knowing when to stop is the skill. The honest measure is the
**residual**: take the current guess, put it back into the equation, and see
how badly it is violated. For "every cell equals the average of its four neighbours" that
is the 5-point Laplacian — `left + right + up + down − 4·centre` — the same
stencil Reaction–Diffusion uses to diffuse, borrowed here as a scorecard. Zero where the
equation holds, large where it does not.

Boundary cells have no equation to violate: they are given, not solved. Their residual
is `0` by definition rather than by accident, and saying so in the kernel keeps
the edge from polluting the score forever.

A grid of numbers is not a progress report, so collapse it to one: the root-mean-square
over the grid. Totalling a grid *on* the GPU is the halving ladder Reductions
builds; at 1,024 cells the read-back is cheaper than the ladder, so this sum happens in
plain JavaScript.

## Goal

**Goal:** complete the `residual` kernel and the
`rmsOf` helper. The sweep loop is already wired and will print the residual
every 10 sweeps — you should watch it fall by a factor of about 36.

## Requirements

- The kernel returns `0` for boundary cells
- Interior cells return `left + right + up + down − 4·centre` — a sum, with no division
- `rmsOf(grid)` returns `Math.sqrt(sum of every cell squared / (32 × 32))`

## Hint 1 — the kernel is the stencil, unaveraged

Same five reads as the Jacobi sweep, assembled differently: the update divides
the neighbour sum by 4, the residual subtracts `4 × centre` from it.

```js
return u[y][x - 1] + u[y][x + 1] + u[y - 1][x]
     + u[y + 1][x] - 4 * u[y][x];
```

## Hint 2 — the reduction is ordinary JavaScript

`residual(u)` hands back a plain 2D array of numbers, so:

```js
let sum = 0;
for (let y = 0; y < 32; y++) {
  for (let x = 0; x < 32; x++) sum += grid[y][x] * grid[y][x];
}
return Math.sqrt(sum / (32 * 32));
```

Square, mean, root — in that order. Divide by every cell in the grid, not just the
interior ones.

## Same idea elsewhere

Every production solver stops on a residual, not on a sweep count, and every one
of them argues about how often to measure it: the norm needs a reduction across the whole
device and then a read-back to the host, which is a synchronisation point in CUDA, WebGPU
and MPI alike. Checking every iteration can cost more than the iterations do — the usual
answer is exactly what this task does, sample it every few sweeps.

## Starter code

```js
// How wrong is the current guess? Plug it back into the equation.
const gpu = new GPU({ mode });
const SWEEPS = 60;

const sweep = gpu.createKernel(function (u) {
  const x = this.thread.x;
  const y = this.thread.y;
  if (x === 0 || y === 0 || x === this.constants.size - 1 || y === this.constants.size - 1) {
    return u[y][x];
  }
  return (u[y][x - 1] + u[y][x + 1] + u[y - 1][x] + u[y + 1][x]) / 4;
}, { output: [32, 32], constants: { size: 32 } });

const residual = gpu.createKernel(function (u) {
  const x = this.thread.x;
  const y = this.thread.y;
  // TODO 1: a boundary cell has no equation to violate — return 0.
  // TODO 2: every other cell returns
  //         left + right + up + down - 4 * centre.
  return 1;
}, { output: [32, 32], constants: { size: 32 } });

function rmsOf(grid) {
  // TODO 3: square every cell, take the mean over all 32 × 32 of them,
  //         then the square root.
  return 0;
}

let u = plate;
for (let k = 0; k <= SWEEPS; k++) {
  if (k % 10 === 0) console.log('sweep', k, '— RMS residual', rmsOf(await residual(u)));
  u = await sweep(u);
}
```

---

Interactive version: https://gpu.rocks/learn/iterative-solvers-e73b8e1f/2

[Previous task](https://gpu.rocks/learn/iterative-solvers-e73b8e1f/1.md) · [Next task](https://gpu.rocks/learn/iterative-solvers-e73b8e1f/3.md)