xref: /petsc/src/tao/constrained/tutorials/ex1.c (revision 48a46eb9bd028bec07ec0f396b1a3abb43f14558)
1aad13602SShrirang Abhyankar /* Program usage: mpiexec -n 2 ./ex1 [-help] [all TAO options] */
2aad13602SShrirang Abhyankar 
3aad13602SShrirang Abhyankar /* ----------------------------------------------------------------------
4628da978Sresundermann min f(x) = (x0 - 2)^2 + (x1 - 2)^2 - 2*(x0 + x1)
5aad13602SShrirang Abhyankar s.t.  x0^2 + x1 - 2 = 0
6aad13602SShrirang Abhyankar       0  <= x0^2 - x1 <= 1
7aad13602SShrirang Abhyankar       -1 <= x0, x1 <= 2
8628da978Sresundermann -->
9628da978Sresundermann       g(x)  = 0
10628da978Sresundermann       h(x) >= 0
11628da978Sresundermann       -1 <= x0, x1 <= 2
12628da978Sresundermann where
13628da978Sresundermann       g(x) = x0^2 + x1 - 2
14628da978Sresundermann       h(x) = [x0^2 - x1
15628da978Sresundermann               1 -(x0^2 - x1)]
16aad13602SShrirang Abhyankar ---------------------------------------------------------------------- */
17aad13602SShrirang Abhyankar 
18aad13602SShrirang Abhyankar #include <petsctao.h>
19aad13602SShrirang Abhyankar 
20aad13602SShrirang Abhyankar static char help[] = "Solves constrained optimiztion problem using pdipm.\n\
21aad13602SShrirang Abhyankar Input parameters include:\n\
22aad13602SShrirang Abhyankar   -tao_type pdipm    : sets Tao solver\n\
238e58fa1dSresundermann   -no_eq             : removes the equaility constraints from the problem\n\
24f560b561SHong Zhang   -init_view         : view initial object setup\n\
25aad13602SShrirang Abhyankar   -snes_fd           : snes with finite difference Jacobian (needed for pdipm)\n\
26628da978Sresundermann   -snes_compare_explicit : compare user Jacobian with finite difference Jacobian \n\
27aad13602SShrirang Abhyankar   -tao_cmonitor      : convergence monitor with constraint norm \n\
28a5b23f4aSJose E. Roman   -tao_view_solution : view exact solution at each iteration\n\
29f560b561SHong Zhang   Note: external package MUMPS is required to run pdipm in parallel. This is designed for a maximum of 2 processors, the code will error if size > 2.\n";
30aad13602SShrirang Abhyankar 
31aad13602SShrirang Abhyankar /*
32628da978Sresundermann    User-defined application context - contains data needed by the application
33aad13602SShrirang Abhyankar */
34aad13602SShrirang Abhyankar typedef struct {
35aad13602SShrirang Abhyankar   PetscInt   n;  /* Global length of x */
36aad13602SShrirang Abhyankar   PetscInt   ne; /* Global number of equality constraints */
37aad13602SShrirang Abhyankar   PetscInt   ni; /* Global number of inequality constraints */
38f560b561SHong Zhang   PetscBool  noeqflag, initview;
39aad13602SShrirang Abhyankar   Vec        x, xl, xu;
40aad13602SShrirang Abhyankar   Vec        ce, ci, bl, bu, Xseq;
41aad13602SShrirang Abhyankar   Mat        Ae, Ai, H;
42aad13602SShrirang Abhyankar   VecScatter scat;
43aad13602SShrirang Abhyankar } AppCtx;
44aad13602SShrirang Abhyankar 
45aad13602SShrirang Abhyankar /* -------- User-defined Routines --------- */
46aad13602SShrirang Abhyankar PetscErrorCode InitializeProblem(AppCtx *);
47aad13602SShrirang Abhyankar PetscErrorCode DestroyProblem(AppCtx *);
48aad13602SShrirang Abhyankar PetscErrorCode FormFunctionGradient(Tao, Vec, PetscReal *, Vec, void *);
49aad13602SShrirang Abhyankar PetscErrorCode FormHessian(Tao, Vec, Mat, Mat, void *);
50aad13602SShrirang Abhyankar PetscErrorCode FormInequalityConstraints(Tao, Vec, Vec, void *);
51aad13602SShrirang Abhyankar PetscErrorCode FormEqualityConstraints(Tao, Vec, Vec, void *);
52aad13602SShrirang Abhyankar PetscErrorCode FormInequalityJacobian(Tao, Vec, Mat, Mat, void *);
53aad13602SShrirang Abhyankar PetscErrorCode FormEqualityJacobian(Tao, Vec, Mat, Mat, void *);
54aad13602SShrirang Abhyankar 
559371c9d4SSatish Balay PetscErrorCode main(int argc, char **argv) {
56aad13602SShrirang Abhyankar   Tao         tao;
57aad13602SShrirang Abhyankar   KSP         ksp;
58aad13602SShrirang Abhyankar   PC          pc;
59aad13602SShrirang Abhyankar   AppCtx      user; /* application context */
60f560b561SHong Zhang   Vec         x, G, CI, CE;
61f560b561SHong Zhang   PetscMPIInt size;
62661095bbSAlp Dener   TaoType     type;
63f560b561SHong Zhang   PetscReal   f;
64661095bbSAlp Dener   PetscBool   pdipm;
65aad13602SShrirang Abhyankar 
66327415f7SBarry Smith   PetscFunctionBeginUser;
679566063dSJacob Faibussowitsch   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
689566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
693c859ba3SBarry Smith   PetscCheck(size <= 2, PETSC_COMM_SELF, PETSC_ERR_WRONG_MPI_SIZE, "More than 2 processors detected. Example written to use max of 2 processors.");
70aad13602SShrirang Abhyankar 
719566063dSJacob Faibussowitsch   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "---- Constrained Problem -----\n"));
729566063dSJacob Faibussowitsch   PetscCall(InitializeProblem(&user)); /* sets up problem, function below */
739566063dSJacob Faibussowitsch   PetscCall(TaoCreate(PETSC_COMM_WORLD, &tao));
749566063dSJacob Faibussowitsch   PetscCall(TaoSetType(tao, TAOPDIPM));
759566063dSJacob Faibussowitsch   PetscCall(TaoSetSolution(tao, user.x));                 /* gets solution vector from problem */
769566063dSJacob Faibussowitsch   PetscCall(TaoSetVariableBounds(tao, user.xl, user.xu)); /* sets lower upper bounds from given solution */
779566063dSJacob Faibussowitsch   PetscCall(TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, (void *)&user));
78aad13602SShrirang Abhyankar 
79*48a46eb9SPierre Jolivet   if (!user.noeqflag) PetscCall(TaoSetEqualityConstraintsRoutine(tao, user.ce, FormEqualityConstraints, (void *)&user));
809566063dSJacob Faibussowitsch   PetscCall(TaoSetInequalityConstraintsRoutine(tao, user.ci, FormInequalityConstraints, (void *)&user));
819371c9d4SSatish Balay   if (!user.noeqflag) { PetscCall(TaoSetJacobianEqualityRoutine(tao, user.Ae, user.Ae, FormEqualityJacobian, (void *)&user)); /* equality jacobian */ }
829566063dSJacob Faibussowitsch   PetscCall(TaoSetJacobianInequalityRoutine(tao, user.Ai, user.Ai, FormInequalityJacobian, (void *)&user)); /* inequality jacobian */
839566063dSJacob Faibussowitsch   PetscCall(TaoSetTolerances(tao, 1.e-6, 1.e-6, 1.e-6));
849566063dSJacob Faibussowitsch   PetscCall(TaoSetConstraintTolerances(tao, 1.e-6, 1.e-6));
85aad13602SShrirang Abhyankar 
869566063dSJacob Faibussowitsch   PetscCall(TaoGetKSP(tao, &ksp));
879566063dSJacob Faibussowitsch   PetscCall(KSPGetPC(ksp, &pc));
889566063dSJacob Faibussowitsch   PetscCall(PCSetType(pc, PCCHOLESKY));
89aad13602SShrirang Abhyankar   /*
90aad13602SShrirang Abhyankar       This algorithm produces matrices with zeros along the diagonal therefore we use
9112d688e0SRylee Sundermann     MUMPS which provides solver for indefinite matrices
92aad13602SShrirang Abhyankar   */
93f560b561SHong Zhang #if defined(PETSC_HAVE_MUMPS)
949566063dSJacob Faibussowitsch   PetscCall(PCFactorSetMatSolverType(pc, MATSOLVERMUMPS)); /* requires mumps to solve pdipm */
95f560b561SHong Zhang #else
963c859ba3SBarry Smith   PetscCheck(size == 1, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Requires an external package that supports parallel PCCHOLESKY, e.g., MUMPS.");
97f560b561SHong Zhang #endif
989566063dSJacob Faibussowitsch   PetscCall(KSPSetType(ksp, KSPPREONLY));
999566063dSJacob Faibussowitsch   PetscCall(KSPSetFromOptions(ksp));
100aad13602SShrirang Abhyankar 
1019566063dSJacob Faibussowitsch   PetscCall(TaoSetFromOptions(tao));
1029566063dSJacob Faibussowitsch   PetscCall(TaoGetType(tao, &type));
1039566063dSJacob Faibussowitsch   PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOPDIPM, &pdipm));
104661095bbSAlp Dener   if (pdipm) {
1059566063dSJacob Faibussowitsch     PetscCall(TaoSetHessian(tao, user.H, user.H, FormHessian, (void *)&user));
106f560b561SHong Zhang     if (user.initview) {
1079566063dSJacob Faibussowitsch       PetscCall(TaoSetUp(tao));
1089566063dSJacob Faibussowitsch       PetscCall(VecDuplicate(user.x, &G));
1099566063dSJacob Faibussowitsch       PetscCall(FormFunctionGradient(tao, user.x, &f, G, (void *)&user));
1109566063dSJacob Faibussowitsch       PetscCall(FormHessian(tao, user.x, user.H, user.H, (void *)&user));
1119566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_WORLD));
11263a3b9bcSJacob Faibussowitsch       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nInitial point X:\n"));
1139566063dSJacob Faibussowitsch       PetscCall(VecView(user.x, PETSC_VIEWER_STDOUT_WORLD));
11463a3b9bcSJacob Faibussowitsch       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nInitial objective f(x) = %g\n", (double)f));
11563a3b9bcSJacob Faibussowitsch       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nInitial gradient and Hessian:\n"));
1169566063dSJacob Faibussowitsch       PetscCall(VecView(G, PETSC_VIEWER_STDOUT_WORLD));
1179566063dSJacob Faibussowitsch       PetscCall(MatView(user.H, PETSC_VIEWER_STDOUT_WORLD));
1189566063dSJacob Faibussowitsch       PetscCall(VecDestroy(&G));
1199566063dSJacob Faibussowitsch       PetscCall(FormInequalityJacobian(tao, user.x, user.Ai, user.Ai, (void *)&user));
1209566063dSJacob Faibussowitsch       PetscCall(MatCreateVecs(user.Ai, NULL, &CI));
1219566063dSJacob Faibussowitsch       PetscCall(FormInequalityConstraints(tao, user.x, CI, (void *)&user));
12263a3b9bcSJacob Faibussowitsch       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nInitial inequality constraints and Jacobian:\n"));
1239566063dSJacob Faibussowitsch       PetscCall(VecView(CI, PETSC_VIEWER_STDOUT_WORLD));
1249566063dSJacob Faibussowitsch       PetscCall(MatView(user.Ai, PETSC_VIEWER_STDOUT_WORLD));
1259566063dSJacob Faibussowitsch       PetscCall(VecDestroy(&CI));
126f560b561SHong Zhang       if (!user.noeqflag) {
1279566063dSJacob Faibussowitsch         PetscCall(FormEqualityJacobian(tao, user.x, user.Ae, user.Ae, (void *)&user));
1289566063dSJacob Faibussowitsch         PetscCall(MatCreateVecs(user.Ae, NULL, &CE));
1299566063dSJacob Faibussowitsch         PetscCall(FormEqualityConstraints(tao, user.x, CE, (void *)&user));
13063a3b9bcSJacob Faibussowitsch         PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nInitial equality constraints and Jacobian:\n"));
1319566063dSJacob Faibussowitsch         PetscCall(VecView(CE, PETSC_VIEWER_STDOUT_WORLD));
1329566063dSJacob Faibussowitsch         PetscCall(MatView(user.Ae, PETSC_VIEWER_STDOUT_WORLD));
1339566063dSJacob Faibussowitsch         PetscCall(VecDestroy(&CE));
134f560b561SHong Zhang       }
1359566063dSJacob Faibussowitsch       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n"));
1369566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_WORLD));
137f560b561SHong Zhang     }
138661095bbSAlp Dener   }
139aad13602SShrirang Abhyankar 
1409566063dSJacob Faibussowitsch   PetscCall(TaoSolve(tao));
1419566063dSJacob Faibussowitsch   PetscCall(TaoGetSolution(tao, &x));
1429566063dSJacob Faibussowitsch   PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD));
143aad13602SShrirang Abhyankar 
144aad13602SShrirang Abhyankar   /* Free objects */
1459566063dSJacob Faibussowitsch   PetscCall(DestroyProblem(&user));
1469566063dSJacob Faibussowitsch   PetscCall(TaoDestroy(&tao));
1479566063dSJacob Faibussowitsch   PetscCall(PetscFinalize());
148b122ec5aSJacob Faibussowitsch   return 0;
149aad13602SShrirang Abhyankar }
150aad13602SShrirang Abhyankar 
1519371c9d4SSatish Balay PetscErrorCode InitializeProblem(AppCtx *user) {
152aad13602SShrirang Abhyankar   PetscMPIInt size;
153aad13602SShrirang Abhyankar   PetscMPIInt rank;
154aad13602SShrirang Abhyankar   PetscInt    nloc, neloc, niloc;
155aad13602SShrirang Abhyankar 
156aad13602SShrirang Abhyankar   PetscFunctionBegin;
1579566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
1589566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
1598e58fa1dSresundermann   user->noeqflag = PETSC_FALSE;
1609566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL, NULL, "-no_eq", &user->noeqflag, NULL));
161f560b561SHong Zhang   user->initview = PETSC_FALSE;
1629566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL, NULL, "-init_view", &user->initview, NULL));
163628da978Sresundermann 
1648e58fa1dSresundermann   if (!user->noeqflag) {
165f560b561SHong Zhang     /* Tell user the correct solution, not an error checking */
1669566063dSJacob Faibussowitsch     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Solution should be f(1,1)=-2\n"));
1678e58fa1dSresundermann   }
168aad13602SShrirang Abhyankar 
1692d4ee042Sprj-   /* create vector x and set initial values */
170aad13602SShrirang Abhyankar   user->n = 2; /* global length */
171f560b561SHong Zhang   nloc    = (size == 1) ? user->n : 1;
1729566063dSJacob Faibussowitsch   PetscCall(VecCreate(PETSC_COMM_WORLD, &user->x));
1739566063dSJacob Faibussowitsch   PetscCall(VecSetSizes(user->x, nloc, user->n));
1749566063dSJacob Faibussowitsch   PetscCall(VecSetFromOptions(user->x));
1759566063dSJacob Faibussowitsch   PetscCall(VecSet(user->x, 0));
176aad13602SShrirang Abhyankar 
177aad13602SShrirang Abhyankar   /* create and set lower and upper bound vectors */
1789566063dSJacob Faibussowitsch   PetscCall(VecDuplicate(user->x, &user->xl));
1799566063dSJacob Faibussowitsch   PetscCall(VecDuplicate(user->x, &user->xu));
1809566063dSJacob Faibussowitsch   PetscCall(VecSet(user->xl, -1.0));
1819566063dSJacob Faibussowitsch   PetscCall(VecSet(user->xu, 2.0));
182aad13602SShrirang Abhyankar 
183aad13602SShrirang Abhyankar   /* create scater to zero */
1849566063dSJacob Faibussowitsch   PetscCall(VecScatterCreateToZero(user->x, &user->scat, &user->Xseq));
185aad13602SShrirang Abhyankar 
186aad13602SShrirang Abhyankar   user->ne = 1;
187628da978Sresundermann   user->ni = 2;
188aad13602SShrirang Abhyankar   neloc    = (rank == 0) ? user->ne : 0;
189f560b561SHong Zhang   niloc    = (size == 1) ? user->ni : 1;
190aad13602SShrirang Abhyankar 
1918e58fa1dSresundermann   if (!user->noeqflag) {
1929566063dSJacob Faibussowitsch     PetscCall(VecCreate(PETSC_COMM_WORLD, &user->ce)); /* a 1x1 vec for equality constraints */
1939566063dSJacob Faibussowitsch     PetscCall(VecSetSizes(user->ce, neloc, user->ne));
1949566063dSJacob Faibussowitsch     PetscCall(VecSetFromOptions(user->ce));
1959566063dSJacob Faibussowitsch     PetscCall(VecSetUp(user->ce));
1968e58fa1dSresundermann   }
197628da978Sresundermann 
1989566063dSJacob Faibussowitsch   PetscCall(VecCreate(PETSC_COMM_WORLD, &user->ci)); /* a 2x1 vec for inequality constraints */
1999566063dSJacob Faibussowitsch   PetscCall(VecSetSizes(user->ci, niloc, user->ni));
2009566063dSJacob Faibussowitsch   PetscCall(VecSetFromOptions(user->ci));
2019566063dSJacob Faibussowitsch   PetscCall(VecSetUp(user->ci));
202aad13602SShrirang Abhyankar 
203a5b23f4aSJose E. Roman   /* nexn & nixn matricies for equally and inequalty constraints */
2048e58fa1dSresundermann   if (!user->noeqflag) {
2059566063dSJacob Faibussowitsch     PetscCall(MatCreate(PETSC_COMM_WORLD, &user->Ae));
2069566063dSJacob Faibussowitsch     PetscCall(MatSetSizes(user->Ae, neloc, nloc, user->ne, user->n));
2079566063dSJacob Faibussowitsch     PetscCall(MatSetFromOptions(user->Ae));
2089566063dSJacob Faibussowitsch     PetscCall(MatSetUp(user->Ae));
2098e58fa1dSresundermann   }
2108e58fa1dSresundermann 
2119566063dSJacob Faibussowitsch   PetscCall(MatCreate(PETSC_COMM_WORLD, &user->Ai));
2129566063dSJacob Faibussowitsch   PetscCall(MatSetSizes(user->Ai, niloc, nloc, user->ni, user->n));
2139566063dSJacob Faibussowitsch   PetscCall(MatSetFromOptions(user->Ai));
2149566063dSJacob Faibussowitsch   PetscCall(MatSetUp(user->Ai));
215628da978Sresundermann 
2169566063dSJacob Faibussowitsch   PetscCall(MatCreate(PETSC_COMM_WORLD, &user->H));
2179566063dSJacob Faibussowitsch   PetscCall(MatSetSizes(user->H, nloc, nloc, user->n, user->n));
2189566063dSJacob Faibussowitsch   PetscCall(MatSetFromOptions(user->H));
2199566063dSJacob Faibussowitsch   PetscCall(MatSetUp(user->H));
220aad13602SShrirang Abhyankar   PetscFunctionReturn(0);
221aad13602SShrirang Abhyankar }
222aad13602SShrirang Abhyankar 
2239371c9d4SSatish Balay PetscErrorCode DestroyProblem(AppCtx *user) {
224aad13602SShrirang Abhyankar   PetscFunctionBegin;
225*48a46eb9SPierre Jolivet   if (!user->noeqflag) PetscCall(MatDestroy(&user->Ae));
2269566063dSJacob Faibussowitsch   PetscCall(MatDestroy(&user->Ai));
2279566063dSJacob Faibussowitsch   PetscCall(MatDestroy(&user->H));
228aad13602SShrirang Abhyankar 
2299566063dSJacob Faibussowitsch   PetscCall(VecDestroy(&user->x));
230*48a46eb9SPierre Jolivet   if (!user->noeqflag) PetscCall(VecDestroy(&user->ce));
2319566063dSJacob Faibussowitsch   PetscCall(VecDestroy(&user->ci));
2329566063dSJacob Faibussowitsch   PetscCall(VecDestroy(&user->xl));
2339566063dSJacob Faibussowitsch   PetscCall(VecDestroy(&user->xu));
2349566063dSJacob Faibussowitsch   PetscCall(VecDestroy(&user->Xseq));
2359566063dSJacob Faibussowitsch   PetscCall(VecScatterDestroy(&user->scat));
236aad13602SShrirang Abhyankar   PetscFunctionReturn(0);
237aad13602SShrirang Abhyankar }
238aad13602SShrirang Abhyankar 
239628da978Sresundermann /* Evaluate
240628da978Sresundermann    f(x) = (x0 - 2)^2 + (x1 - 2)^2 - 2*(x0 + x1)
241628da978Sresundermann    G = grad f = [2*(x0 - 2) - 2;
242628da978Sresundermann                  2*(x1 - 2) - 2]
243aad13602SShrirang Abhyankar */
2449371c9d4SSatish Balay PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ctx) {
245aad13602SShrirang Abhyankar   PetscScalar        g;
246aad13602SShrirang Abhyankar   const PetscScalar *x;
247aad13602SShrirang Abhyankar   MPI_Comm           comm;
248aad13602SShrirang Abhyankar   PetscMPIInt        rank;
249aad13602SShrirang Abhyankar   PetscReal          fin;
250aad13602SShrirang Abhyankar   AppCtx            *user = (AppCtx *)ctx;
251aad13602SShrirang Abhyankar   Vec                Xseq = user->Xseq;
252aad13602SShrirang Abhyankar   VecScatter         scat = user->scat;
253aad13602SShrirang Abhyankar 
254aad13602SShrirang Abhyankar   PetscFunctionBegin;
2559566063dSJacob Faibussowitsch   PetscCall(PetscObjectGetComm((PetscObject)tao, &comm));
2569566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Comm_rank(comm, &rank));
257aad13602SShrirang Abhyankar 
2589566063dSJacob Faibussowitsch   PetscCall(VecScatterBegin(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD));
2599566063dSJacob Faibussowitsch   PetscCall(VecScatterEnd(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD));
260aad13602SShrirang Abhyankar 
261aad13602SShrirang Abhyankar   fin = 0.0;
262dd400576SPatrick Sanan   if (rank == 0) {
2639566063dSJacob Faibussowitsch     PetscCall(VecGetArrayRead(Xseq, &x));
264aad13602SShrirang Abhyankar     fin = (x[0] - 2.0) * (x[0] - 2.0) + (x[1] - 2.0) * (x[1] - 2.0) - 2.0 * (x[0] + x[1]);
265aad13602SShrirang Abhyankar     g   = 2.0 * (x[0] - 2.0) - 2.0;
2669566063dSJacob Faibussowitsch     PetscCall(VecSetValue(G, 0, g, INSERT_VALUES));
267aad13602SShrirang Abhyankar     g = 2.0 * (x[1] - 2.0) - 2.0;
2689566063dSJacob Faibussowitsch     PetscCall(VecSetValue(G, 1, g, INSERT_VALUES));
2699566063dSJacob Faibussowitsch     PetscCall(VecRestoreArrayRead(Xseq, &x));
270aad13602SShrirang Abhyankar   }
2719566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Allreduce(&fin, f, 1, MPIU_REAL, MPIU_SUM, comm));
2729566063dSJacob Faibussowitsch   PetscCall(VecAssemblyBegin(G));
2739566063dSJacob Faibussowitsch   PetscCall(VecAssemblyEnd(G));
274aad13602SShrirang Abhyankar   PetscFunctionReturn(0);
275aad13602SShrirang Abhyankar }
276aad13602SShrirang Abhyankar 
277628da978Sresundermann /* Evaluate
278628da978Sresundermann    H = fxx + grad (grad g^T*DI) - grad (grad h^T*DE)]
279628da978Sresundermann      = [ 2*(1+de[0]-di[0]+di[1]), 0;
280628da978Sresundermann                    0,             2]
281628da978Sresundermann */
2829371c9d4SSatish Balay PetscErrorCode FormHessian(Tao tao, Vec x, Mat H, Mat Hpre, void *ctx) {
283aad13602SShrirang Abhyankar   AppCtx            *user = (AppCtx *)ctx;
284aad13602SShrirang Abhyankar   Vec                DE, DI;
285aad13602SShrirang Abhyankar   const PetscScalar *de, *di;
286aad13602SShrirang Abhyankar   PetscInt           zero = 0, one = 1;
287aad13602SShrirang Abhyankar   PetscScalar        two = 2.0;
288aad13602SShrirang Abhyankar   PetscScalar        val = 0.0;
289f560b561SHong Zhang   Vec                Deseq, Diseq;
290f560b561SHong Zhang   VecScatter         Descat, Discat;
291aad13602SShrirang Abhyankar   PetscMPIInt        rank;
292aad13602SShrirang Abhyankar   MPI_Comm           comm;
293aad13602SShrirang Abhyankar 
294aad13602SShrirang Abhyankar   PetscFunctionBegin;
2959566063dSJacob Faibussowitsch   PetscCall(TaoGetDualVariables(tao, &DE, &DI));
296aad13602SShrirang Abhyankar 
2979566063dSJacob Faibussowitsch   PetscCall(PetscObjectGetComm((PetscObject)tao, &comm));
2989566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Comm_rank(comm, &rank));
299aad13602SShrirang Abhyankar 
3008e58fa1dSresundermann   if (!user->noeqflag) {
3019566063dSJacob Faibussowitsch     PetscCall(VecScatterCreateToZero(DE, &Descat, &Deseq));
3029566063dSJacob Faibussowitsch     PetscCall(VecScatterBegin(Descat, DE, Deseq, INSERT_VALUES, SCATTER_FORWARD));
3039566063dSJacob Faibussowitsch     PetscCall(VecScatterEnd(Descat, DE, Deseq, INSERT_VALUES, SCATTER_FORWARD));
3048e58fa1dSresundermann   }
3059566063dSJacob Faibussowitsch   PetscCall(VecScatterCreateToZero(DI, &Discat, &Diseq));
3069566063dSJacob Faibussowitsch   PetscCall(VecScatterBegin(Discat, DI, Diseq, INSERT_VALUES, SCATTER_FORWARD));
3079566063dSJacob Faibussowitsch   PetscCall(VecScatterEnd(Discat, DI, Diseq, INSERT_VALUES, SCATTER_FORWARD));
308aad13602SShrirang Abhyankar 
309dd400576SPatrick Sanan   if (rank == 0) {
3109371c9d4SSatish Balay     if (!user->noeqflag) { PetscCall(VecGetArrayRead(Deseq, &de)); /* places equality constraint dual into array */ }
3119566063dSJacob Faibussowitsch     PetscCall(VecGetArrayRead(Diseq, &di)); /* places inequality constraint dual into array */
312628da978Sresundermann 
3138e58fa1dSresundermann     if (!user->noeqflag) {
314628da978Sresundermann       val = 2.0 * (1 + de[0] - di[0] + di[1]);
3159566063dSJacob Faibussowitsch       PetscCall(VecRestoreArrayRead(Deseq, &de));
3169566063dSJacob Faibussowitsch       PetscCall(VecRestoreArrayRead(Diseq, &di));
3178e58fa1dSresundermann     } else {
318628da978Sresundermann       val = 2.0 * (1 - di[0] + di[1]);
3198e58fa1dSresundermann     }
3209566063dSJacob Faibussowitsch     PetscCall(VecRestoreArrayRead(Diseq, &di));
3219566063dSJacob Faibussowitsch     PetscCall(MatSetValues(H, 1, &zero, 1, &zero, &val, INSERT_VALUES));
3229566063dSJacob Faibussowitsch     PetscCall(MatSetValues(H, 1, &one, 1, &one, &two, INSERT_VALUES));
323aad13602SShrirang Abhyankar   }
3249566063dSJacob Faibussowitsch   PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY));
3259566063dSJacob Faibussowitsch   PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY));
3268e58fa1dSresundermann   if (!user->noeqflag) {
3279566063dSJacob Faibussowitsch     PetscCall(VecScatterDestroy(&Descat));
3289566063dSJacob Faibussowitsch     PetscCall(VecDestroy(&Deseq));
3298e58fa1dSresundermann   }
3309566063dSJacob Faibussowitsch   PetscCall(VecScatterDestroy(&Discat));
3319566063dSJacob Faibussowitsch   PetscCall(VecDestroy(&Diseq));
332aad13602SShrirang Abhyankar   PetscFunctionReturn(0);
333aad13602SShrirang Abhyankar }
334aad13602SShrirang Abhyankar 
335628da978Sresundermann /* Evaluate
336628da978Sresundermann    h = [ x0^2 - x1;
337628da978Sresundermann          1 -(x0^2 - x1)]
338628da978Sresundermann */
3399371c9d4SSatish Balay PetscErrorCode FormInequalityConstraints(Tao tao, Vec X, Vec CI, void *ctx) {
340aad13602SShrirang Abhyankar   const PetscScalar *x;
341aad13602SShrirang Abhyankar   PetscScalar        ci;
342aad13602SShrirang Abhyankar   MPI_Comm           comm;
343aad13602SShrirang Abhyankar   PetscMPIInt        rank;
344aad13602SShrirang Abhyankar   AppCtx            *user = (AppCtx *)ctx;
345aad13602SShrirang Abhyankar   Vec                Xseq = user->Xseq;
346aad13602SShrirang Abhyankar   VecScatter         scat = user->scat;
347aad13602SShrirang Abhyankar 
348aad13602SShrirang Abhyankar   PetscFunctionBegin;
3499566063dSJacob Faibussowitsch   PetscCall(PetscObjectGetComm((PetscObject)tao, &comm));
3509566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Comm_rank(comm, &rank));
351aad13602SShrirang Abhyankar 
3529566063dSJacob Faibussowitsch   PetscCall(VecScatterBegin(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD));
3539566063dSJacob Faibussowitsch   PetscCall(VecScatterEnd(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD));
354aad13602SShrirang Abhyankar 
355dd400576SPatrick Sanan   if (rank == 0) {
3569566063dSJacob Faibussowitsch     PetscCall(VecGetArrayRead(Xseq, &x));
357aad13602SShrirang Abhyankar     ci = x[0] * x[0] - x[1];
3589566063dSJacob Faibussowitsch     PetscCall(VecSetValue(CI, 0, ci, INSERT_VALUES));
359aad13602SShrirang Abhyankar     ci = -x[0] * x[0] + x[1] + 1.0;
3609566063dSJacob Faibussowitsch     PetscCall(VecSetValue(CI, 1, ci, INSERT_VALUES));
3619566063dSJacob Faibussowitsch     PetscCall(VecRestoreArrayRead(Xseq, &x));
362aad13602SShrirang Abhyankar   }
3639566063dSJacob Faibussowitsch   PetscCall(VecAssemblyBegin(CI));
3649566063dSJacob Faibussowitsch   PetscCall(VecAssemblyEnd(CI));
365aad13602SShrirang Abhyankar   PetscFunctionReturn(0);
366aad13602SShrirang Abhyankar }
367aad13602SShrirang Abhyankar 
368628da978Sresundermann /* Evaluate
369628da978Sresundermann    g = [ x0^2 + x1 - 2]
370628da978Sresundermann */
3719371c9d4SSatish Balay PetscErrorCode FormEqualityConstraints(Tao tao, Vec X, Vec CE, void *ctx) {
372aad13602SShrirang Abhyankar   const PetscScalar *x;
373aad13602SShrirang Abhyankar   PetscScalar        ce;
374aad13602SShrirang Abhyankar   MPI_Comm           comm;
375aad13602SShrirang Abhyankar   PetscMPIInt        rank;
376aad13602SShrirang Abhyankar   AppCtx            *user = (AppCtx *)ctx;
377aad13602SShrirang Abhyankar   Vec                Xseq = user->Xseq;
378aad13602SShrirang Abhyankar   VecScatter         scat = user->scat;
379aad13602SShrirang Abhyankar 
380aad13602SShrirang Abhyankar   PetscFunctionBegin;
3819566063dSJacob Faibussowitsch   PetscCall(PetscObjectGetComm((PetscObject)tao, &comm));
3829566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Comm_rank(comm, &rank));
383aad13602SShrirang Abhyankar 
3849566063dSJacob Faibussowitsch   PetscCall(VecScatterBegin(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD));
3859566063dSJacob Faibussowitsch   PetscCall(VecScatterEnd(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD));
386aad13602SShrirang Abhyankar 
387dd400576SPatrick Sanan   if (rank == 0) {
3889566063dSJacob Faibussowitsch     PetscCall(VecGetArrayRead(Xseq, &x));
389aad13602SShrirang Abhyankar     ce = x[0] * x[0] + x[1] - 2.0;
3909566063dSJacob Faibussowitsch     PetscCall(VecSetValue(CE, 0, ce, INSERT_VALUES));
3919566063dSJacob Faibussowitsch     PetscCall(VecRestoreArrayRead(Xseq, &x));
392aad13602SShrirang Abhyankar   }
3939566063dSJacob Faibussowitsch   PetscCall(VecAssemblyBegin(CE));
3949566063dSJacob Faibussowitsch   PetscCall(VecAssemblyEnd(CE));
395aad13602SShrirang Abhyankar   PetscFunctionReturn(0);
396aad13602SShrirang Abhyankar }
397aad13602SShrirang Abhyankar 
398628da978Sresundermann /*
399628da978Sresundermann   grad h = [  2*x0, -1;
400628da978Sresundermann              -2*x0,  1]
401628da978Sresundermann */
4029371c9d4SSatish Balay PetscErrorCode FormInequalityJacobian(Tao tao, Vec X, Mat JI, Mat JIpre, void *ctx) {
403aad13602SShrirang Abhyankar   AppCtx            *user = (AppCtx *)ctx;
404f560b561SHong Zhang   PetscInt           zero = 0, one = 1, cols[2];
405aad13602SShrirang Abhyankar   PetscScalar        vals[2];
406aad13602SShrirang Abhyankar   const PetscScalar *x;
407aad13602SShrirang Abhyankar   Vec                Xseq = user->Xseq;
408aad13602SShrirang Abhyankar   VecScatter         scat = user->scat;
409aad13602SShrirang Abhyankar   MPI_Comm           comm;
410aad13602SShrirang Abhyankar   PetscMPIInt        rank;
411aad13602SShrirang Abhyankar 
412aad13602SShrirang Abhyankar   PetscFunctionBegin;
4139566063dSJacob Faibussowitsch   PetscCall(PetscObjectGetComm((PetscObject)tao, &comm));
4149566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Comm_rank(comm, &rank));
4159566063dSJacob Faibussowitsch   PetscCall(VecScatterBegin(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD));
4169566063dSJacob Faibussowitsch   PetscCall(VecScatterEnd(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD));
417aad13602SShrirang Abhyankar 
4189566063dSJacob Faibussowitsch   PetscCall(VecGetArrayRead(Xseq, &x));
419c5853193SPierre Jolivet   if (rank == 0) {
4209371c9d4SSatish Balay     cols[0] = 0;
4219371c9d4SSatish Balay     cols[1] = 1;
4229371c9d4SSatish Balay     vals[0] = 2 * x[0];
4239371c9d4SSatish Balay     vals[1] = -1.0;
4249566063dSJacob Faibussowitsch     PetscCall(MatSetValues(JI, 1, &zero, 2, cols, vals, INSERT_VALUES));
4259371c9d4SSatish Balay     vals[0] = -2 * x[0];
4269371c9d4SSatish Balay     vals[1] = 1.0;
4279566063dSJacob Faibussowitsch     PetscCall(MatSetValues(JI, 1, &one, 2, cols, vals, INSERT_VALUES));
428aad13602SShrirang Abhyankar   }
4299566063dSJacob Faibussowitsch   PetscCall(VecRestoreArrayRead(Xseq, &x));
4309566063dSJacob Faibussowitsch   PetscCall(MatAssemblyBegin(JI, MAT_FINAL_ASSEMBLY));
4319566063dSJacob Faibussowitsch   PetscCall(MatAssemblyEnd(JI, MAT_FINAL_ASSEMBLY));
432aad13602SShrirang Abhyankar   PetscFunctionReturn(0);
433aad13602SShrirang Abhyankar }
434aad13602SShrirang Abhyankar 
435628da978Sresundermann /*
436628da978Sresundermann   grad g = [2*x0
437628da978Sresundermann              1.0 ]
438628da978Sresundermann */
4399371c9d4SSatish Balay PetscErrorCode FormEqualityJacobian(Tao tao, Vec X, Mat JE, Mat JEpre, void *ctx) {
440f560b561SHong Zhang   PetscInt           zero = 0, cols[2];
441aad13602SShrirang Abhyankar   PetscScalar        vals[2];
442aad13602SShrirang Abhyankar   const PetscScalar *x;
443aad13602SShrirang Abhyankar   PetscMPIInt        rank;
444aad13602SShrirang Abhyankar   MPI_Comm           comm;
445aad13602SShrirang Abhyankar 
446aad13602SShrirang Abhyankar   PetscFunctionBegin;
4479566063dSJacob Faibussowitsch   PetscCall(PetscObjectGetComm((PetscObject)tao, &comm));
4489566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Comm_rank(comm, &rank));
449aad13602SShrirang Abhyankar 
450dd400576SPatrick Sanan   if (rank == 0) {
4519566063dSJacob Faibussowitsch     PetscCall(VecGetArrayRead(X, &x));
4529371c9d4SSatish Balay     cols[0] = 0;
4539371c9d4SSatish Balay     cols[1] = 1;
4549371c9d4SSatish Balay     vals[0] = 2 * x[0];
4559371c9d4SSatish Balay     vals[1] = 1.0;
4569566063dSJacob Faibussowitsch     PetscCall(MatSetValues(JE, 1, &zero, 2, cols, vals, INSERT_VALUES));
4579566063dSJacob Faibussowitsch     PetscCall(VecRestoreArrayRead(X, &x));
458aad13602SShrirang Abhyankar   }
4599566063dSJacob Faibussowitsch   PetscCall(MatAssemblyBegin(JE, MAT_FINAL_ASSEMBLY));
4609566063dSJacob Faibussowitsch   PetscCall(MatAssemblyEnd(JE, MAT_FINAL_ASSEMBLY));
461aad13602SShrirang Abhyankar   PetscFunctionReturn(0);
462aad13602SShrirang Abhyankar }
463aad13602SShrirang Abhyankar 
464aad13602SShrirang Abhyankar /*TEST
465aad13602SShrirang Abhyankar 
466aad13602SShrirang Abhyankar    build:
467f560b561SHong Zhang       requires: !complex !defined(PETSC_USE_CXX)
468aad13602SShrirang Abhyankar 
469aad13602SShrirang Abhyankar    test:
47009ee8bb0SRylee Sundermann       args: -tao_converged_reason -tao_pdipm_kkt_shift_pd
471f560b561SHong Zhang       requires: mumps
472aad13602SShrirang Abhyankar 
473aad13602SShrirang Abhyankar    test:
474aad13602SShrirang Abhyankar       suffix: 2
475aad13602SShrirang Abhyankar       nsize: 2
47609ee8bb0SRylee Sundermann       args: -tao_converged_reason -tao_pdipm_kkt_shift_pd
477f560b561SHong Zhang       requires: mumps
478aad13602SShrirang Abhyankar 
4798e58fa1dSresundermann    test:
4808e58fa1dSresundermann       suffix: 3
4818e58fa1dSresundermann       args: -tao_converged_reason -no_eq
482f560b561SHong Zhang       requires: mumps
4838e58fa1dSresundermann 
4848e58fa1dSresundermann    test:
4858e58fa1dSresundermann       suffix: 4
4868e58fa1dSresundermann       nsize: 2
4878e58fa1dSresundermann       args: -tao_converged_reason -no_eq
488f560b561SHong Zhang       requires: mumps
4898e58fa1dSresundermann 
490661095bbSAlp Dener    test:
491661095bbSAlp Dener       suffix: 5
492661095bbSAlp Dener       args: -tao_cmonitor -tao_type almm
493f560b561SHong Zhang       requires: mumps
494661095bbSAlp Dener 
495661095bbSAlp Dener    test:
496661095bbSAlp Dener       suffix: 6
497661095bbSAlp Dener       args: -tao_cmonitor -tao_type almm -tao_almm_type phr
498f560b561SHong Zhang       requires: mumps
499661095bbSAlp Dener 
500661095bbSAlp Dener    test:
501661095bbSAlp Dener       suffix: 7
502661095bbSAlp Dener       nsize: 2
503f560b561SHong Zhang       requires: mumps
504661095bbSAlp Dener       args: -tao_cmonitor -tao_type almm
505661095bbSAlp Dener 
506661095bbSAlp Dener    test:
507661095bbSAlp Dener       suffix: 8
508661095bbSAlp Dener       nsize: 2
509f560b561SHong Zhang       requires: cuda mumps
510661095bbSAlp Dener       args: -tao_cmonitor -tao_type almm -vec_type cuda -mat_type aijcusparse
511661095bbSAlp Dener 
512661095bbSAlp Dener    test:
513661095bbSAlp Dener       suffix: 9
514661095bbSAlp Dener       nsize: 2
515661095bbSAlp Dener       args: -tao_cmonitor -tao_type almm -no_eq
516f560b561SHong Zhang       requires: mumps
517661095bbSAlp Dener 
518661095bbSAlp Dener    test:
519661095bbSAlp Dener       suffix: 10
520661095bbSAlp Dener       nsize: 2
521661095bbSAlp Dener       args: -tao_cmonitor -tao_type almm -tao_almm_type phr -no_eq
522f560b561SHong Zhang       requires: mumps
523661095bbSAlp Dener 
524aad13602SShrirang Abhyankar TEST*/
525