Well, what would make sense from a gufunc broadcasting point of view is to have a single function that solved things like `c = (.., M)` and `b = (..., M)`, and if you wanted to solve for multiple RHS pass in something like `c = (M)` and `b = (N, M)`, if you wanted multiple LHS pass in `c = (P, M)` and `b = (M)`, and if you wanted both then go with `c = (P, 1, M)` and `b = (N, M)`, which would give you a return `x = (P, N, M)` where item `x[p, n]` would be the solution with the `p`-th LHS and the `n`-th RHS. But GESV follows the opposite convention for the RHS, and so the efficient handling of multiple RHS requires that `b = (M, N)`, and that's what `solve` has been doing for ages, and sticking to the convention basically ruins having a nice consistent interface. If you look at [the Python source code of `np.linalg.solve`](https://github.com/jaimefrio/numpy/blob/master/numpy/linalg/linalg.py#L359), there are two different gufuncs called from `solve`, one for single RHS (`solve1`), other for multiple(`solve`), and the logic to select one or the other is not entirely satisfactory.