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:   -init_view         : view initial object setup\n\
 25:   -snes_fd           : snes with finite difference Jacobian (needed for pdipm)\n\
 26:   -snes_compare_explicit : compare user Jacobian with finite difference Jacobian \n\
 27:   -tao_monitor_constraint_norm : convergence monitor with constraint norm \n\
 28:   -tao_monitor_solution : view exact solution at each iteration\n\
 29:   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";

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

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

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

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

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

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

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

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

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

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

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

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

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

170:   /* Tell user the correct solution, not an error checking */
171:   if (!user->noeqflag) {
172:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Expected solution: f(1, 1) = -2\n"));
173:   } else {
174:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Expected solution (-no_eq): f(1.73205, 2) = -7.3923\n"));
175:   }

177:   /* create vector x and set initial values */
178:   user->n = 2; /* global length */
179:   nloc    = (size == 1) ? user->n : 1;
180:   PetscCall(VecCreate(PETSC_COMM_WORLD, &user->x));
181:   PetscCall(VecSetSizes(user->x, nloc, user->n));
182:   PetscCall(VecSetFromOptions(user->x));
183:   PetscCall(VecSet(user->x, 0.0));

185:   /* create and set lower and upper bound vectors */
186:   PetscCall(VecDuplicate(user->x, &user->xl));
187:   PetscCall(VecDuplicate(user->x, &user->xu));
188:   PetscCall(VecSet(user->xl, -1.0));
189:   PetscCall(VecSet(user->xu, 2.0));

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

194:   user->ne = 1;
195:   user->ni = 2;
196:   neloc    = (rank == 0) ? user->ne : 0;
197:   niloc    = (size == 1) ? user->ni : 1;

199:   if (!user->noeqflag) {
200:     PetscCall(VecCreate(PETSC_COMM_WORLD, &user->ce)); /* a 1x1 vec for equality constraints */
201:     PetscCall(VecSetSizes(user->ce, neloc, user->ne));
202:     PetscCall(VecSetFromOptions(user->ce));
203:     PetscCall(VecSetUp(user->ce));
204:   }

206:   PetscCall(VecCreate(PETSC_COMM_WORLD, &user->ci)); /* a 2x1 vec for inequality constraints */
207:   PetscCall(VecSetSizes(user->ci, niloc, user->ni));
208:   PetscCall(VecSetFromOptions(user->ci));
209:   PetscCall(VecSetUp(user->ci));

211:   /* nexn & nixn matrices for equally and inequalty constraints */
212:   if (!user->noeqflag) {
213:     PetscCall(MatCreate(PETSC_COMM_WORLD, &user->Ae));
214:     PetscCall(MatSetSizes(user->Ae, neloc, nloc, user->ne, user->n));
215:     PetscCall(MatSetFromOptions(user->Ae));
216:     PetscCall(MatSetUp(user->Ae));
217:   }

219:   PetscCall(MatCreate(PETSC_COMM_WORLD, &user->Ai));
220:   PetscCall(MatSetSizes(user->Ai, niloc, nloc, user->ni, user->n));
221:   PetscCall(MatSetFromOptions(user->Ai));
222:   PetscCall(MatSetUp(user->Ai));

224:   PetscCall(MatCreate(PETSC_COMM_WORLD, &user->H));
225:   PetscCall(MatSetSizes(user->H, nloc, nloc, user->n, user->n));
226:   PetscCall(MatSetFromOptions(user->H));
227:   PetscCall(MatSetUp(user->H));
228:   PetscFunctionReturn(PETSC_SUCCESS);
229: }

231: PetscErrorCode DestroyProblem(AppCtx *user)
232: {
233:   PetscFunctionBegin;
234:   if (!user->noeqflag) PetscCall(MatDestroy(&user->Ae));
235:   PetscCall(MatDestroy(&user->Ai));
236:   PetscCall(MatDestroy(&user->H));

238:   PetscCall(VecDestroy(&user->x));
239:   if (!user->noeqflag) PetscCall(VecDestroy(&user->ce));
240:   PetscCall(VecDestroy(&user->ci));
241:   PetscCall(VecDestroy(&user->xl));
242:   PetscCall(VecDestroy(&user->xu));
243:   PetscCall(VecDestroy(&user->Xseq));
244:   PetscCall(VecScatterDestroy(&user->scat));
245:   PetscFunctionReturn(PETSC_SUCCESS);
246: }

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

262:   PetscFunctionBegin;
263:   /* f = (x0 - 2)^2 + (x1 - 2)^2 - 2*(x0 + x1) */
264:   PetscCall(PetscObjectGetComm((PetscObject)tao, &comm));
265:   PetscCallMPI(MPI_Comm_rank(comm, &rank));

267:   PetscCall(VecScatterBegin(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD));
268:   PetscCall(VecScatterEnd(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD));

270:   fin = 0.0;
271:   if (rank == 0) {
272:     PetscCall(VecGetArrayRead(Xseq, &x));
273:     fin = (x[0] - 2.0) * (x[0] - 2.0) + (x[1] - 2.0) * (x[1] - 2.0) - 2.0 * (x[0] + x[1]);
274:     PetscCall(VecRestoreArrayRead(Xseq, &x));
275:   }
276:   PetscCallMPI(MPIU_Allreduce(&fin, f, 1, MPIU_REAL, MPIU_SUM, comm));

278:   /* G = 2*X - 6 */
279:   PetscCall(VecSet(G, -6.0));
280:   PetscCall(VecAXPY(G, 2.0, X));
281:   PetscFunctionReturn(PETSC_SUCCESS);
282: }

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

304:   PetscFunctionBegin;
305:   PetscCall(TaoGetDualVariables(tao, &DE, &DI));

307:   PetscCall(PetscObjectGetComm((PetscObject)tao, &comm));
308:   PetscCallMPI(MPI_Comm_rank(comm, &rank));

310:   if (!user->noeqflag) {
311:     PetscCall(VecScatterCreateToZero(DE, &Descat, &Deseq));
312:     PetscCall(VecScatterBegin(Descat, DE, Deseq, INSERT_VALUES, SCATTER_FORWARD));
313:     PetscCall(VecScatterEnd(Descat, DE, Deseq, INSERT_VALUES, SCATTER_FORWARD));
314:   }
315:   PetscCall(VecScatterCreateToZero(DI, &Discat, &Diseq));
316:   PetscCall(VecScatterBegin(Discat, DI, Diseq, INSERT_VALUES, SCATTER_FORWARD));
317:   PetscCall(VecScatterEnd(Discat, DI, Diseq, INSERT_VALUES, SCATTER_FORWARD));

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

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

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

359:   PetscFunctionBegin;
360:   PetscCall(PetscObjectGetComm((PetscObject)tao, &comm));
361:   PetscCallMPI(MPI_Comm_rank(comm, &rank));

363:   PetscCall(VecScatterBegin(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD));
364:   PetscCall(VecScatterEnd(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD));

366:   if (rank == 0) {
367:     PetscCall(VecGetArrayRead(Xseq, &x));
368:     ci = x[0] * x[0] - x[1];
369:     PetscCall(VecSetValue(CI, 0, ci, INSERT_VALUES));
370:     ci = -x[0] * x[0] + x[1] + 1.0;
371:     PetscCall(VecSetValue(CI, 1, ci, INSERT_VALUES));
372:     PetscCall(VecRestoreArrayRead(Xseq, &x));
373:   }
374:   PetscCall(VecAssemblyBegin(CI));
375:   PetscCall(VecAssemblyEnd(CI));
376:   PetscFunctionReturn(PETSC_SUCCESS);
377: }

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

392:   PetscFunctionBegin;
393:   PetscCall(PetscObjectGetComm((PetscObject)tao, &comm));
394:   PetscCallMPI(MPI_Comm_rank(comm, &rank));

396:   PetscCall(VecScatterBegin(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD));
397:   PetscCall(VecScatterEnd(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD));

399:   if (rank == 0) {
400:     PetscCall(VecGetArrayRead(Xseq, &x));
401:     ce = x[0] * x[0] + x[1] - 2.0;
402:     PetscCall(VecSetValue(CE, 0, ce, INSERT_VALUES));
403:     PetscCall(VecRestoreArrayRead(Xseq, &x));
404:   }
405:   PetscCall(VecAssemblyBegin(CE));
406:   PetscCall(VecAssemblyEnd(CE));
407:   PetscFunctionReturn(PETSC_SUCCESS);
408: }

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

425:   PetscFunctionBegin;
426:   PetscCall(PetscObjectGetComm((PetscObject)tao, &comm));
427:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
428:   PetscCall(VecScatterBegin(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD));
429:   PetscCall(VecScatterEnd(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD));

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

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

459:   PetscFunctionBegin;
460:   PetscCall(PetscObjectGetComm((PetscObject)tao, &comm));
461:   PetscCallMPI(MPI_Comm_rank(comm, &rank));

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

477: /*TEST

479:    build:
480:       requires: !complex !defined(PETSC_USE_CXX)

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

487:    test:
488:       suffix: pdipm_2
489:       requires: mumps
490:       nsize: 2
491:       args: -tao_converged_reason -tao_gatol 1.e-6 -tao_type pdipm
492:       filter: sed  -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"

494:    test:
495:       suffix: 2
496:       args: -tao_converged_reason
497:       filter: sed  -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"

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

504:    test:
505:       suffix: 4
506:       args: -tao_converged_reason -tao_almm_type classic
507:       requires: !single
508:       filter: sed  -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"

510:    test:
511:       suffix: 5
512:       args: -tao_converged_reason -tao_almm_type classic -no_eq -tao_almm_subsolver_tao_max_it 100
513:       requires: !single !defined(PETSCTEST_VALGRIND)
514:       filter: sed  -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"

516:    test:
517:       suffix: 6
518:       args: -tao_converged_reason -tao_almm_subsolver_tao_type bqnktr
519:       filter: sed  -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"

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

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

532:    test:
533:       suffix: 9
534:       nsize: 2
535:       args: -tao_converged_reason -vec_type cuda -mat_type aijcusparse
536:       requires: cuda
537:       filter: sed  -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"

539: TEST*/