Actual source code: ex1.c

  1: /* Program usage: mpiexec -n 2 ./ex1 [-help] [all TAO options] */

  3: /* ----------------------------------------------------------------------
  4: min f(x) = (x0 - 2)^2 + (x1 - 2)^2 - 2*(x0 + x1)
  5: s.t.  x0^2 + x1 - 2 = 0
  6:       0  <= x0^2 - x1 <= 1
  7:       -1 <= x0, x1 <= 2
  8: -->
  9:       g(x)  = 0
 10:       h(x) >= 0
 11:       -1 <= x0, x1 <= 2
 12: where
 13:       g(x) = x0^2 + x1 - 2
 14:       h(x) = [x0^2 - x1
 15:               1 -(x0^2 - x1)]
 16: ---------------------------------------------------------------------- */

 18: #include <petsctao.h>

 20: static char help[] = "Solves constrained optimization problem using pdipm.\n\
 21: Input parameters include:\n\
 22:   -tao_type pdipm    : sets Tao solver\n\
 23:   -no_eq             : removes the equaility constraints from the problem\n\
 24:   -no_bound          : removes the bound constraints from the problem\n\
 25:   -init_view         : view initial object setup\n\
 26:   -snes_fd           : snes with finite difference Jacobian (needed for pdipm)\n\
 27:   -snes_compare_explicit : compare user Jacobian with finite difference Jacobian \n\
 28:   -tao_monitor_constraint_norm : convergence monitor with constraint norm \n\
 29:   -tao_monitor_solution : view exact solution at each iteration\n\
 30:   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";

 32: /*
 33:    User-defined application context - contains data needed by the application
 34: */
 35: typedef struct {
 36:   PetscInt   n;  /* Global length of x */
 37:   PetscInt   ne; /* Global number of equality constraints */
 38:   PetscInt   ni; /* Global number of inequality constraints */
 39:   PetscBool  noeqflag, noboundflag, initview;
 40:   Vec        x, xl, xu;
 41:   Vec        ce, ci, bl, bu, Xseq;
 42:   Mat        Ae, Ai, H;
 43:   VecScatter scat;
 44: } AppCtx;

 46: /* -------- User-defined Routines --------- */
 47: PetscErrorCode InitializeProblem(AppCtx *);
 48: PetscErrorCode DestroyProblem(AppCtx *);
 49: PetscErrorCode FormFunctionGradient(Tao, Vec, PetscReal *, Vec, void *);
 50: PetscErrorCode FormPDIPMHessian(Tao, Vec, Mat, Mat, void *);
 51: PetscErrorCode FormInequalityConstraints(Tao, Vec, Vec, void *);
 52: PetscErrorCode FormEqualityConstraints(Tao, Vec, Vec, void *);
 53: PetscErrorCode FormInequalityJacobian(Tao, Vec, Mat, Mat, void *);
 54: PetscErrorCode FormEqualityJacobian(Tao, Vec, Mat, Mat, void *);

 56: int main(int argc, char **argv)
 57: {
 58:   Tao         tao;
 59:   KSP         ksp;
 60:   PC          pc;
 61:   AppCtx      user; /* application context */
 62:   Vec         x, G, CI, CE;
 63:   PetscMPIInt size;
 64:   TaoType     type;
 65:   PetscReal   f;
 66:   PetscBool   pdipm, mumps;

 68:   PetscFunctionBeginUser;
 69:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
 70:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 71:   PetscCheck(size <= 2, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "More than 2 processors detected. Example written to use max of 2 processors.");

 73:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "---- Constrained Problem -----\n"));
 74:   PetscCall(InitializeProblem(&user)); /* sets up problem, function below */

 76:   PetscCall(TaoCreate(PETSC_COMM_WORLD, &tao));
 77:   PetscCall(TaoSetType(tao, TAOALMM));
 78:   PetscCall(TaoSetSolution(tao, user.x));
 79:   if (!user.noboundflag) PetscCall(TaoSetVariableBounds(tao, user.xl, user.xu));
 80:   PetscCall(TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, (void *)&user));
 81:   PetscCall(TaoSetTolerances(tao, 1.e-4, 0.0, 0.0));
 82:   PetscCall(TaoSetConstraintTolerances(tao, 1.e-4, 0.0));
 83:   PetscCall(TaoSetMaximumFunctionEvaluations(tao, 1e6));
 84:   PetscCall(TaoSetFromOptions(tao));

 86:   if (!user.noeqflag) {
 87:     PetscCall(TaoSetEqualityConstraintsRoutine(tao, user.ce, FormEqualityConstraints, (void *)&user));
 88:     PetscCall(TaoSetJacobianEqualityRoutine(tao, user.Ae, user.Ae, FormEqualityJacobian, (void *)&user));
 89:   }
 90:   PetscCall(TaoSetInequalityConstraintsRoutine(tao, user.ci, FormInequalityConstraints, (void *)&user));
 91:   PetscCall(TaoSetJacobianInequalityRoutine(tao, user.Ai, user.Ai, FormInequalityJacobian, (void *)&user));

 93:   PetscCall(TaoGetType(tao, &type));
 94:   PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOPDIPM, &pdipm));
 95:   if (pdipm) {
 96:     /*
 97:       PDIPM produces an indefinite KKT matrix with some zeros along the diagonal
 98:       Inverting this indefinite matrix requires PETSc to be configured with MUMPS
 99:     */
100:     PetscCall(PetscHasExternalPackage("mumps", &mumps));
101:     PetscCheck(mumps, PetscObjectComm((PetscObject)tao), PETSC_ERR_SUP, "TAOPDIPM requires PETSc to be configured with MUMPS (--download-mumps)");
102:     PetscCall(TaoGetKSP(tao, &ksp));
103:     PetscCall(KSPGetPC(ksp, &pc));
104:     PetscCall(PCFactorSetMatSolverType(pc, MATSOLVERMUMPS));
105:     PetscCall(KSPSetType(ksp, KSPPREONLY));
106:     PetscCall(PCSetType(pc, PCCHOLESKY));
107:     PetscCall(KSPSetFromOptions(ksp));
108:     PetscCall(TaoSetHessian(tao, user.H, user.H, FormPDIPMHessian, (void *)&user));
109:   }

111:   /* Print out an initial view of the problem */
112:   if (user.initview) {
113:     PetscCall(TaoSetUp(tao));
114:     PetscCall(VecDuplicate(user.x, &G));
115:     PetscCall(FormFunctionGradient(tao, user.x, &f, G, (void *)&user));
116:     if (pdipm) PetscCall(FormPDIPMHessian(tao, user.x, user.H, user.H, (void *)&user));
117:     PetscCall(PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_WORLD));
118:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nInitial point X:\n"));
119:     PetscCall(VecView(user.x, PETSC_VIEWER_STDOUT_WORLD));
120:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nInitial objective f(x) = %g\n", (double)f));
121:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nInitial gradient and Hessian:\n"));
122:     PetscCall(VecView(G, PETSC_VIEWER_STDOUT_WORLD));
123:     if (pdipm) PetscCall(MatView(user.H, PETSC_VIEWER_STDOUT_WORLD));
124:     PetscCall(VecDestroy(&G));
125:     PetscCall(FormInequalityJacobian(tao, user.x, user.Ai, user.Ai, (void *)&user));
126:     PetscCall(MatCreateVecs(user.Ai, NULL, &CI));
127:     PetscCall(FormInequalityConstraints(tao, user.x, CI, (void *)&user));
128:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nInitial inequality constraints and Jacobian:\n"));
129:     PetscCall(VecView(CI, PETSC_VIEWER_STDOUT_WORLD));
130:     PetscCall(MatView(user.Ai, PETSC_VIEWER_STDOUT_WORLD));
131:     PetscCall(VecDestroy(&CI));
132:     if (!user.noeqflag) {
133:       PetscCall(FormEqualityJacobian(tao, user.x, user.Ae, user.Ae, (void *)&user));
134:       PetscCall(MatCreateVecs(user.Ae, NULL, &CE));
135:       PetscCall(FormEqualityConstraints(tao, user.x, CE, (void *)&user));
136:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nInitial equality constraints and Jacobian:\n"));
137:       PetscCall(VecView(CE, PETSC_VIEWER_STDOUT_WORLD));
138:       PetscCall(MatView(user.Ae, PETSC_VIEWER_STDOUT_WORLD));
139:       PetscCall(VecDestroy(&CE));
140:     }
141:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n"));
142:     PetscCall(PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_WORLD));
143:   }

145:   PetscCall(TaoSolve(tao));
146:   PetscCall(TaoGetSolution(tao, &x));
147:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Found solution:\n"));
148:   PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD));

150:   /* Free objects */
151:   PetscCall(DestroyProblem(&user));
152:   PetscCall(TaoDestroy(&tao));
153:   PetscCall(PetscFinalize());
154:   return PETSC_SUCCESS;
155: }

157: PetscErrorCode InitializeProblem(AppCtx *user)
158: {
159:   PetscMPIInt size;
160:   PetscMPIInt rank;
161:   PetscInt    nloc, neloc, niloc;

163:   PetscFunctionBegin;
164:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
165:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
166:   user->noeqflag    = PETSC_FALSE;
167:   user->noboundflag = PETSC_FALSE;
168:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-no_eq", &user->noeqflag, NULL));
169:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-no_bound", &user->noboundflag, NULL));
170:   user->initview = PETSC_FALSE;
171:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-init_view", &user->initview, NULL));

173:   /* Tell user the correct solution, not an error checking */
174:   if (!user->noeqflag) {
175:     /* With equality constraint, bound or no bound makes no different to the solution */
176:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Expected solution: f(1, 1) = -2\n"));
177:   } else {
178:     if (user->noboundflag) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Expected solution (-no_eq, -no_bound): f(2.05655, 3.22938) = -9.05728\n"));
179:     else PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Expected solution (-no_eq): f(1.73205, 2) = -7.3923\n"));
180:   }

182:   /* create vector x and set initial values */
183:   user->n = 2; /* global length */
184:   nloc    = (size == 1) ? user->n : 1;
185:   PetscCall(VecCreate(PETSC_COMM_WORLD, &user->x));
186:   PetscCall(VecSetSizes(user->x, nloc, user->n));
187:   PetscCall(VecSetFromOptions(user->x));
188:   PetscCall(VecSet(user->x, 0.0));

190:   /* create and set lower and upper bound vectors */
191:   PetscCall(VecDuplicate(user->x, &user->xl));
192:   PetscCall(VecDuplicate(user->x, &user->xu));
193:   PetscCall(VecSet(user->xl, -1.0));
194:   PetscCall(VecSet(user->xu, 2.0));

196:   /* create scater to zero */
197:   PetscCall(VecScatterCreateToZero(user->x, &user->scat, &user->Xseq));

199:   user->ne = 1;
200:   user->ni = 2;
201:   neloc    = (rank == 0) ? user->ne : 0;
202:   niloc    = (size == 1) ? user->ni : 1;

204:   if (!user->noeqflag) {
205:     PetscCall(VecCreate(PETSC_COMM_WORLD, &user->ce)); /* a 1x1 vec for equality constraints */
206:     PetscCall(VecSetSizes(user->ce, neloc, user->ne));
207:     PetscCall(VecSetFromOptions(user->ce));
208:     PetscCall(VecSetUp(user->ce));
209:   }

211:   PetscCall(VecCreate(PETSC_COMM_WORLD, &user->ci)); /* a 2x1 vec for inequality constraints */
212:   PetscCall(VecSetSizes(user->ci, niloc, user->ni));
213:   PetscCall(VecSetFromOptions(user->ci));
214:   PetscCall(VecSetUp(user->ci));

216:   /* nexn & nixn matrices for equally and inequalty constraints */
217:   if (!user->noeqflag) {
218:     PetscCall(MatCreate(PETSC_COMM_WORLD, &user->Ae));
219:     PetscCall(MatSetSizes(user->Ae, neloc, nloc, user->ne, user->n));
220:     PetscCall(MatSetFromOptions(user->Ae));
221:     PetscCall(MatSetUp(user->Ae));
222:   }

224:   PetscCall(MatCreate(PETSC_COMM_WORLD, &user->Ai));
225:   PetscCall(MatSetSizes(user->Ai, niloc, nloc, user->ni, user->n));
226:   PetscCall(MatSetFromOptions(user->Ai));
227:   PetscCall(MatSetUp(user->Ai));

229:   PetscCall(MatCreate(PETSC_COMM_WORLD, &user->H));
230:   PetscCall(MatSetSizes(user->H, nloc, nloc, user->n, user->n));
231:   PetscCall(MatSetFromOptions(user->H));
232:   PetscCall(MatSetUp(user->H));
233:   PetscFunctionReturn(PETSC_SUCCESS);
234: }

236: PetscErrorCode DestroyProblem(AppCtx *user)
237: {
238:   PetscFunctionBegin;
239:   if (!user->noeqflag) PetscCall(MatDestroy(&user->Ae));
240:   PetscCall(MatDestroy(&user->Ai));
241:   PetscCall(MatDestroy(&user->H));

243:   PetscCall(VecDestroy(&user->x));
244:   if (!user->noeqflag) PetscCall(VecDestroy(&user->ce));
245:   PetscCall(VecDestroy(&user->ci));
246:   PetscCall(VecDestroy(&user->xl));
247:   PetscCall(VecDestroy(&user->xu));
248:   PetscCall(VecDestroy(&user->Xseq));
249:   PetscCall(VecScatterDestroy(&user->scat));
250:   PetscFunctionReturn(PETSC_SUCCESS);
251: }

253: /* Evaluate
254:    f(x) = (x0 - 2)^2 + (x1 - 2)^2 - 2*(x0 + x1)
255:    G = grad f = [2*(x0 - 2) - 2, 2*(x1 - 2) - 2] = 2*X - 6
256: */
257: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ctx)
258: {
259:   const PetscScalar *x;
260:   MPI_Comm           comm;
261:   PetscMPIInt        rank;
262:   PetscReal          fin;
263:   AppCtx            *user = (AppCtx *)ctx;
264:   Vec                Xseq = user->Xseq;
265:   VecScatter         scat = user->scat;

267:   PetscFunctionBegin;
268:   /* f = (x0 - 2)^2 + (x1 - 2)^2 - 2*(x0 + x1) */
269:   PetscCall(PetscObjectGetComm((PetscObject)tao, &comm));
270:   PetscCallMPI(MPI_Comm_rank(comm, &rank));

272:   PetscCall(VecScatterBegin(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD));
273:   PetscCall(VecScatterEnd(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD));

275:   fin = 0.0;
276:   if (rank == 0) {
277:     PetscCall(VecGetArrayRead(Xseq, &x));
278:     fin = (x[0] - 2.0) * (x[0] - 2.0) + (x[1] - 2.0) * (x[1] - 2.0) - 2.0 * (x[0] + x[1]);
279:     PetscCall(VecRestoreArrayRead(Xseq, &x));
280:   }
281:   PetscCallMPI(MPIU_Allreduce(&fin, f, 1, MPIU_REAL, MPIU_SUM, comm));

283:   /* G = 2*X - 6 */
284:   PetscCall(VecSet(G, -6.0));
285:   PetscCall(VecAXPY(G, 2.0, X));
286:   PetscFunctionReturn(PETSC_SUCCESS);
287: }

289: /* Evaluate PDIPM Hessian, see Eqn(22) in http://doi.org/10.1049/gtd2.12708
290:    H = Wxx = fxx        + Jacobian(grad g^T*DE)   - Jacobian(grad h^T*DI)]
291:            = fxx        + Jacobin([2*x0; 1]de[0]) - Jacobian([2*x0, -2*x0; -1, 1][di[0] di[1]]^T)
292:            = [2 0; 0 2] + [2*de[0]  0;      0  0] - [2*di[0]-2*di[1], 0; 0, 0]
293:            = [ 2*(1+de[0]-di[0]+di[1]), 0;
294:                           0,            2]
295: */
296: PetscErrorCode FormPDIPMHessian(Tao tao, Vec x, Mat H, Mat Hpre, void *ctx)
297: {
298:   AppCtx            *user = (AppCtx *)ctx;
299:   Vec                DE, DI;
300:   const PetscScalar *de, *di;
301:   PetscInt           zero = 0, one = 1;
302:   PetscScalar        two = 2.0;
303:   PetscScalar        val = 0.0;
304:   Vec                Deseq, Diseq;
305:   VecScatter         Descat, Discat;
306:   PetscMPIInt        rank;
307:   MPI_Comm           comm;

309:   PetscFunctionBegin;
310:   PetscCall(TaoGetDualVariables(tao, &DE, &DI));

312:   PetscCall(PetscObjectGetComm((PetscObject)tao, &comm));
313:   PetscCallMPI(MPI_Comm_rank(comm, &rank));

315:   if (!user->noeqflag) {
316:     PetscCall(VecScatterCreateToZero(DE, &Descat, &Deseq));
317:     PetscCall(VecScatterBegin(Descat, DE, Deseq, INSERT_VALUES, SCATTER_FORWARD));
318:     PetscCall(VecScatterEnd(Descat, DE, Deseq, INSERT_VALUES, SCATTER_FORWARD));
319:   }
320:   PetscCall(VecScatterCreateToZero(DI, &Discat, &Diseq));
321:   PetscCall(VecScatterBegin(Discat, DI, Diseq, INSERT_VALUES, SCATTER_FORWARD));
322:   PetscCall(VecScatterEnd(Discat, DI, Diseq, INSERT_VALUES, SCATTER_FORWARD));

324:   if (rank == 0) {
325:     if (!user->noeqflag) { PetscCall(VecGetArrayRead(Deseq, &de)); /* places equality constraint dual into array */ }
326:     PetscCall(VecGetArrayRead(Diseq, &di)); /* places inequality constraint dual into array */

328:     if (!user->noeqflag) {
329:       val = 2.0 * (1 + de[0] - di[0] + di[1]);
330:       PetscCall(VecRestoreArrayRead(Deseq, &de));
331:       PetscCall(VecRestoreArrayRead(Diseq, &di));
332:     } else {
333:       val = 2.0 * (1 - di[0] + di[1]);
334:     }
335:     PetscCall(VecRestoreArrayRead(Diseq, &di));
336:     PetscCall(MatSetValues(H, 1, &zero, 1, &zero, &val, INSERT_VALUES));
337:     PetscCall(MatSetValues(H, 1, &one, 1, &one, &two, INSERT_VALUES));
338:   }
339:   PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY));
340:   PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY));
341:   if (!user->noeqflag) {
342:     PetscCall(VecScatterDestroy(&Descat));
343:     PetscCall(VecDestroy(&Deseq));
344:   }
345:   PetscCall(VecScatterDestroy(&Discat));
346:   PetscCall(VecDestroy(&Diseq));
347:   PetscFunctionReturn(PETSC_SUCCESS);
348: }

350: /* Evaluate
351:    h = [ x0^2 - x1;
352:          1 -(x0^2 - x1)]
353: */
354: PetscErrorCode FormInequalityConstraints(Tao tao, Vec X, Vec CI, void *ctx)
355: {
356:   const PetscScalar *x;
357:   PetscScalar        ci;
358:   MPI_Comm           comm;
359:   PetscMPIInt        rank;
360:   AppCtx            *user = (AppCtx *)ctx;
361:   Vec                Xseq = user->Xseq;
362:   VecScatter         scat = user->scat;

364:   PetscFunctionBegin;
365:   PetscCall(PetscObjectGetComm((PetscObject)tao, &comm));
366:   PetscCallMPI(MPI_Comm_rank(comm, &rank));

368:   PetscCall(VecScatterBegin(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD));
369:   PetscCall(VecScatterEnd(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD));

371:   if (rank == 0) {
372:     PetscCall(VecGetArrayRead(Xseq, &x));
373:     ci = x[0] * x[0] - x[1];
374:     PetscCall(VecSetValue(CI, 0, ci, INSERT_VALUES));
375:     ci = -x[0] * x[0] + x[1] + 1.0;
376:     PetscCall(VecSetValue(CI, 1, ci, INSERT_VALUES));
377:     PetscCall(VecRestoreArrayRead(Xseq, &x));
378:   }
379:   PetscCall(VecAssemblyBegin(CI));
380:   PetscCall(VecAssemblyEnd(CI));
381:   PetscFunctionReturn(PETSC_SUCCESS);
382: }

384: /* Evaluate
385:    g = [ x0^2 + x1 - 2]
386: */
387: PetscErrorCode FormEqualityConstraints(Tao tao, Vec X, Vec CE, void *ctx)
388: {
389:   const PetscScalar *x;
390:   PetscScalar        ce;
391:   MPI_Comm           comm;
392:   PetscMPIInt        rank;
393:   AppCtx            *user = (AppCtx *)ctx;
394:   Vec                Xseq = user->Xseq;
395:   VecScatter         scat = user->scat;

397:   PetscFunctionBegin;
398:   PetscCall(PetscObjectGetComm((PetscObject)tao, &comm));
399:   PetscCallMPI(MPI_Comm_rank(comm, &rank));

401:   PetscCall(VecScatterBegin(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD));
402:   PetscCall(VecScatterEnd(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD));

404:   if (rank == 0) {
405:     PetscCall(VecGetArrayRead(Xseq, &x));
406:     ce = x[0] * x[0] + x[1] - 2.0;
407:     PetscCall(VecSetValue(CE, 0, ce, INSERT_VALUES));
408:     PetscCall(VecRestoreArrayRead(Xseq, &x));
409:   }
410:   PetscCall(VecAssemblyBegin(CE));
411:   PetscCall(VecAssemblyEnd(CE));
412:   PetscFunctionReturn(PETSC_SUCCESS);
413: }

415: /*
416:   grad h = [  2*x0, -1;
417:              -2*x0,  1]
418: */
419: PetscErrorCode FormInequalityJacobian(Tao tao, Vec X, Mat JI, Mat JIpre, void *ctx)
420: {
421:   AppCtx            *user = (AppCtx *)ctx;
422:   PetscInt           zero = 0, one = 1, cols[2];
423:   PetscScalar        vals[2];
424:   const PetscScalar *x;
425:   Vec                Xseq = user->Xseq;
426:   VecScatter         scat = user->scat;
427:   MPI_Comm           comm;
428:   PetscMPIInt        rank;

430:   PetscFunctionBegin;
431:   PetscCall(PetscObjectGetComm((PetscObject)tao, &comm));
432:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
433:   PetscCall(VecScatterBegin(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD));
434:   PetscCall(VecScatterEnd(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD));

436:   PetscCall(VecGetArrayRead(Xseq, &x));
437:   if (rank == 0) {
438:     cols[0] = 0;
439:     cols[1] = 1;
440:     vals[0] = 2 * x[0];
441:     vals[1] = -1.0;
442:     PetscCall(MatSetValues(JI, 1, &zero, 2, cols, vals, INSERT_VALUES));
443:     vals[0] = -2 * x[0];
444:     vals[1] = 1.0;
445:     PetscCall(MatSetValues(JI, 1, &one, 2, cols, vals, INSERT_VALUES));
446:   }
447:   PetscCall(VecRestoreArrayRead(Xseq, &x));
448:   PetscCall(MatAssemblyBegin(JI, MAT_FINAL_ASSEMBLY));
449:   PetscCall(MatAssemblyEnd(JI, MAT_FINAL_ASSEMBLY));
450:   PetscFunctionReturn(PETSC_SUCCESS);
451: }

453: /*
454:   grad g = [2*x0, 1.0]
455: */
456: PetscErrorCode FormEqualityJacobian(Tao tao, Vec X, Mat JE, Mat JEpre, void *ctx)
457: {
458:   PetscInt           zero = 0, cols[2];
459:   PetscScalar        vals[2];
460:   const PetscScalar *x;
461:   PetscMPIInt        rank;
462:   MPI_Comm           comm;

464:   PetscFunctionBegin;
465:   PetscCall(PetscObjectGetComm((PetscObject)tao, &comm));
466:   PetscCallMPI(MPI_Comm_rank(comm, &rank));

468:   if (rank == 0) {
469:     PetscCall(VecGetArrayRead(X, &x));
470:     cols[0] = 0;
471:     cols[1] = 1;
472:     vals[0] = 2 * x[0];
473:     vals[1] = 1.0;
474:     PetscCall(MatSetValues(JE, 1, &zero, 2, cols, vals, INSERT_VALUES));
475:     PetscCall(VecRestoreArrayRead(X, &x));
476:   }
477:   PetscCall(MatAssemblyBegin(JE, MAT_FINAL_ASSEMBLY));
478:   PetscCall(MatAssemblyEnd(JE, MAT_FINAL_ASSEMBLY));
479:   PetscFunctionReturn(PETSC_SUCCESS);
480: }

482: /*TEST

484:    build:
485:       requires: !complex !defined(PETSC_USE_CXX)

487:    test:
488:       args: -tao_converged_reason -tao_gatol 1.e-6 -tao_type pdipm -tao_pdipm_kkt_shift_pd
489:       requires: mumps !single
490:       filter: sed  -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"

492:    test:
493:       suffix: pdipm_2
494:       requires: mumps
495:       nsize: 2
496:       args: -tao_converged_reason -tao_gatol 1.e-6 -tao_type pdipm
497:       filter: sed  -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"

499:    test:
500:       suffix: 2
501:       args: -tao_converged_reason
502:       filter: sed  -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"

504:    test:
505:       suffix: 3
506:       args: -tao_converged_reason -no_eq
507:       filter: sed  -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"

509:    test:
510:       suffix: 4
511:       args: -tao_converged_reason -tao_almm_type classic
512:       requires: !single
513:       filter: sed  -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"

515:    test:
516:       suffix: 5
517:       args: -tao_converged_reason -tao_almm_type classic -no_eq
518:       requires: !single
519:       filter: sed  -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"

521:    test:
522:       suffix: 6
523:       args: -tao_converged_reason -tao_almm_subsolver_tao_type bqnktr
524:       filter: sed  -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"

526:    test:
527:       suffix: 7
528:       args: -tao_converged_reason -tao_almm_subsolver_tao_type bncg
529:       filter: sed  -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"

531:    test:
532:       suffix: 8
533:       nsize: 2
534:       args: -tao_converged_reason
535:       filter: sed  -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"

537:    test:
538:       suffix: 9
539:       nsize: 2
540:       args: -tao_converged_reason -vec_type cuda -mat_type aijcusparse
541:       requires: cuda
542:       filter: sed  -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"

544:    test:
545:       suffix: 10
546:       args: -tao_converged_reason -tao_almm_type classic -no_eq -no_bound
547:       requires: !single
548:       filter: sed  -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"

550:    test:
551:       suffix: 11
552:       args: -tao_converged_reason -tao_almm_type classic -no_bound
553:       requires: !single
554:       filter: sed  -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"

556:    test:
557:       suffix: 12
558:       args: -tao_converged_reason -tao_almm_type phr -no_eq -no_bound
559:       requires: !single
560:       filter: sed  -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"

562:    test:
563:       suffix: 13
564:       args: -tao_converged_reason -tao_almm_type phr -no_bound
565:       requires: !single
566:       filter: sed  -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"

568: TEST*/