def solve_factorized_schur(z, Fval, LU, G, A):
    M, N=G.shape
    P, N=A.shape

    """Total number of inequality constraints"""
    m = M    

    """Primal variable"""
    x = z[0:N]

    """Multiplier for equality constraints"""
    nu = z[N:N+P]

    """Multiplier for inequality constraints"""
    l = z[N+P:N+P+M]

    """Slacks"""
    s = z[N+P+M:]

    """Dual infeasibility"""
    rd = Fval[0:N]
    
    """Primal infeasibility"""
    rp1 = Fval[N:N+P]
    rp2 = Fval[N+P:N+P+M]

    """Centrality"""
    rc = Fval[N+P+M:]

    """Sigma matrix"""
    SIG = diags(l/s, 0)

    """Assemble right hand side of augmented system"""
    r1 = rd + mydot(G.T, mydot(SIG, rp2)) - mydot(G.T, rc/s)
    r2 = rp1

    """Unpack LU-factors"""
    LU_S, LU_H = LU

    """Assemble right hand side for normal equation"""
    b = r2 - mydot(A, mysolve(LU_H, r1))  

    """Solve for dnu"""
    dnu = mysolve(LU_S, b)
       
    """Solve for dx"""
    dx = mysolve(LU_H, -(r1 + mydot(A.T, dnu)))    
    
    """Obtain search directions for l and s"""
    ds = -rp2 - mydot(G, dx)
    dl = -mydot(SIG, ds) - rc/s

    dz = np.hstack((dx, dnu, dl, ds))
    return dz
Total number of inequality constraints