docs: def solve_full(z, Fval, DPhival, 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 = np.diag(l/s)
  
      """Condensed system"""
      if issparse(DPhival):
          if not issparse(A):
              A = csr_matrix(A)        
          H = DPhival + mydot(G.T, mydot(SIG, G))
          J = bmat([[H, A.T], [A, None]])
      else:
          if issparse(A):
              A = A.toarray()
          J = np.zeros((N+P, N+P))
          J[0:N, 0:N] = DPhival + mydot(G.T, mydot(SIG, G))            
          J[0:N, N:] = A.T
          J[N:, 0:N] = A
  
      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
      S = diags(sign, 0)
      J_new = mydot(S, csr_matrix(J))
      b_new = mydot(S, b)
  
      dJ_new = np.abs(J_new.diagonal())
      dPc = np.ones(J_new.shape[0])
      ind = (dJ_new > 0.0)
      dPc[ind] = 1.0/dJ_new[ind]
      Pc = diags(dPc, 0)    
      dxnu, info = minres(J_new, b_new, tol=1e-8, 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)
      dl = -mydot(SIG, ds) - rc/s
  
      dz = np.hstack((dx, dnu, dl, ds))
      return dz