  We introduce an efficient first-order primal-dual method for the solution of
nonsmooth PDE-constrained optimization problems. We achieve this efficiency
through not solving the PDE or its linearisation on each iteration of the
optimization method. Instead, we run the method interwoven with a simple
conventional linear system solver (Jacobi, Gauss-Seidel, conjugate gradients),
always taking only one step of the linear system solver for each step of the
optimization method. The control parameter is updated on each iteration as
determined by the optimization method. We prove linear convergence under a
second-order growth condition, and numerically demonstrate the performance on a
variety of PDEs related to inverse problems involving boundary measurements.
