docs: def primal_dual_solve(func, x0, Dfunc, A, b, G, h, args=(), tol=1e-10,
                        maxiter=100, show_progress=True, full_output=False):
      """Wrap calls to function and Jacobian"""
      fcalls, func = wrap_function(func, args)
      Dfcalls, Dfunc = wrap_function(Dfunc, args)
  
      M, N = G.shape
      P, N = A.shape
  
      """Total number of inequality constraints"""
      m = M
  
      def gap(z):
          r"""Gap-function"""
          l = z[N+P:N+P+M]
          s = z[N+P+M:]
          return mydot(l, s)/m
  
      def centrality(z):
          r"""Centrality function"""
          l = z[N+P:N+P+M]
          s = z[N+P+M:]
          return np.min(l*s)
  
      def KKT(z, sigma=0.0):
          r"""KKT system (possible perturbed)."""
  
          """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:]
  
          """Evaluate objective function"""
          F = func(x)
  
          """Dual infeasibility"""
          rdual = F+mydot(A.transpose(), nu)+mydot(G.transpose(), l)
  
          """Primal infeasibilities"""
          rprim1 = mydot(A, x)-b
          rprim2 = mydot(G, x)-h+s
  
          """Complementary slackness (perturbed)"""
          mu = gap(z)
          rslack = l*s-sigma*mu
  
          return np.hstack((rdual, rprim1, rprim2, rslack))
  
      def step_fast(z, KKTval, LU, G, A, mu, beta, gamma, alpha0):
          r"""Affine scaling step."""
          dz = solve_factorized(z, KKTval, LU, G, A)
  
          """Reduce step length until slacks s and multipliers l are positive"""
          alpha = 1.0*alpha_0
          k = 0
          for k in range(10):
              z_new = z + alpha*dz
              if np.all( z_new[N+P:] > 0.0 ):
                  break
              alpha *= 0.5
              k += 1
          if k == 10 - 1:
              raise RuntimeError("Maximum steplength reduction reached")
  
          """Reduce step length until iterates lie in correct neighborhood"""
          for k in range(10):
              z_new = z + alpha*dz
              KKTval_new = KKT(z_new)
              dual = mynorm(KKTval_new[0:N])
              prim = mynorm(KKTval_new[N:N+P+M])
              mu_new = gap(z_new)
              cent_new = centrality(z_new)
  
              if (dual <= beta*mu_new and prim <= beta*mu_new and
                  cent_new >= gamma*mu_new):
                  break
              alpha *= 0.5
              # alpha *= 0.95
          if k == 10 - 1:
              raise RuntimeError("Maximum steplength reduction reached")
          return z_new, mu_new
  
      def step_safe(z, KKTval, LU, G, A, mu, beta, gamma):
          r"""Centering step."""
          dz = solve_factorized(z, KKTval, LU, G, A)
  
          """Reduce step length until slacks s and multipliers l are positive"""
          alpha = 1.0
          k = 0
          for k in range(10):
              z_new = z + alpha*dz
              if np.all( z_new[N+P:] > 0.0 ):
                  break
              alpha *= 0.5
              k += 1
          if k == 10 - 1:
              raise RuntimeError("Maximum steplength reduction (pos.) reached")
  
          """Reduce step length until iterates lie in correct neighborhood
          and mu fulfills Armijo condition"""
          k = 0
          for k in range(10):
              z_new = z+alpha*dz
              KKTval_new = KKT(z_new, sigma=SIGMA)
              dual = mynorm(KKTval_new[0:N])
              prim = mynorm(KKTval_new[N:N+P+M])
              mu_new = gap(z_new)
              cent_new = centrality(z_new)
  
              if (dual <= beta*mu_new and prim <= beta*mu_new and
                  cent_new >= gamma*mu_new and
                  mu_new<=(1.0-KAPPA*alpha*(1.0-SIGMA))*mu):
                  break
              alpha *= 0.5
          if k == 10 - 1:
              raise RuntimeError("Maximum steplength reduction reached")
          return z_new, mu_new
  
      """INITIALIZATION"""
  
      """Initial Slacks for inequality constraints"""
      s0 = -1.0*(mydot(G, x0)-h)
  
      """Initial multipliers for inequality constraints"""
      l0 = 1.0*np.ones(M)
  
      """Initial multipliers for equality constraints"""
      nu0 = np.zeros(P)
  
      """Initial point"""
      z0 = np.hstack((x0, nu0, l0, s0))
  
      """Initial KKT-values"""
      KKTval0 = KKT(z0, sigma=0.0)
      mu0 = gap(z0)
      dual0 = mynorm(KKTval0[0:N])
      prim0 = mynorm(KKTval0[N:N+P+M])
  
      """Initial neighborhood"""
      beta = BETA * np.sqrt(dual0**2 + prim0**2)/mu0
      gamma = 1.0 * GAMMA_MAX
  
      """Number of fast steps"""
      t = 0
  
      """Number of iterations"""
      n = 0
  
      """Dummy variable for step type"""
      step_type = " "
  
      if show_progress:
          print("%s %s %s %s %s %s %s" %("iter", "gap", "dual", "primal",
                                         "min", "max", "step"))
      """MAIN LOOP"""
      z = z0
      x = z0[0:N]
      KKTval = KKTval0
      dual = dual0
      prim = prim0
      mu = mu0
      Dfunc_val = Dfunc(x)
      LU = factor(z, Dfunc_val, G, A)
  
      if full_output:
          info = {'z': []}
          info['z'].append(z)
  
      for n in range(maxiter):
          if show_progress:
              l=z[N+P:N+P+M]
              s=z[N+P+M:]
              print("%i %.6e %.6e %.6e %.6e %.6e %s" %(n+1, mu, dual, prim,
                                                       np.min(l*s), np.max(l*s),
                                                       step_type))
          """Attempt fast step"""
          beta_new = (1.0 + GAMMA_BAR**(t+1)) * beta
          gamma_new = GAMMA_MIN + GAMMA_BAR**(t+1)*(GAMMA_MAX - GAMMA_MIN)
          alpha_0 = 1.0 - np.sqrt(mu)/GAMMA_BAR**t
  
          if alpha_0 > 0.0:
              z_new, mu_new = step_fast(z, KKTval, LU, G, A, mu,
                                        beta_new, gamma_new, alpha_0)
              if mu_new < RHO * mu:
                  """Fast successful"""
                  z = z_new
                  mu = mu_new
                  beta = beta_new
                  gamma = gamma_new
                  t += 1
                  step_type = "f"
              else:
                  """Perturbed right-had side"""
                  KKTval_pert = 1.0*KKTval
                  KKTval_pert[N+P+M:] -= SIGMA * mu
                  z, mu = step_safe(z, KKTval_pert, LU, G, A, mu, beta, gamma)
                  step_type = "s"
          else:
              """Perturbed right-hand side"""
              KKTval_pert = 1.0*KKTval
              KKTval_pert[N+P+M:] -= SIGMA * mu
              z, mu = step_safe(z, KKTval_pert, LU, G, A, mu, beta, gamma)
              step_type = "s"
  
          """Compute new iterates"""
          KKTval = KKT(z, sigma=0.0)
          dual = mynorm(KKTval[0:N])
          prim = mynorm(KKTval[N:N+P+M])
          x = z[0:N]
          Dfunc_val = Dfunc(x)
          LU = factor(z, Dfunc_val, G, A)
          if full_output:
              info['z'].append(z)
          if mu < tol and dual < tol and prim < tol:
              break
      if n == maxiter - 1:
          raise RuntimeError("Maximum number of iterations reached")
  
      if show_progress:
          l=z[N+P:N+P+M]
          s=z[N+P+M:]
          print("%i %.6e %.6e %.6e %.6e %.6e %s" %(n+1, mu, dual, prim,
                                                   np.min(l*s), np.max(l*s),
                                                   step_type))
      if full_output:
          return z[0:N], info
      else:
          return z[0:N]