from scipy.sparse.linalg import bicg

def bicg_solve(M, B):
    X, info = zip(*(bicg(M, b) for b in B.T))
return np.transpose(X), info
Iterative solving of sparse systems of linear equations with (M, N) right-hand size matrix