docs: def factor_schur(z, DPhival, G, A):
      M, N = G.shape
      P, N = A.shape
      """Multiplier for inequality constraints"""
      l = z[N+P:N+P+M]
  
      """Slacks"""
      s = z[N+P+M:]
  
      """Sigma matrix"""
      SIG = diags(l/s, 0)
  
      """Augmented Jacobian"""
      H = DPhival + mydot(G.T, mydot(SIG, G))
  
      """Factor H"""
      LU_H = myfactor(H)
  
      """Compute H^{-1}A^{T}"""
      HinvAt = mysolve(LU_H, A.T)
  
      """Compute Schur complement AH^{-1}A^{T}"""
      S = mydot(A, HinvAt)
  
      """Factor Schur complement"""
      LU_S = myfactor(S)
  
      LU = (LU_S, LU_H)
      return LU