xref: /petsc/src/tao/constrained/tutorials/ex1.c (revision b9e5a4c7d30bdcab0a753290c79ea6b79adb512b)
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 
20da81f932SPierre Jolivet static char help[] = "Solves constrained optimization 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 
55d71ae5a4SJacob Faibussowitsch PetscErrorCode main(int argc, char **argv)
56d71ae5a4SJacob Faibussowitsch {
57aad13602SShrirang Abhyankar   Tao         tao;
58aad13602SShrirang Abhyankar   KSP         ksp;
59aad13602SShrirang Abhyankar   PC          pc;
60aad13602SShrirang Abhyankar   AppCtx      user; /* application context */
61f560b561SHong Zhang   Vec         x, G, CI, CE;
62f560b561SHong Zhang   PetscMPIInt size;
63661095bbSAlp Dener   TaoType     type;
64f560b561SHong Zhang   PetscReal   f;
659a09327dSAlp Dener   PetscBool   pdipm, mumps;
66aad13602SShrirang Abhyankar 
67327415f7SBarry Smith   PetscFunctionBeginUser;
689566063dSJacob Faibussowitsch   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
699566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
709a09327dSAlp Dener   PetscCheck(size <= 2, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "More than 2 processors detected. Example written to use max of 2 processors.");
71aad13602SShrirang Abhyankar 
729566063dSJacob Faibussowitsch   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "---- Constrained Problem -----\n"));
739566063dSJacob Faibussowitsch   PetscCall(InitializeProblem(&user)); /* sets up problem, function below */
749a09327dSAlp Dener 
759566063dSJacob Faibussowitsch   PetscCall(TaoCreate(PETSC_COMM_WORLD, &tao));
769a09327dSAlp Dener   PetscCall(TaoSetType(tao, TAOALMM));
779a09327dSAlp Dener   PetscCall(TaoSetSolution(tao, user.x));
789a09327dSAlp Dener   PetscCall(TaoSetVariableBounds(tao, user.xl, user.xu));
799566063dSJacob Faibussowitsch   PetscCall(TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, (void *)&user));
809a09327dSAlp Dener   PetscCall(TaoSetTolerances(tao, 1.e-4, 0.0, 0.0));
819a09327dSAlp Dener   PetscCall(TaoSetConstraintTolerances(tao, 1.e-4, 0.0));
829a09327dSAlp Dener   PetscCall(TaoSetMaximumFunctionEvaluations(tao, 1e6));
839566063dSJacob Faibussowitsch   PetscCall(TaoSetFromOptions(tao));
849a09327dSAlp Dener 
859a09327dSAlp Dener   if (!user.noeqflag) {
869a09327dSAlp Dener     PetscCall(TaoSetEqualityConstraintsRoutine(tao, user.ce, FormEqualityConstraints, (void *)&user));
879a09327dSAlp Dener     PetscCall(TaoSetJacobianEqualityRoutine(tao, user.Ae, user.Ae, FormEqualityJacobian, (void *)&user));
889a09327dSAlp Dener   }
899a09327dSAlp Dener   PetscCall(TaoSetInequalityConstraintsRoutine(tao, user.ci, FormInequalityConstraints, (void *)&user));
909a09327dSAlp Dener   PetscCall(TaoSetJacobianInequalityRoutine(tao, user.Ai, user.Ai, FormInequalityJacobian, (void *)&user));
919a09327dSAlp Dener 
929566063dSJacob Faibussowitsch   PetscCall(TaoGetType(tao, &type));
939566063dSJacob Faibussowitsch   PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOPDIPM, &pdipm));
94661095bbSAlp Dener   if (pdipm) {
959a09327dSAlp Dener     /*
969a09327dSAlp Dener       PDIPM produces an indefinite KKT matrix with some zeros along the diagonal
979a09327dSAlp Dener       Inverting this indefinite matrix requires PETSc to be configured with MUMPS
989a09327dSAlp Dener     */
999a09327dSAlp Dener     PetscCall(PetscHasExternalPackage("mumps", &mumps));
1009a09327dSAlp Dener     PetscCheck(mumps, PetscObjectComm((PetscObject)tao), PETSC_ERR_SUP, "TAOPDIPM requires PETSc to be configured with MUMPS (--download-mumps)");
1019a09327dSAlp Dener     PetscCall(TaoGetKSP(tao, &ksp));
1029a09327dSAlp Dener     PetscCall(KSPGetPC(ksp, &pc));
1039a09327dSAlp Dener     PetscCall(PCFactorSetMatSolverType(pc, MATSOLVERMUMPS));
1049a09327dSAlp Dener     PetscCall(KSPSetType(ksp, KSPPREONLY));
1059a09327dSAlp Dener     PetscCall(KSPSetFromOptions(ksp));
1069566063dSJacob Faibussowitsch     PetscCall(TaoSetHessian(tao, user.H, user.H, FormHessian, (void *)&user));
1079a09327dSAlp Dener   }
1089a09327dSAlp Dener 
1099a09327dSAlp Dener   /* Print out an initial view of the problem */
110f560b561SHong Zhang   if (user.initview) {
1119566063dSJacob Faibussowitsch     PetscCall(TaoSetUp(tao));
1129566063dSJacob Faibussowitsch     PetscCall(VecDuplicate(user.x, &G));
1139566063dSJacob Faibussowitsch     PetscCall(FormFunctionGradient(tao, user.x, &f, G, (void *)&user));
1149566063dSJacob Faibussowitsch     PetscCall(FormHessian(tao, user.x, user.H, user.H, (void *)&user));
1159566063dSJacob Faibussowitsch     PetscCall(PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_WORLD));
11663a3b9bcSJacob Faibussowitsch     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nInitial point X:\n"));
1179566063dSJacob Faibussowitsch     PetscCall(VecView(user.x, PETSC_VIEWER_STDOUT_WORLD));
11863a3b9bcSJacob Faibussowitsch     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nInitial objective f(x) = %g\n", (double)f));
11963a3b9bcSJacob Faibussowitsch     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nInitial gradient and Hessian:\n"));
1209566063dSJacob Faibussowitsch     PetscCall(VecView(G, PETSC_VIEWER_STDOUT_WORLD));
1219566063dSJacob Faibussowitsch     PetscCall(MatView(user.H, PETSC_VIEWER_STDOUT_WORLD));
1229566063dSJacob Faibussowitsch     PetscCall(VecDestroy(&G));
1239566063dSJacob Faibussowitsch     PetscCall(FormInequalityJacobian(tao, user.x, user.Ai, user.Ai, (void *)&user));
1249566063dSJacob Faibussowitsch     PetscCall(MatCreateVecs(user.Ai, NULL, &CI));
1259566063dSJacob Faibussowitsch     PetscCall(FormInequalityConstraints(tao, user.x, CI, (void *)&user));
12663a3b9bcSJacob Faibussowitsch     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nInitial inequality constraints and Jacobian:\n"));
1279566063dSJacob Faibussowitsch     PetscCall(VecView(CI, PETSC_VIEWER_STDOUT_WORLD));
1289566063dSJacob Faibussowitsch     PetscCall(MatView(user.Ai, PETSC_VIEWER_STDOUT_WORLD));
1299566063dSJacob Faibussowitsch     PetscCall(VecDestroy(&CI));
130f560b561SHong Zhang     if (!user.noeqflag) {
1319566063dSJacob Faibussowitsch       PetscCall(FormEqualityJacobian(tao, user.x, user.Ae, user.Ae, (void *)&user));
1329566063dSJacob Faibussowitsch       PetscCall(MatCreateVecs(user.Ae, NULL, &CE));
1339566063dSJacob Faibussowitsch       PetscCall(FormEqualityConstraints(tao, user.x, CE, (void *)&user));
13463a3b9bcSJacob Faibussowitsch       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nInitial equality constraints and Jacobian:\n"));
1359566063dSJacob Faibussowitsch       PetscCall(VecView(CE, PETSC_VIEWER_STDOUT_WORLD));
1369566063dSJacob Faibussowitsch       PetscCall(MatView(user.Ae, PETSC_VIEWER_STDOUT_WORLD));
1379566063dSJacob Faibussowitsch       PetscCall(VecDestroy(&CE));
138f560b561SHong Zhang     }
1399566063dSJacob Faibussowitsch     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n"));
1409566063dSJacob Faibussowitsch     PetscCall(PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_WORLD));
141f560b561SHong Zhang   }
142aad13602SShrirang Abhyankar 
1439566063dSJacob Faibussowitsch   PetscCall(TaoSolve(tao));
1449566063dSJacob Faibussowitsch   PetscCall(TaoGetSolution(tao, &x));
1459a09327dSAlp Dener   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Found solution:\n"));
1469566063dSJacob Faibussowitsch   PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD));
147aad13602SShrirang Abhyankar 
148aad13602SShrirang Abhyankar   /* Free objects */
1499566063dSJacob Faibussowitsch   PetscCall(DestroyProblem(&user));
1509566063dSJacob Faibussowitsch   PetscCall(TaoDestroy(&tao));
1519566063dSJacob Faibussowitsch   PetscCall(PetscFinalize());
152b122ec5aSJacob Faibussowitsch   return 0;
153aad13602SShrirang Abhyankar }
154aad13602SShrirang Abhyankar 
155d71ae5a4SJacob Faibussowitsch PetscErrorCode InitializeProblem(AppCtx *user)
156d71ae5a4SJacob Faibussowitsch {
157aad13602SShrirang Abhyankar   PetscMPIInt size;
158aad13602SShrirang Abhyankar   PetscMPIInt rank;
159aad13602SShrirang Abhyankar   PetscInt    nloc, neloc, niloc;
160aad13602SShrirang Abhyankar 
161aad13602SShrirang Abhyankar   PetscFunctionBegin;
1629566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
1639566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
1648e58fa1dSresundermann   user->noeqflag = PETSC_FALSE;
1659566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL, NULL, "-no_eq", &user->noeqflag, NULL));
166f560b561SHong Zhang   user->initview = PETSC_FALSE;
1679566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL, NULL, "-init_view", &user->initview, NULL));
168628da978Sresundermann 
169f560b561SHong Zhang   /* Tell user the correct solution, not an error checking */
1709a09327dSAlp Dener   if (!user->noeqflag) {
1719a09327dSAlp Dener     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Expected solution: f(1, 1) = -2\n"));
1729a09327dSAlp Dener   } else {
1739a09327dSAlp Dener     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Expected solution (-no_eq): f(1.73205, 2) = -7.3923\n"));
1748e58fa1dSresundermann   }
175aad13602SShrirang Abhyankar 
1762d4ee042Sprj-   /* create vector x and set initial values */
177aad13602SShrirang Abhyankar   user->n = 2; /* global length */
178f560b561SHong Zhang   nloc    = (size == 1) ? user->n : 1;
1799566063dSJacob Faibussowitsch   PetscCall(VecCreate(PETSC_COMM_WORLD, &user->x));
1809566063dSJacob Faibussowitsch   PetscCall(VecSetSizes(user->x, nloc, user->n));
1819566063dSJacob Faibussowitsch   PetscCall(VecSetFromOptions(user->x));
1829a09327dSAlp Dener   PetscCall(VecSet(user->x, 0.0));
183aad13602SShrirang Abhyankar 
184aad13602SShrirang Abhyankar   /* create and set lower and upper bound vectors */
1859566063dSJacob Faibussowitsch   PetscCall(VecDuplicate(user->x, &user->xl));
1869566063dSJacob Faibussowitsch   PetscCall(VecDuplicate(user->x, &user->xu));
1879566063dSJacob Faibussowitsch   PetscCall(VecSet(user->xl, -1.0));
1889566063dSJacob Faibussowitsch   PetscCall(VecSet(user->xu, 2.0));
189aad13602SShrirang Abhyankar 
190aad13602SShrirang Abhyankar   /* create scater to zero */
1919566063dSJacob Faibussowitsch   PetscCall(VecScatterCreateToZero(user->x, &user->scat, &user->Xseq));
192aad13602SShrirang Abhyankar 
193aad13602SShrirang Abhyankar   user->ne = 1;
194628da978Sresundermann   user->ni = 2;
195aad13602SShrirang Abhyankar   neloc    = (rank == 0) ? user->ne : 0;
196f560b561SHong Zhang   niloc    = (size == 1) ? user->ni : 1;
197aad13602SShrirang Abhyankar 
1988e58fa1dSresundermann   if (!user->noeqflag) {
1999566063dSJacob Faibussowitsch     PetscCall(VecCreate(PETSC_COMM_WORLD, &user->ce)); /* a 1x1 vec for equality constraints */
2009566063dSJacob Faibussowitsch     PetscCall(VecSetSizes(user->ce, neloc, user->ne));
2019566063dSJacob Faibussowitsch     PetscCall(VecSetFromOptions(user->ce));
2029566063dSJacob Faibussowitsch     PetscCall(VecSetUp(user->ce));
2038e58fa1dSresundermann   }
204628da978Sresundermann 
2059566063dSJacob Faibussowitsch   PetscCall(VecCreate(PETSC_COMM_WORLD, &user->ci)); /* a 2x1 vec for inequality constraints */
2069566063dSJacob Faibussowitsch   PetscCall(VecSetSizes(user->ci, niloc, user->ni));
2079566063dSJacob Faibussowitsch   PetscCall(VecSetFromOptions(user->ci));
2089566063dSJacob Faibussowitsch   PetscCall(VecSetUp(user->ci));
209aad13602SShrirang Abhyankar 
210a5b23f4aSJose E. Roman   /* nexn & nixn matricies for equally and inequalty constraints */
2118e58fa1dSresundermann   if (!user->noeqflag) {
2129566063dSJacob Faibussowitsch     PetscCall(MatCreate(PETSC_COMM_WORLD, &user->Ae));
2139566063dSJacob Faibussowitsch     PetscCall(MatSetSizes(user->Ae, neloc, nloc, user->ne, user->n));
2149566063dSJacob Faibussowitsch     PetscCall(MatSetFromOptions(user->Ae));
2159566063dSJacob Faibussowitsch     PetscCall(MatSetUp(user->Ae));
2168e58fa1dSresundermann   }
2178e58fa1dSresundermann 
2189566063dSJacob Faibussowitsch   PetscCall(MatCreate(PETSC_COMM_WORLD, &user->Ai));
2199566063dSJacob Faibussowitsch   PetscCall(MatSetSizes(user->Ai, niloc, nloc, user->ni, user->n));
2209566063dSJacob Faibussowitsch   PetscCall(MatSetFromOptions(user->Ai));
2219566063dSJacob Faibussowitsch   PetscCall(MatSetUp(user->Ai));
222628da978Sresundermann 
2239566063dSJacob Faibussowitsch   PetscCall(MatCreate(PETSC_COMM_WORLD, &user->H));
2249566063dSJacob Faibussowitsch   PetscCall(MatSetSizes(user->H, nloc, nloc, user->n, user->n));
2259566063dSJacob Faibussowitsch   PetscCall(MatSetFromOptions(user->H));
2269566063dSJacob Faibussowitsch   PetscCall(MatSetUp(user->H));
227aad13602SShrirang Abhyankar   PetscFunctionReturn(0);
228aad13602SShrirang Abhyankar }
229aad13602SShrirang Abhyankar 
230d71ae5a4SJacob Faibussowitsch PetscErrorCode DestroyProblem(AppCtx *user)
231d71ae5a4SJacob Faibussowitsch {
232aad13602SShrirang Abhyankar   PetscFunctionBegin;
23348a46eb9SPierre Jolivet   if (!user->noeqflag) PetscCall(MatDestroy(&user->Ae));
2349566063dSJacob Faibussowitsch   PetscCall(MatDestroy(&user->Ai));
2359566063dSJacob Faibussowitsch   PetscCall(MatDestroy(&user->H));
236aad13602SShrirang Abhyankar 
2379566063dSJacob Faibussowitsch   PetscCall(VecDestroy(&user->x));
23848a46eb9SPierre Jolivet   if (!user->noeqflag) PetscCall(VecDestroy(&user->ce));
2399566063dSJacob Faibussowitsch   PetscCall(VecDestroy(&user->ci));
2409566063dSJacob Faibussowitsch   PetscCall(VecDestroy(&user->xl));
2419566063dSJacob Faibussowitsch   PetscCall(VecDestroy(&user->xu));
2429566063dSJacob Faibussowitsch   PetscCall(VecDestroy(&user->Xseq));
2439566063dSJacob Faibussowitsch   PetscCall(VecScatterDestroy(&user->scat));
244aad13602SShrirang Abhyankar   PetscFunctionReturn(0);
245aad13602SShrirang Abhyankar }
246aad13602SShrirang Abhyankar 
247628da978Sresundermann /* Evaluate
248628da978Sresundermann    f(x) = (x0 - 2)^2 + (x1 - 2)^2 - 2*(x0 + x1)
249628da978Sresundermann    G = grad f = [2*(x0 - 2) - 2;
250628da978Sresundermann                  2*(x1 - 2) - 2]
251aad13602SShrirang Abhyankar */
252d71ae5a4SJacob Faibussowitsch PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ctx)
253d71ae5a4SJacob Faibussowitsch {
254aad13602SShrirang Abhyankar   PetscScalar        g;
255aad13602SShrirang Abhyankar   const PetscScalar *x;
256aad13602SShrirang Abhyankar   MPI_Comm           comm;
257aad13602SShrirang Abhyankar   PetscMPIInt        rank;
258aad13602SShrirang Abhyankar   PetscReal          fin;
259aad13602SShrirang Abhyankar   AppCtx            *user = (AppCtx *)ctx;
260aad13602SShrirang Abhyankar   Vec                Xseq = user->Xseq;
261aad13602SShrirang Abhyankar   VecScatter         scat = user->scat;
262aad13602SShrirang Abhyankar 
263aad13602SShrirang Abhyankar   PetscFunctionBegin;
2649566063dSJacob Faibussowitsch   PetscCall(PetscObjectGetComm((PetscObject)tao, &comm));
2659566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Comm_rank(comm, &rank));
266aad13602SShrirang Abhyankar 
2679566063dSJacob Faibussowitsch   PetscCall(VecScatterBegin(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD));
2689566063dSJacob Faibussowitsch   PetscCall(VecScatterEnd(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD));
269aad13602SShrirang Abhyankar 
270aad13602SShrirang Abhyankar   fin = 0.0;
271dd400576SPatrick Sanan   if (rank == 0) {
2729566063dSJacob Faibussowitsch     PetscCall(VecGetArrayRead(Xseq, &x));
273aad13602SShrirang Abhyankar     fin = (x[0] - 2.0) * (x[0] - 2.0) + (x[1] - 2.0) * (x[1] - 2.0) - 2.0 * (x[0] + x[1]);
274aad13602SShrirang Abhyankar     g   = 2.0 * (x[0] - 2.0) - 2.0;
2759566063dSJacob Faibussowitsch     PetscCall(VecSetValue(G, 0, g, INSERT_VALUES));
276aad13602SShrirang Abhyankar     g = 2.0 * (x[1] - 2.0) - 2.0;
2779566063dSJacob Faibussowitsch     PetscCall(VecSetValue(G, 1, g, INSERT_VALUES));
2789566063dSJacob Faibussowitsch     PetscCall(VecRestoreArrayRead(Xseq, &x));
279aad13602SShrirang Abhyankar   }
2809566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Allreduce(&fin, f, 1, MPIU_REAL, MPIU_SUM, comm));
2819566063dSJacob Faibussowitsch   PetscCall(VecAssemblyBegin(G));
2829566063dSJacob Faibussowitsch   PetscCall(VecAssemblyEnd(G));
283aad13602SShrirang Abhyankar   PetscFunctionReturn(0);
284aad13602SShrirang Abhyankar }
285aad13602SShrirang Abhyankar 
286628da978Sresundermann /* Evaluate
287628da978Sresundermann    H = fxx + grad (grad g^T*DI) - grad (grad h^T*DE)]
288628da978Sresundermann      = [ 2*(1+de[0]-di[0]+di[1]), 0;
289628da978Sresundermann                    0,             2]
290628da978Sresundermann */
291d71ae5a4SJacob Faibussowitsch PetscErrorCode FormHessian(Tao tao, Vec x, Mat H, Mat Hpre, void *ctx)
292d71ae5a4SJacob Faibussowitsch {
293aad13602SShrirang Abhyankar   AppCtx            *user = (AppCtx *)ctx;
294aad13602SShrirang Abhyankar   Vec                DE, DI;
295aad13602SShrirang Abhyankar   const PetscScalar *de, *di;
296aad13602SShrirang Abhyankar   PetscInt           zero = 0, one = 1;
297aad13602SShrirang Abhyankar   PetscScalar        two = 2.0;
298aad13602SShrirang Abhyankar   PetscScalar        val = 0.0;
299f560b561SHong Zhang   Vec                Deseq, Diseq;
300f560b561SHong Zhang   VecScatter         Descat, Discat;
301aad13602SShrirang Abhyankar   PetscMPIInt        rank;
302aad13602SShrirang Abhyankar   MPI_Comm           comm;
303aad13602SShrirang Abhyankar 
304aad13602SShrirang Abhyankar   PetscFunctionBegin;
3059566063dSJacob Faibussowitsch   PetscCall(TaoGetDualVariables(tao, &DE, &DI));
306aad13602SShrirang Abhyankar 
3079566063dSJacob Faibussowitsch   PetscCall(PetscObjectGetComm((PetscObject)tao, &comm));
3089566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Comm_rank(comm, &rank));
309aad13602SShrirang Abhyankar 
3108e58fa1dSresundermann   if (!user->noeqflag) {
3119566063dSJacob Faibussowitsch     PetscCall(VecScatterCreateToZero(DE, &Descat, &Deseq));
3129566063dSJacob Faibussowitsch     PetscCall(VecScatterBegin(Descat, DE, Deseq, INSERT_VALUES, SCATTER_FORWARD));
3139566063dSJacob Faibussowitsch     PetscCall(VecScatterEnd(Descat, DE, Deseq, INSERT_VALUES, SCATTER_FORWARD));
3148e58fa1dSresundermann   }
3159566063dSJacob Faibussowitsch   PetscCall(VecScatterCreateToZero(DI, &Discat, &Diseq));
3169566063dSJacob Faibussowitsch   PetscCall(VecScatterBegin(Discat, DI, Diseq, INSERT_VALUES, SCATTER_FORWARD));
3179566063dSJacob Faibussowitsch   PetscCall(VecScatterEnd(Discat, DI, Diseq, INSERT_VALUES, SCATTER_FORWARD));
318aad13602SShrirang Abhyankar 
319dd400576SPatrick Sanan   if (rank == 0) {
3209371c9d4SSatish Balay     if (!user->noeqflag) { PetscCall(VecGetArrayRead(Deseq, &de)); /* places equality constraint dual into array */ }
3219566063dSJacob Faibussowitsch     PetscCall(VecGetArrayRead(Diseq, &di)); /* places inequality constraint dual into array */
322628da978Sresundermann 
3238e58fa1dSresundermann     if (!user->noeqflag) {
324628da978Sresundermann       val = 2.0 * (1 + de[0] - di[0] + di[1]);
3259566063dSJacob Faibussowitsch       PetscCall(VecRestoreArrayRead(Deseq, &de));
3269566063dSJacob Faibussowitsch       PetscCall(VecRestoreArrayRead(Diseq, &di));
3278e58fa1dSresundermann     } else {
328628da978Sresundermann       val = 2.0 * (1 - di[0] + di[1]);
3298e58fa1dSresundermann     }
3309566063dSJacob Faibussowitsch     PetscCall(VecRestoreArrayRead(Diseq, &di));
3319566063dSJacob Faibussowitsch     PetscCall(MatSetValues(H, 1, &zero, 1, &zero, &val, INSERT_VALUES));
3329566063dSJacob Faibussowitsch     PetscCall(MatSetValues(H, 1, &one, 1, &one, &two, INSERT_VALUES));
333aad13602SShrirang Abhyankar   }
3349566063dSJacob Faibussowitsch   PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY));
3359566063dSJacob Faibussowitsch   PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY));
3368e58fa1dSresundermann   if (!user->noeqflag) {
3379566063dSJacob Faibussowitsch     PetscCall(VecScatterDestroy(&Descat));
3389566063dSJacob Faibussowitsch     PetscCall(VecDestroy(&Deseq));
3398e58fa1dSresundermann   }
3409566063dSJacob Faibussowitsch   PetscCall(VecScatterDestroy(&Discat));
3419566063dSJacob Faibussowitsch   PetscCall(VecDestroy(&Diseq));
342aad13602SShrirang Abhyankar   PetscFunctionReturn(0);
343aad13602SShrirang Abhyankar }
344aad13602SShrirang Abhyankar 
345628da978Sresundermann /* Evaluate
346628da978Sresundermann    h = [ x0^2 - x1;
347628da978Sresundermann          1 -(x0^2 - x1)]
348628da978Sresundermann */
349d71ae5a4SJacob Faibussowitsch PetscErrorCode FormInequalityConstraints(Tao tao, Vec X, Vec CI, void *ctx)
350d71ae5a4SJacob Faibussowitsch {
351aad13602SShrirang Abhyankar   const PetscScalar *x;
352aad13602SShrirang Abhyankar   PetscScalar        ci;
353aad13602SShrirang Abhyankar   MPI_Comm           comm;
354aad13602SShrirang Abhyankar   PetscMPIInt        rank;
355aad13602SShrirang Abhyankar   AppCtx            *user = (AppCtx *)ctx;
356aad13602SShrirang Abhyankar   Vec                Xseq = user->Xseq;
357aad13602SShrirang Abhyankar   VecScatter         scat = user->scat;
358aad13602SShrirang Abhyankar 
359aad13602SShrirang Abhyankar   PetscFunctionBegin;
3609566063dSJacob Faibussowitsch   PetscCall(PetscObjectGetComm((PetscObject)tao, &comm));
3619566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Comm_rank(comm, &rank));
362aad13602SShrirang Abhyankar 
3639566063dSJacob Faibussowitsch   PetscCall(VecScatterBegin(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD));
3649566063dSJacob Faibussowitsch   PetscCall(VecScatterEnd(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD));
365aad13602SShrirang Abhyankar 
366dd400576SPatrick Sanan   if (rank == 0) {
3679566063dSJacob Faibussowitsch     PetscCall(VecGetArrayRead(Xseq, &x));
368aad13602SShrirang Abhyankar     ci = x[0] * x[0] - x[1];
3699566063dSJacob Faibussowitsch     PetscCall(VecSetValue(CI, 0, ci, INSERT_VALUES));
370aad13602SShrirang Abhyankar     ci = -x[0] * x[0] + x[1] + 1.0;
3719566063dSJacob Faibussowitsch     PetscCall(VecSetValue(CI, 1, ci, INSERT_VALUES));
3729566063dSJacob Faibussowitsch     PetscCall(VecRestoreArrayRead(Xseq, &x));
373aad13602SShrirang Abhyankar   }
3749566063dSJacob Faibussowitsch   PetscCall(VecAssemblyBegin(CI));
3759566063dSJacob Faibussowitsch   PetscCall(VecAssemblyEnd(CI));
376aad13602SShrirang Abhyankar   PetscFunctionReturn(0);
377aad13602SShrirang Abhyankar }
378aad13602SShrirang Abhyankar 
379628da978Sresundermann /* Evaluate
380628da978Sresundermann    g = [ x0^2 + x1 - 2]
381628da978Sresundermann */
382d71ae5a4SJacob Faibussowitsch PetscErrorCode FormEqualityConstraints(Tao tao, Vec X, Vec CE, void *ctx)
383d71ae5a4SJacob Faibussowitsch {
384aad13602SShrirang Abhyankar   const PetscScalar *x;
385aad13602SShrirang Abhyankar   PetscScalar        ce;
386aad13602SShrirang Abhyankar   MPI_Comm           comm;
387aad13602SShrirang Abhyankar   PetscMPIInt        rank;
388aad13602SShrirang Abhyankar   AppCtx            *user = (AppCtx *)ctx;
389aad13602SShrirang Abhyankar   Vec                Xseq = user->Xseq;
390aad13602SShrirang Abhyankar   VecScatter         scat = user->scat;
391aad13602SShrirang Abhyankar 
392aad13602SShrirang Abhyankar   PetscFunctionBegin;
3939566063dSJacob Faibussowitsch   PetscCall(PetscObjectGetComm((PetscObject)tao, &comm));
3949566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Comm_rank(comm, &rank));
395aad13602SShrirang Abhyankar 
3969566063dSJacob Faibussowitsch   PetscCall(VecScatterBegin(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD));
3979566063dSJacob Faibussowitsch   PetscCall(VecScatterEnd(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD));
398aad13602SShrirang Abhyankar 
399dd400576SPatrick Sanan   if (rank == 0) {
4009566063dSJacob Faibussowitsch     PetscCall(VecGetArrayRead(Xseq, &x));
401aad13602SShrirang Abhyankar     ce = x[0] * x[0] + x[1] - 2.0;
4029566063dSJacob Faibussowitsch     PetscCall(VecSetValue(CE, 0, ce, INSERT_VALUES));
4039566063dSJacob Faibussowitsch     PetscCall(VecRestoreArrayRead(Xseq, &x));
404aad13602SShrirang Abhyankar   }
4059566063dSJacob Faibussowitsch   PetscCall(VecAssemblyBegin(CE));
4069566063dSJacob Faibussowitsch   PetscCall(VecAssemblyEnd(CE));
407aad13602SShrirang Abhyankar   PetscFunctionReturn(0);
408aad13602SShrirang Abhyankar }
409aad13602SShrirang Abhyankar 
410628da978Sresundermann /*
411628da978Sresundermann   grad h = [  2*x0, -1;
412628da978Sresundermann              -2*x0,  1]
413628da978Sresundermann */
414d71ae5a4SJacob Faibussowitsch PetscErrorCode FormInequalityJacobian(Tao tao, Vec X, Mat JI, Mat JIpre, void *ctx)
415d71ae5a4SJacob Faibussowitsch {
416aad13602SShrirang Abhyankar   AppCtx            *user = (AppCtx *)ctx;
417f560b561SHong Zhang   PetscInt           zero = 0, one = 1, cols[2];
418aad13602SShrirang Abhyankar   PetscScalar        vals[2];
419aad13602SShrirang Abhyankar   const PetscScalar *x;
420aad13602SShrirang Abhyankar   Vec                Xseq = user->Xseq;
421aad13602SShrirang Abhyankar   VecScatter         scat = user->scat;
422aad13602SShrirang Abhyankar   MPI_Comm           comm;
423aad13602SShrirang Abhyankar   PetscMPIInt        rank;
424aad13602SShrirang Abhyankar 
425aad13602SShrirang Abhyankar   PetscFunctionBegin;
4269566063dSJacob Faibussowitsch   PetscCall(PetscObjectGetComm((PetscObject)tao, &comm));
4279566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Comm_rank(comm, &rank));
4289566063dSJacob Faibussowitsch   PetscCall(VecScatterBegin(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD));
4299566063dSJacob Faibussowitsch   PetscCall(VecScatterEnd(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD));
430aad13602SShrirang Abhyankar 
4319566063dSJacob Faibussowitsch   PetscCall(VecGetArrayRead(Xseq, &x));
432c5853193SPierre Jolivet   if (rank == 0) {
4339371c9d4SSatish Balay     cols[0] = 0;
4349371c9d4SSatish Balay     cols[1] = 1;
4359371c9d4SSatish Balay     vals[0] = 2 * x[0];
4369371c9d4SSatish Balay     vals[1] = -1.0;
4379566063dSJacob Faibussowitsch     PetscCall(MatSetValues(JI, 1, &zero, 2, cols, vals, INSERT_VALUES));
4389371c9d4SSatish Balay     vals[0] = -2 * x[0];
4399371c9d4SSatish Balay     vals[1] = 1.0;
4409566063dSJacob Faibussowitsch     PetscCall(MatSetValues(JI, 1, &one, 2, cols, vals, INSERT_VALUES));
441aad13602SShrirang Abhyankar   }
4429566063dSJacob Faibussowitsch   PetscCall(VecRestoreArrayRead(Xseq, &x));
4439566063dSJacob Faibussowitsch   PetscCall(MatAssemblyBegin(JI, MAT_FINAL_ASSEMBLY));
4449566063dSJacob Faibussowitsch   PetscCall(MatAssemblyEnd(JI, MAT_FINAL_ASSEMBLY));
445aad13602SShrirang Abhyankar   PetscFunctionReturn(0);
446aad13602SShrirang Abhyankar }
447aad13602SShrirang Abhyankar 
448628da978Sresundermann /*
449628da978Sresundermann   grad g = [2*x0
450628da978Sresundermann              1.0 ]
451628da978Sresundermann */
452d71ae5a4SJacob Faibussowitsch PetscErrorCode FormEqualityJacobian(Tao tao, Vec X, Mat JE, Mat JEpre, void *ctx)
453d71ae5a4SJacob Faibussowitsch {
454f560b561SHong Zhang   PetscInt           zero = 0, cols[2];
455aad13602SShrirang Abhyankar   PetscScalar        vals[2];
456aad13602SShrirang Abhyankar   const PetscScalar *x;
457aad13602SShrirang Abhyankar   PetscMPIInt        rank;
458aad13602SShrirang Abhyankar   MPI_Comm           comm;
459aad13602SShrirang Abhyankar 
460aad13602SShrirang Abhyankar   PetscFunctionBegin;
4619566063dSJacob Faibussowitsch   PetscCall(PetscObjectGetComm((PetscObject)tao, &comm));
4629566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Comm_rank(comm, &rank));
463aad13602SShrirang Abhyankar 
464dd400576SPatrick Sanan   if (rank == 0) {
4659566063dSJacob Faibussowitsch     PetscCall(VecGetArrayRead(X, &x));
4669371c9d4SSatish Balay     cols[0] = 0;
4679371c9d4SSatish Balay     cols[1] = 1;
4689371c9d4SSatish Balay     vals[0] = 2 * x[0];
4699371c9d4SSatish Balay     vals[1] = 1.0;
4709566063dSJacob Faibussowitsch     PetscCall(MatSetValues(JE, 1, &zero, 2, cols, vals, INSERT_VALUES));
4719566063dSJacob Faibussowitsch     PetscCall(VecRestoreArrayRead(X, &x));
472aad13602SShrirang Abhyankar   }
4739566063dSJacob Faibussowitsch   PetscCall(MatAssemblyBegin(JE, MAT_FINAL_ASSEMBLY));
4749566063dSJacob Faibussowitsch   PetscCall(MatAssemblyEnd(JE, MAT_FINAL_ASSEMBLY));
475aad13602SShrirang Abhyankar   PetscFunctionReturn(0);
476aad13602SShrirang Abhyankar }
477aad13602SShrirang Abhyankar 
478aad13602SShrirang Abhyankar /*TEST
479aad13602SShrirang Abhyankar 
480aad13602SShrirang Abhyankar    build:
481f560b561SHong Zhang       requires: !complex !defined(PETSC_USE_CXX)
482aad13602SShrirang Abhyankar 
483aad13602SShrirang Abhyankar    test:
4849a09327dSAlp Dener       args: -tao_converged_reason -tao_gatol 1.e-6 -tao_type pdipm -tao_pdipm_kkt_shift_pd
485f560b561SHong Zhang       requires: mumps
4869a09327dSAlp Dener       filter: sed  -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"
487aad13602SShrirang Abhyankar 
488aad13602SShrirang Abhyankar    test:
489aad13602SShrirang Abhyankar       suffix: 2
4909a09327dSAlp Dener       args: -tao_converged_reason
4919a09327dSAlp Dener       filter: sed  -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"
492aad13602SShrirang Abhyankar 
4938e58fa1dSresundermann    test:
4948e58fa1dSresundermann       suffix: 3
4958e58fa1dSresundermann       args: -tao_converged_reason -no_eq
4969a09327dSAlp Dener       filter: sed  -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"
4978e58fa1dSresundermann 
4988e58fa1dSresundermann    test:
4998e58fa1dSresundermann       suffix: 4
5009a09327dSAlp Dener       args: -tao_converged_reason -tao_almm_type classic
5019a09327dSAlp Dener       requires: !single
5029a09327dSAlp Dener       filter: sed  -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"
5038e58fa1dSresundermann 
504661095bbSAlp Dener    test:
505661095bbSAlp Dener       suffix: 5
5069a09327dSAlp Dener       args: -tao_converged_reason -tao_almm_type classic -no_eq
507*b9e5a4c7SSatish Balay       requires: !single !defined(PETSCTEST_VALGRIND)
5089a09327dSAlp Dener       filter: sed  -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"
509661095bbSAlp Dener 
510661095bbSAlp Dener    test:
511661095bbSAlp Dener       suffix: 6
5129a09327dSAlp Dener       args: -tao_converged_reason -tao_almm_subsolver_tao_type bqnktr
5139a09327dSAlp Dener       filter: sed  -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"
514661095bbSAlp Dener 
515661095bbSAlp Dener    test:
516661095bbSAlp Dener       suffix: 7
5179a09327dSAlp Dener       args: -tao_converged_reason -tao_almm_subsolver_tao_type bncg
5189a09327dSAlp Dener       filter: sed  -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"
519661095bbSAlp Dener 
520661095bbSAlp Dener    test:
521661095bbSAlp Dener       suffix: 8
522661095bbSAlp Dener       nsize: 2
5239a09327dSAlp Dener       args: -tao_converged_reason
5249a09327dSAlp Dener       filter: sed  -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"
525661095bbSAlp Dener 
526661095bbSAlp Dener    test:
527661095bbSAlp Dener       suffix: 9
528661095bbSAlp Dener       nsize: 2
5299a09327dSAlp Dener       args: -tao_converged_reason -vec_type cuda -mat_type aijcusparse
5309a09327dSAlp Dener       requires: cuda
5319a09327dSAlp Dener       filter: sed  -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"
532661095bbSAlp Dener 
533aad13602SShrirang Abhyankar TEST*/
534