Describe your issue.

When calling iterative methods for linear systems from the scipy.sparse.linalg module (like cg or gmres) with a right hand side b of shape (N,1) the returned solution x is of the wrong shape (N,).

I tested for cg, cgs and gmres but it most likely affects some other iterative solvers as well.

I've tracked down the error to scipy.sparse.linalg._isolve.utils.make_system. There the function postprocess is defined which should ensure the correct shape of the returned solution.
While in the original call of make_system (and hence the definition of postprocess) b has the right shape (N,1). But at the time when postprocess is called the shape of b has magically changed to (N,). Hence, the solution is returned as (N,) vector.

Reproducing Code Example

import numpy as np
from scipy.sparse.linalg import cg

A = np.eye(10)
b = np.zeros((10,1))
x0 = np.zeros((10,1))

print("Without x0:")
print(f"\tb.shape  = {b.shape}")
x,flag = cg(A,b)
print(f"\tx.shape  = {x.shape}")

print("With x0:")
print(f"\tb.shape  = {b.shape}")
print(f"\tx0.shape = {x0.shape}")
x,flag = cg(A,b,x0)
print(f"\tx.shape  = {x.shape}")

Error message

Without x0:
        b.shape  = (10, 1)
        x.shape  = (10,)
With x0:
        b.shape  = (10, 1)
        x0.shape = (10, 1)
        x.shape  = (10,)

SciPy/NumPy/Python version information

1.8.0 1.22.3 sys.version_info(major=3, minor=9, micro=7, releaselevel='final', serial=0)