code: def solve_factorized_aug(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)
  
      """RHS for condensed system"""
      b1 = -rd - mydot(G.T, mydot(SIG, rp2)) + mydot(G.T, rc/s)
      b2 = -rp1
      b = np.hstack((b1, b2))
      dxnu = mysolve(LU, b)
      dx = dxnu[0:N]
      dnu = dxnu[N:]
  
      """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