[Python] Solve indefinite linear systems using NumPy only

Suppose we are going to solve for variables $ \bf x$ $\in {\mathbb R}^{n} $
$${\bf Ax}={\bf b}$$
where $ \bf A$ $\in {\mathbb R}^{n_{\rm eq} \times n} $ is a matrix and $ \bf b$ $\in {\mathbb R}^{n_{\rm eq}} $ is a vector. For simplicity, ${\rm rank}( {\bf A}) = n_{\rm eq} $ is assumed hereafter.

import numpy as np

def solve_base(A,b,tol=1.0e-8):
    '''
    Solve Ax=b to obtain a base solution.
    (input): A[neq,nv]<float>
             b[neq]<float>
    (output): x[nv,nv-neq]: base solution
    '''
    E = np.concatenate([A,b.reshape([len(b),1]),np.zeros_like(A)],axis=1)
    x_order = np.arange(A.shape[1],dtype=int)
    for i in range(A.shape[0]):
        if np.all(np.abs(E[i:,i])<tol) and np.all(np.abs(E[i,i:])<tol):
            E[i,i] = 1
            E[i,A.shape[1]+1] = 1
        else:
            if np.all(np.abs(E[i:,i])<tol): # swap columns
                swap_index = np.where(np.abs(E[i,i+1:])>tol)[0][0]+i+1
                E[:,[i,swap_index]] = E[:,[swap_index,i]]
                x_order[[i,swap_index]] = x_order[[swap_index,i]]
            if abs(E[i,i])<tol: # swap rows
                swap_index = np.where(np.abs(E[:,i])>tol)[0][0]
                E[[i,swap_index]] = E[[swap_index,i]]
            E[i] /= E[i,i]
            for j in range(i+1,A.shape[0]):
                E[j] -= E[j,i] * E[i]

    E = np.r_[E,np.zeros([A.shape[1]-A.shape[0],E.shape[1]])]
    for i in range(A.shape[1]-A.shape[0]):
        E[A.shape[0]+i,A.shape[0]+i] = 1
        E[A.shape[0]+i,A.shape[0]+A.shape[1]+1+i] = 1
    for i in range(1,E.shape[0]):
        for j in range(i):
            E[j,i:] -= E[i,i:]*E[j,i]

    x = E[x_order,A.shape[1]:]
        
    return x

x = solve_base(np.array([[1,1],[2,2]],dtype=float),np.array([2,4],dtype=float))
'''
x = [[2,-1,0],
      0, 1,0]]
'''