This issue appears to be distinct from #2900, as #2900 does not cause problems on my machine. Also, this problem is that the result is incorrect (and the termination code is undocumented); no error or warning is raised.

Unfortunately, I don't have a minimal example; just an example 107 x 107 system. The json for the matrix M and vector r is available at:
https://www.dropbox.com/s/nxcytfeqf5py7km/bad_Mr.json?dl=0
but I can paste it here if that would help.

The data can be loaded with the commands:

with open("bad_Mr.json","r") as f:
    d = json.load(f)
M = np.array(d["M"])
r = np.array(d["r"])

and the problem exposed with:

v = sp.linalg.solve(M,r)
v2 = sp.linalg.lstsq(M,r)[0]
v3 = sp.sparse.linalg.spsolve(M,r)
v4 = sp.sparse.linalg.lsqr(M,r)[0]

print(np.max(np.abs(M.dot(v)-r)))
print(np.max(np.abs(M.dot(v2)-r)))
print(np.max(np.abs(M.dot(v3)-r)))
print(np.max(np.abs(M.dot(v4)-r)))

which prints:

2.54658516496e-11
8.0035533756e-11
1.09139364213e-10
0.244709761602

The matrix is not of full row rank, so solve complains, but the solution it produces is fine.

Why can't lsqr handle it?