Add a `Axpy` that implements y <- y + ax.

This is very useful for accumulations to compute, e.g. centers.