1: #pragma once
  2: #include <petsc/private/taoimpl.h>

  4: /*
  5:  Context for Primal-Dual Interior-Point Method
  6:  See the document pdipm.pdf
  7: */

  9: typedef struct {
 10:   /* Sizes (n = local, N = global) */
 11:   PetscInt nx, Nx;           /* Decision variables nx = nxfixed + nxub + nxlb + nxbox + nxfree */
 12:   PetscInt nxfixed, Nxfixed; /* Fixed decision variables */
 13:   PetscInt nxlb, Nxlb;       /* Decision variables with lower bounds only */
 14:   PetscInt nxub, Nxub;       /* Decision variables with upper bounds only */
 15:   PetscInt nxbox, Nxbox;     /* Decision variables with box constraints */
 16:   PetscInt nxfree, Nxfree;   /* Free variables */
 17:   PetscInt ng, Ng;           /* user equality constraints g(x) = 0. */
 18:   PetscInt nh, Nh;           /* user inequality constraints h(x) >= 0. */
 19:   PetscInt nce, Nce;         /* total equality constraints. nce = ng + nxfixed */
 20:   PetscInt nci, Nci;         /* total inequality constraints nci = nh + nxlb + nxub + 2*nxbox */
 21:   PetscInt n, N;             /* Big KKT system size n = nx + nce + 2*nci */

 23:   /* Vectors */
 24:   Vec X;       /* R^n   - Big KKT system vector [x; lambdae; lambdai; z] */
 25:   Vec x;       /* R^nx - work vector, same layout as tao->solution */
 26:   Vec lambdae; /* R^nce - vector, shares local arrays with X */
 27:   Vec lambdai; /* R^nci - vector, shares local arrays with X */
 28:   Vec z;       /* R^nci - vector, shares local arrays with X */

 30:   /* Work vectors */
 31:   Vec lambdae_xfixed; /* Equality constraints lagrangian multiplier vector for fixed variables */
 32:   Vec lambdai_xb;     /* User inequality constraints lagrangian multiplier vector */

 34:   /* Lagrangian equality and inequality Vec */
 35:   Vec ce, ci; /* equality and inequality constraints */

 37:   /* Offsets for subvectors */
 38:   PetscInt off_lambdae, off_lambdai, off_z;

 40:   /* Scalars */
 41:   PetscReal L;     /* Lagrangian = f(x) - lambdae^T*ce(x) - lambdai^T*(ci(x) - z) - mu*sum_{i=1}^{Nci}(log(z_i)) */
 42:   PetscReal gradL; /* gradient of L w.r.t. x */

 44:   /* Matrices */
 45:   Mat Jce_xfixed; /* Jacobian of equality constraints cebound(x) = J(nxfixed) */
 46:   Mat Jci_xb;     /* Jacobian of inequality constraints Jci = [tao->jacobian_inequality ; J(nxub); J(nxlb); J(nxbx)] */
 47:   Mat K;          /* KKT matrix */

 49:   /* Parameters */
 50:   PetscReal mu;               /* Barrier parameter */
 51:   PetscReal mu_update_factor; /* Multiplier for mu update */
 52:   PetscReal deltaw;
 53:   PetscReal lastdeltaw;
 54:   PetscReal deltac;

 56:   /* Tolerances */

 58:   /* Index sets for types of bounds on variables */
 59:   IS isxub;    /* Finite upper bound only -inf < x < ub   */
 60:   IS isxlb;    /* Finite lower bound only  lb <= x < inf  */
 61:   IS isxfixed; /* Fixed variables          lb  = x = ub   */
 62:   IS isxbox;   /* Boxed variables          lb <= x <= ub  */
 63:   IS isxfree;  /* Free variables         -inf <= x <= inf */

 65:   /* Index sets for PC fieldsplit */
 66:   IS is1, is2;

 68:   /* Options */
 69:   PetscBool monitorkkt;          /* Monitor KKT */
 70:   PetscReal push_init_slack;     /* Push initial slack variables (z) away from bounds */
 71:   PetscReal push_init_lambdai;   /* Push initial inequality variables (lambdai) away from bounds */
 72:   PetscBool solve_reduced_kkt;   /* Solve Reduced KKT with fieldsplit */
 73:   PetscBool solve_symmetric_kkt; /* Solve non-reduced symmetric KKT system */
 74:   PetscBool kkt_pd;              /* Add deltaw and deltac shifts to make KKT matrix positive definite */

 76:   SNES snes;                                     /* Nonlinear solver */
 77:   Mat  jac_equality_trans, jac_inequality_trans; /* working matrices */

 79:   PetscReal obj; /* Objective function */

 81:   /* Offsets for parallel assembly */
 82:   PetscInt *nce_all;
 83: } TAO_PDIPM;