code: def solve_factorized_schur(z, Fval, LU, G, A):
      M, N=G.shape
      P, N=A.shape
  
      
      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
docs: Total number of inequality constraints