Reductions in an interior-point solve: a worked solution

The worked solution to the first exercise of tutorial 1. Attempt it there first.

Three per-iteration quantities of an interior-point solver, each a single mapreduce over the device, with no temporary array: the primal infeasibility maximum |c|, the complementarity max |x_i z_i|, and the fraction-to-boundary step, the largest α keeping x + α dx positive.

using CUDA

c = CUDA.randn(Float64, 10^5)
xv = CUDA.rand(Float64, 10^5) .+ 0.1
zv = CUDA.rand(Float64, 10^5) .+ 0.1
dx = CUDA.randn(Float64, 10^5)

inf_pr = mapreduce(abs, max, c)

# mapreduce takes several arrays and walks them together
compl = mapreduce((x, z) -> abs(x * z), max, xv, zv)

# Inf for entries that do not constrain the step: the identity for min
alpha_max = mapreduce((x, d) -> d < 0 ? -x / d : Inf, min, xv, dx)

(inf_pr = inf_pr, complementarity = compl, alpha_max = alpha_max)

(inf_pr = 4.790118053726513, complementarity = 1.2048985717299163, alpha_max = 0.033370771696594755)

Each is one pass over the device, and none of them builds an intermediate array. This is how a GPU-resident interior-point method keeps its work on the device: the iterates never come back to the host, and the scalars that steer the algorithm are computed where the data already is.

---

This notebook was generated using Literate.jl.