docs: def solve_factorized_aug(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)
  
      """LU is actually the augmented system W"""
      W = LU
  
      b1 = -rd - mydot(G.T, mydot(SIG, rp2)) + mydot(G.T, rc/s)
      b2 = -rp1
      b = np.hstack((b1, b2))
  
      """Prepare iterative solve via MINRES"""
      sign = np.zeros(N+P)
      sign[0:N//2] = 1.0
      sign[N//2:] = -1.0
      T = diags(sign, 0)        
     
      """Change rhs"""
      b_new = mydot(T, b)
  
      dW = np.abs(W.diagonal())
      dPc = np.ones(W.shape[0])
      ind = (dW > 0.0)
      dPc[ind] = 1.0/dW[ind]
      Pc = diags(dPc, 0)    
      dxnu, info = minres(W, b_new, tol=1e-10, M=Pc)
      
      # dxnu = solve(J, b)
      dx = dxnu[0:N]
      dnu = dxnu[N:]
  
      """Obtain search directions for l and s"""
      ds = -rp2 - mydot(G, dx)
      # ds = s*ds
      # SIG = np.diag(l/s)
      dl = -mydot(SIG, ds) - rc/s
  
      dz = np.hstack((dx, dnu, dl, ds))
      return dz