Bonus if you can use Subspace-BB (Barzilai-Borwein) on what is left after activating your constraints. That stuff is magical. Also because you can throw it at linear systems that are actually FFT-based time/space variant convolutions or similar things where you never, ever want so much as touch the sparse matrix of that linear system.

Y not touch the sparse matrix of that system?

Because that'd be equivalent to computing the convolution explicitly. The whole approach is based on tricking Overlapp-Add to support compactly supported interpolation kernels between representative convolution kernels. Easiest would be bilinear interpolation between between the surrounding representative kernels, which mostly translated to using the bilinear weights before using the fft convolution, and deferring the summation of the interpolation to be fused with the addition done due to overlap-add.

Yea...if you're an OR/Applied Math guru, there is some magical stuff you can do.