splu does not work for non-vector inputs

[pv]

splu does not actually work when called with a non-vector rhs:

import numpy as np
from scipy.sparse.linalg import splu
A = np.random.rand(2,2)
b = np.random.rand(2,2)
x2 = splu(A).solve(b)
x1 = np.linalg.solve(A, b)
print(abs(x1 - x2).max())
# -> 2.2361768860074456

works for vectors:

x2 = np.array([splu(A).solve(b[:,j]) for j in range(2)]).T
print(abs(x1 - x2).max())
# -> 0

Raises an error for some rhs shapes:

b = np.random.rand(2,1)
x2 = splu(A).solve(b)
# -> SystemError

[pv]

Not a recent regression, also 0.7.2 does the same.