  This paper proposes a versatile high-performance execution model, inspired by
systolic arrays, for memory-bound regular kernels running on CUDA-enabled GPUs.
We formulate a systolic model that shifts partial sums by CUDA warp primitives
for the computation. We also employ register files as a cache resource in order
to operate the entire model efficiently. We demonstrate the effectiveness and
versatility of the proposed model for a wide variety of stencil kernels that
appear commonly in HPC, and also convolution kernels (increasingly important in
deep learning workloads). Our algorithm outperforms the top reported
state-of-the-art stencil implementations, including implementations with
sophisticated temporal and spatial blocking techniques, on the two latest
Nvidia architectures: Tesla V100 and P100. For 2D convolution of general filter
sizes and shapes, our algorithm is on average 2.5x faster than Nvidia's NPP on
V100 and P100 GPUs.
