xref: /petsc/src/tao/constrained/tutorials/ex1.c (revision 910b42cbfb1b6524ef4594118dd5e013aee7da5d)
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\
2710978b7dSBarry Smith   -tao_monitor_constraint_norm : convergence monitor with constraint norm \n\
2810978b7dSBarry Smith   -tao_monitor_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 *);
4983f04556SHong Zhang PetscErrorCode FormPDIPMHessian(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 
553ba16761SJacob Faibussowitsch int 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));
10583f04556SHong Zhang     PetscCall(PCSetType(pc, PCCHOLESKY));
1069a09327dSAlp Dener     PetscCall(KSPSetFromOptions(ksp));
10783f04556SHong Zhang     PetscCall(TaoSetHessian(tao, user.H, user.H, FormPDIPMHessian, (void *)&user));
1089a09327dSAlp Dener   }
1099a09327dSAlp Dener 
1109a09327dSAlp Dener   /* Print out an initial view of the problem */
111f560b561SHong Zhang   if (user.initview) {
1129566063dSJacob Faibussowitsch     PetscCall(TaoSetUp(tao));
1139566063dSJacob Faibussowitsch     PetscCall(VecDuplicate(user.x, &G));
1149566063dSJacob Faibussowitsch     PetscCall(FormFunctionGradient(tao, user.x, &f, G, (void *)&user));
11583f04556SHong Zhang     PetscCall(FormPDIPMHessian(tao, user.x, user.H, user.H, (void *)&user));
1169566063dSJacob Faibussowitsch     PetscCall(PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_WORLD));
11763a3b9bcSJacob Faibussowitsch     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nInitial point X:\n"));
1189566063dSJacob Faibussowitsch     PetscCall(VecView(user.x, PETSC_VIEWER_STDOUT_WORLD));
11963a3b9bcSJacob Faibussowitsch     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nInitial objective f(x) = %g\n", (double)f));
12063a3b9bcSJacob Faibussowitsch     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nInitial gradient and Hessian:\n"));
1219566063dSJacob Faibussowitsch     PetscCall(VecView(G, PETSC_VIEWER_STDOUT_WORLD));
1229566063dSJacob Faibussowitsch     PetscCall(MatView(user.H, PETSC_VIEWER_STDOUT_WORLD));
1239566063dSJacob Faibussowitsch     PetscCall(VecDestroy(&G));
1249566063dSJacob Faibussowitsch     PetscCall(FormInequalityJacobian(tao, user.x, user.Ai, user.Ai, (void *)&user));
1259566063dSJacob Faibussowitsch     PetscCall(MatCreateVecs(user.Ai, NULL, &CI));
1269566063dSJacob Faibussowitsch     PetscCall(FormInequalityConstraints(tao, user.x, CI, (void *)&user));
12763a3b9bcSJacob Faibussowitsch     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nInitial inequality constraints and Jacobian:\n"));
1289566063dSJacob Faibussowitsch     PetscCall(VecView(CI, PETSC_VIEWER_STDOUT_WORLD));
1299566063dSJacob Faibussowitsch     PetscCall(MatView(user.Ai, PETSC_VIEWER_STDOUT_WORLD));
1309566063dSJacob Faibussowitsch     PetscCall(VecDestroy(&CI));
131f560b561SHong Zhang     if (!user.noeqflag) {
1329566063dSJacob Faibussowitsch       PetscCall(FormEqualityJacobian(tao, user.x, user.Ae, user.Ae, (void *)&user));
1339566063dSJacob Faibussowitsch       PetscCall(MatCreateVecs(user.Ae, NULL, &CE));
1349566063dSJacob Faibussowitsch       PetscCall(FormEqualityConstraints(tao, user.x, CE, (void *)&user));
13563a3b9bcSJacob Faibussowitsch       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nInitial equality constraints and Jacobian:\n"));
1369566063dSJacob Faibussowitsch       PetscCall(VecView(CE, PETSC_VIEWER_STDOUT_WORLD));
1379566063dSJacob Faibussowitsch       PetscCall(MatView(user.Ae, PETSC_VIEWER_STDOUT_WORLD));
1389566063dSJacob Faibussowitsch       PetscCall(VecDestroy(&CE));
139f560b561SHong Zhang     }
1409566063dSJacob Faibussowitsch     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n"));
1419566063dSJacob Faibussowitsch     PetscCall(PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_WORLD));
142f560b561SHong Zhang   }
143aad13602SShrirang Abhyankar 
1449566063dSJacob Faibussowitsch   PetscCall(TaoSolve(tao));
1459566063dSJacob Faibussowitsch   PetscCall(TaoGetSolution(tao, &x));
1469a09327dSAlp Dener   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Found solution:\n"));
1479566063dSJacob Faibussowitsch   PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD));
148aad13602SShrirang Abhyankar 
149aad13602SShrirang Abhyankar   /* Free objects */
1509566063dSJacob Faibussowitsch   PetscCall(DestroyProblem(&user));
1519566063dSJacob Faibussowitsch   PetscCall(TaoDestroy(&tao));
1529566063dSJacob Faibussowitsch   PetscCall(PetscFinalize());
1533ba16761SJacob Faibussowitsch   return PETSC_SUCCESS;
154aad13602SShrirang Abhyankar }
155aad13602SShrirang Abhyankar 
156d71ae5a4SJacob Faibussowitsch PetscErrorCode InitializeProblem(AppCtx *user)
157d71ae5a4SJacob Faibussowitsch {
158aad13602SShrirang Abhyankar   PetscMPIInt size;
159aad13602SShrirang Abhyankar   PetscMPIInt rank;
160aad13602SShrirang Abhyankar   PetscInt    nloc, neloc, niloc;
161aad13602SShrirang Abhyankar 
162aad13602SShrirang Abhyankar   PetscFunctionBegin;
1639566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
1649566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
1658e58fa1dSresundermann   user->noeqflag = PETSC_FALSE;
1669566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL, NULL, "-no_eq", &user->noeqflag, NULL));
167f560b561SHong Zhang   user->initview = PETSC_FALSE;
1689566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL, NULL, "-init_view", &user->initview, NULL));
169628da978Sresundermann 
170f560b561SHong Zhang   /* Tell user the correct solution, not an error checking */
1719a09327dSAlp Dener   if (!user->noeqflag) {
1729a09327dSAlp Dener     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Expected solution: f(1, 1) = -2\n"));
1739a09327dSAlp Dener   } else {
1749a09327dSAlp Dener     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Expected solution (-no_eq): f(1.73205, 2) = -7.3923\n"));
1758e58fa1dSresundermann   }
176aad13602SShrirang Abhyankar 
1772d4ee042Sprj-   /* create vector x and set initial values */
178aad13602SShrirang Abhyankar   user->n = 2; /* global length */
179f560b561SHong Zhang   nloc    = (size == 1) ? user->n : 1;
1809566063dSJacob Faibussowitsch   PetscCall(VecCreate(PETSC_COMM_WORLD, &user->x));
1819566063dSJacob Faibussowitsch   PetscCall(VecSetSizes(user->x, nloc, user->n));
1829566063dSJacob Faibussowitsch   PetscCall(VecSetFromOptions(user->x));
1839a09327dSAlp Dener   PetscCall(VecSet(user->x, 0.0));
184aad13602SShrirang Abhyankar 
185aad13602SShrirang Abhyankar   /* create and set lower and upper bound vectors */
1869566063dSJacob Faibussowitsch   PetscCall(VecDuplicate(user->x, &user->xl));
1879566063dSJacob Faibussowitsch   PetscCall(VecDuplicate(user->x, &user->xu));
1889566063dSJacob Faibussowitsch   PetscCall(VecSet(user->xl, -1.0));
1899566063dSJacob Faibussowitsch   PetscCall(VecSet(user->xu, 2.0));
190aad13602SShrirang Abhyankar 
191aad13602SShrirang Abhyankar   /* create scater to zero */
1929566063dSJacob Faibussowitsch   PetscCall(VecScatterCreateToZero(user->x, &user->scat, &user->Xseq));
193aad13602SShrirang Abhyankar 
194aad13602SShrirang Abhyankar   user->ne = 1;
195628da978Sresundermann   user->ni = 2;
196aad13602SShrirang Abhyankar   neloc    = (rank == 0) ? user->ne : 0;
197f560b561SHong Zhang   niloc    = (size == 1) ? user->ni : 1;
198aad13602SShrirang Abhyankar 
1998e58fa1dSresundermann   if (!user->noeqflag) {
2009566063dSJacob Faibussowitsch     PetscCall(VecCreate(PETSC_COMM_WORLD, &user->ce)); /* a 1x1 vec for equality constraints */
2019566063dSJacob Faibussowitsch     PetscCall(VecSetSizes(user->ce, neloc, user->ne));
2029566063dSJacob Faibussowitsch     PetscCall(VecSetFromOptions(user->ce));
2039566063dSJacob Faibussowitsch     PetscCall(VecSetUp(user->ce));
2048e58fa1dSresundermann   }
205628da978Sresundermann 
2069566063dSJacob Faibussowitsch   PetscCall(VecCreate(PETSC_COMM_WORLD, &user->ci)); /* a 2x1 vec for inequality constraints */
2079566063dSJacob Faibussowitsch   PetscCall(VecSetSizes(user->ci, niloc, user->ni));
2089566063dSJacob Faibussowitsch   PetscCall(VecSetFromOptions(user->ci));
2099566063dSJacob Faibussowitsch   PetscCall(VecSetUp(user->ci));
210aad13602SShrirang Abhyankar 
211baca6076SPierre Jolivet   /* nexn & nixn matrices for equally and inequalty constraints */
2128e58fa1dSresundermann   if (!user->noeqflag) {
2139566063dSJacob Faibussowitsch     PetscCall(MatCreate(PETSC_COMM_WORLD, &user->Ae));
2149566063dSJacob Faibussowitsch     PetscCall(MatSetSizes(user->Ae, neloc, nloc, user->ne, user->n));
2159566063dSJacob Faibussowitsch     PetscCall(MatSetFromOptions(user->Ae));
2169566063dSJacob Faibussowitsch     PetscCall(MatSetUp(user->Ae));
2178e58fa1dSresundermann   }
2188e58fa1dSresundermann 
2199566063dSJacob Faibussowitsch   PetscCall(MatCreate(PETSC_COMM_WORLD, &user->Ai));
2209566063dSJacob Faibussowitsch   PetscCall(MatSetSizes(user->Ai, niloc, nloc, user->ni, user->n));
2219566063dSJacob Faibussowitsch   PetscCall(MatSetFromOptions(user->Ai));
2229566063dSJacob Faibussowitsch   PetscCall(MatSetUp(user->Ai));
223628da978Sresundermann 
2249566063dSJacob Faibussowitsch   PetscCall(MatCreate(PETSC_COMM_WORLD, &user->H));
2259566063dSJacob Faibussowitsch   PetscCall(MatSetSizes(user->H, nloc, nloc, user->n, user->n));
2269566063dSJacob Faibussowitsch   PetscCall(MatSetFromOptions(user->H));
2279566063dSJacob Faibussowitsch   PetscCall(MatSetUp(user->H));
2283ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
229aad13602SShrirang Abhyankar }
230aad13602SShrirang Abhyankar 
231d71ae5a4SJacob Faibussowitsch PetscErrorCode DestroyProblem(AppCtx *user)
232d71ae5a4SJacob Faibussowitsch {
233aad13602SShrirang Abhyankar   PetscFunctionBegin;
23448a46eb9SPierre Jolivet   if (!user->noeqflag) PetscCall(MatDestroy(&user->Ae));
2359566063dSJacob Faibussowitsch   PetscCall(MatDestroy(&user->Ai));
2369566063dSJacob Faibussowitsch   PetscCall(MatDestroy(&user->H));
237aad13602SShrirang Abhyankar 
2389566063dSJacob Faibussowitsch   PetscCall(VecDestroy(&user->x));
23948a46eb9SPierre Jolivet   if (!user->noeqflag) PetscCall(VecDestroy(&user->ce));
2409566063dSJacob Faibussowitsch   PetscCall(VecDestroy(&user->ci));
2419566063dSJacob Faibussowitsch   PetscCall(VecDestroy(&user->xl));
2429566063dSJacob Faibussowitsch   PetscCall(VecDestroy(&user->xu));
2439566063dSJacob Faibussowitsch   PetscCall(VecDestroy(&user->Xseq));
2449566063dSJacob Faibussowitsch   PetscCall(VecScatterDestroy(&user->scat));
2453ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
246aad13602SShrirang Abhyankar }
247aad13602SShrirang Abhyankar 
248628da978Sresundermann /* Evaluate
249628da978Sresundermann    f(x) = (x0 - 2)^2 + (x1 - 2)^2 - 2*(x0 + x1)
25083f04556SHong Zhang    G = grad f = [2*(x0 - 2) - 2, 2*(x1 - 2) - 2] = 2*X - 6
251aad13602SShrirang Abhyankar */
252d71ae5a4SJacob Faibussowitsch PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ctx)
253d71ae5a4SJacob Faibussowitsch {
254aad13602SShrirang Abhyankar   const PetscScalar *x;
255aad13602SShrirang Abhyankar   MPI_Comm           comm;
256aad13602SShrirang Abhyankar   PetscMPIInt        rank;
257aad13602SShrirang Abhyankar   PetscReal          fin;
258aad13602SShrirang Abhyankar   AppCtx            *user = (AppCtx *)ctx;
259aad13602SShrirang Abhyankar   Vec                Xseq = user->Xseq;
260aad13602SShrirang Abhyankar   VecScatter         scat = user->scat;
261aad13602SShrirang Abhyankar 
262aad13602SShrirang Abhyankar   PetscFunctionBegin;
26383f04556SHong Zhang   /* f = (x0 - 2)^2 + (x1 - 2)^2 - 2*(x0 + x1) */
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]);
2749566063dSJacob Faibussowitsch     PetscCall(VecRestoreArrayRead(Xseq, &x));
275aad13602SShrirang Abhyankar   }
276712fec58SPierre Jolivet   PetscCall(MPIU_Allreduce(&fin, f, 1, MPIU_REAL, MPIU_SUM, comm));
27783f04556SHong Zhang 
27883f04556SHong Zhang   /* G = 2*X - 6 */
27983f04556SHong Zhang   PetscCall(VecSet(G, -6.0));
28083f04556SHong Zhang   PetscCall(VecAXPY(G, 2.0, X));
2813ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
282aad13602SShrirang Abhyankar }
283aad13602SShrirang Abhyankar 
28483f04556SHong Zhang /* Evaluate PDIPM Hessian, see Eqn(22) in http://doi.org/10.1049/gtd2.12708
28583f04556SHong Zhang    H = Wxx = fxx        + Jacobian(grad g^T*DE)   - Jacobian(grad h^T*DI)]
28683f04556SHong Zhang            = fxx        + Jacobin([2*x0; 1]de[0]) - Jacobian([2*x0, -2*x0; -1, 1][di[0] di[1]]^T)
28783f04556SHong Zhang            = [2 0; 0 2] + [2*de[0]  0;      0  0] - [2*di[0]-2*di[1], 0; 0, 0]
288628da978Sresundermann            = [ 2*(1+de[0]-di[0]+di[1]), 0;
289628da978Sresundermann                           0,            2]
290628da978Sresundermann */
29183f04556SHong Zhang PetscErrorCode FormPDIPMHessian(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));
3423ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
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));
3763ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
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));
4073ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
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));
4453ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
446aad13602SShrirang Abhyankar }
447aad13602SShrirang Abhyankar 
448628da978Sresundermann /*
44983f04556SHong Zhang   grad g = [2*x0, 1.0]
450628da978Sresundermann */
451d71ae5a4SJacob Faibussowitsch PetscErrorCode FormEqualityJacobian(Tao tao, Vec X, Mat JE, Mat JEpre, void *ctx)
452d71ae5a4SJacob Faibussowitsch {
453f560b561SHong Zhang   PetscInt           zero = 0, cols[2];
454aad13602SShrirang Abhyankar   PetscScalar        vals[2];
455aad13602SShrirang Abhyankar   const PetscScalar *x;
456aad13602SShrirang Abhyankar   PetscMPIInt        rank;
457aad13602SShrirang Abhyankar   MPI_Comm           comm;
458aad13602SShrirang Abhyankar 
459aad13602SShrirang Abhyankar   PetscFunctionBegin;
4609566063dSJacob Faibussowitsch   PetscCall(PetscObjectGetComm((PetscObject)tao, &comm));
4619566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Comm_rank(comm, &rank));
462aad13602SShrirang Abhyankar 
463dd400576SPatrick Sanan   if (rank == 0) {
4649566063dSJacob Faibussowitsch     PetscCall(VecGetArrayRead(X, &x));
4659371c9d4SSatish Balay     cols[0] = 0;
4669371c9d4SSatish Balay     cols[1] = 1;
4679371c9d4SSatish Balay     vals[0] = 2 * x[0];
4689371c9d4SSatish Balay     vals[1] = 1.0;
4699566063dSJacob Faibussowitsch     PetscCall(MatSetValues(JE, 1, &zero, 2, cols, vals, INSERT_VALUES));
4709566063dSJacob Faibussowitsch     PetscCall(VecRestoreArrayRead(X, &x));
471aad13602SShrirang Abhyankar   }
4729566063dSJacob Faibussowitsch   PetscCall(MatAssemblyBegin(JE, MAT_FINAL_ASSEMBLY));
4739566063dSJacob Faibussowitsch   PetscCall(MatAssemblyEnd(JE, MAT_FINAL_ASSEMBLY));
4743ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
475aad13602SShrirang Abhyankar }
476aad13602SShrirang Abhyankar 
477aad13602SShrirang Abhyankar /*TEST
478aad13602SShrirang Abhyankar 
479aad13602SShrirang Abhyankar    build:
480f560b561SHong Zhang       requires: !complex !defined(PETSC_USE_CXX)
481aad13602SShrirang Abhyankar 
482aad13602SShrirang Abhyankar    test:
4839a09327dSAlp Dener       args: -tao_converged_reason -tao_gatol 1.e-6 -tao_type pdipm -tao_pdipm_kkt_shift_pd
484*910b42cbSStefano Zampini       requires: mumps !single
4859a09327dSAlp Dener       filter: sed  -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"
486aad13602SShrirang Abhyankar 
487aad13602SShrirang Abhyankar    test:
48883f04556SHong Zhang       suffix: pdipm_2
48983f04556SHong Zhang       requires: mumps
49083f04556SHong Zhang       nsize: 2
49183f04556SHong Zhang       args: -tao_converged_reason -tao_gatol 1.e-6 -tao_type pdipm
49283f04556SHong Zhang       filter: sed  -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"
49383f04556SHong Zhang 
49483f04556SHong Zhang    test:
495aad13602SShrirang Abhyankar       suffix: 2
4969a09327dSAlp Dener       args: -tao_converged_reason
4979a09327dSAlp Dener       filter: sed  -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"
498aad13602SShrirang Abhyankar 
4998e58fa1dSresundermann    test:
5008e58fa1dSresundermann       suffix: 3
5018e58fa1dSresundermann       args: -tao_converged_reason -no_eq
5029a09327dSAlp Dener       filter: sed  -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"
5038e58fa1dSresundermann 
5048e58fa1dSresundermann    test:
5058e58fa1dSresundermann       suffix: 4
5069a09327dSAlp Dener       args: -tao_converged_reason -tao_almm_type classic
5079a09327dSAlp Dener       requires: !single
5089a09327dSAlp Dener       filter: sed  -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"
5098e58fa1dSresundermann 
510661095bbSAlp Dener    test:
511661095bbSAlp Dener       suffix: 5
512910a5171SToby Isaac       args: -tao_converged_reason -tao_almm_type classic -no_eq -tao_almm_subsolver_tao_max_it 100
513b9e5a4c7SSatish Balay       requires: !single !defined(PETSCTEST_VALGRIND)
5149a09327dSAlp Dener       filter: sed  -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"
515661095bbSAlp Dener 
516661095bbSAlp Dener    test:
517661095bbSAlp Dener       suffix: 6
5189a09327dSAlp Dener       args: -tao_converged_reason -tao_almm_subsolver_tao_type bqnktr
5199a09327dSAlp Dener       filter: sed  -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"
520661095bbSAlp Dener 
521661095bbSAlp Dener    test:
522661095bbSAlp Dener       suffix: 7
5239a09327dSAlp Dener       args: -tao_converged_reason -tao_almm_subsolver_tao_type bncg
5249a09327dSAlp Dener       filter: sed  -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"
525661095bbSAlp Dener 
526661095bbSAlp Dener    test:
527661095bbSAlp Dener       suffix: 8
528661095bbSAlp Dener       nsize: 2
5299a09327dSAlp Dener       args: -tao_converged_reason
5309a09327dSAlp Dener       filter: sed  -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"
531661095bbSAlp Dener 
532661095bbSAlp Dener    test:
533661095bbSAlp Dener       suffix: 9
534661095bbSAlp Dener       nsize: 2
5359a09327dSAlp Dener       args: -tao_converged_reason -vec_type cuda -mat_type aijcusparse
5369a09327dSAlp Dener       requires: cuda
5379a09327dSAlp Dener       filter: sed  -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"
538661095bbSAlp Dener 
539aad13602SShrirang Abhyankar TEST*/
540