It may require thinking about the UI a little harder, but I don't think you need to enforce 1-dimensionality. Instead, you could accept `c`s of shape `(..., M)` and `b`s of shape `(..., M)` or `(..., M, N)`. If you can identify where in `c` and `b` is the dimension of shape `M`, you would simply have to take the fft along that axis, and let broadcasting do its magic. That way you could solve not only for `N` different RHS, but also for `P` different `LHS` by passing `c` of shape `(P, M)` and `b` of `(M, N)`... Sort of what happens in `np.linalg.solve` in 1.8.