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 optimiztion 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_cmonitor      : convergence monitor with constraint norm \n\
 28:   -tao_view_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 FormHessian(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: PetscErrorCode 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;

 68:   PetscInitialize(&argc, &argv, (char *)0, help);
 69:   MPI_Comm_size(PETSC_COMM_WORLD, &size);

 72:   PetscPrintf(PETSC_COMM_WORLD, "---- Constrained Problem -----\n");
 73:   InitializeProblem(&user); /* sets up problem, function below */
 74:   TaoCreate(PETSC_COMM_WORLD, &tao);
 75:   TaoSetType(tao, TAOPDIPM);
 76:   TaoSetSolution(tao, user.x);                 /* gets solution vector from problem */
 77:   TaoSetVariableBounds(tao, user.xl, user.xu); /* sets lower upper bounds from given solution */
 78:   TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, (void *)&user);

 80:   if (!user.noeqflag) TaoSetEqualityConstraintsRoutine(tao, user.ce, FormEqualityConstraints, (void *)&user);
 81:   TaoSetInequalityConstraintsRoutine(tao, user.ci, FormInequalityConstraints, (void *)&user);
 82:   if (!user.noeqflag) { TaoSetJacobianEqualityRoutine(tao, user.Ae, user.Ae, FormEqualityJacobian, (void *)&user); /* equality jacobian */ }
 83:   TaoSetJacobianInequalityRoutine(tao, user.Ai, user.Ai, FormInequalityJacobian, (void *)&user); /* inequality jacobian */
 84:   TaoSetTolerances(tao, 1.e-6, 1.e-6, 1.e-6);
 85:   TaoSetConstraintTolerances(tao, 1.e-6, 1.e-6);

 87:   TaoGetKSP(tao, &ksp);
 88:   KSPGetPC(ksp, &pc);
 89:   PCSetType(pc, PCCHOLESKY);
 90:   /*
 91:       This algorithm produces matrices with zeros along the diagonal therefore we use
 92:     MUMPS which provides solver for indefinite matrices
 93:   */
 94: #if defined(PETSC_HAVE_MUMPS)
 95:   PCFactorSetMatSolverType(pc, MATSOLVERMUMPS); /* requires mumps to solve pdipm */
 96: #else
 98: #endif
 99:   KSPSetType(ksp, KSPPREONLY);
100:   KSPSetFromOptions(ksp);

102:   TaoSetFromOptions(tao);
103:   TaoGetType(tao, &type);
104:   PetscObjectTypeCompare((PetscObject)tao, TAOPDIPM, &pdipm);
105:   if (pdipm) {
106:     TaoSetHessian(tao, user.H, user.H, FormHessian, (void *)&user);
107:     if (user.initview) {
108:       TaoSetUp(tao);
109:       VecDuplicate(user.x, &G);
110:       FormFunctionGradient(tao, user.x, &f, G, (void *)&user);
111:       FormHessian(tao, user.x, user.H, user.H, (void *)&user);
112:       PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_WORLD);
113:       PetscPrintf(PETSC_COMM_WORLD, "\nInitial point X:\n");
114:       VecView(user.x, PETSC_VIEWER_STDOUT_WORLD);
115:       PetscPrintf(PETSC_COMM_WORLD, "\nInitial objective f(x) = %g\n", (double)f);
116:       PetscPrintf(PETSC_COMM_WORLD, "\nInitial gradient and Hessian:\n");
117:       VecView(G, PETSC_VIEWER_STDOUT_WORLD);
118:       MatView(user.H, PETSC_VIEWER_STDOUT_WORLD);
119:       VecDestroy(&G);
120:       FormInequalityJacobian(tao, user.x, user.Ai, user.Ai, (void *)&user);
121:       MatCreateVecs(user.Ai, NULL, &CI);
122:       FormInequalityConstraints(tao, user.x, CI, (void *)&user);
123:       PetscPrintf(PETSC_COMM_WORLD, "\nInitial inequality constraints and Jacobian:\n");
124:       VecView(CI, PETSC_VIEWER_STDOUT_WORLD);
125:       MatView(user.Ai, PETSC_VIEWER_STDOUT_WORLD);
126:       VecDestroy(&CI);
127:       if (!user.noeqflag) {
128:         FormEqualityJacobian(tao, user.x, user.Ae, user.Ae, (void *)&user);
129:         MatCreateVecs(user.Ae, NULL, &CE);
130:         FormEqualityConstraints(tao, user.x, CE, (void *)&user);
131:         PetscPrintf(PETSC_COMM_WORLD, "\nInitial equality constraints and Jacobian:\n");
132:         VecView(CE, PETSC_VIEWER_STDOUT_WORLD);
133:         MatView(user.Ae, PETSC_VIEWER_STDOUT_WORLD);
134:         VecDestroy(&CE);
135:       }
136:       PetscPrintf(PETSC_COMM_WORLD, "\n");
137:       PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_WORLD);
138:     }
139:   }

141:   TaoSolve(tao);
142:   TaoGetSolution(tao, &x);
143:   VecView(x, PETSC_VIEWER_STDOUT_WORLD);

145:   /* Free objects */
146:   DestroyProblem(&user);
147:   TaoDestroy(&tao);
148:   PetscFinalize();
149:   return 0;
150: }

152: PetscErrorCode InitializeProblem(AppCtx *user)
153: {
154:   PetscMPIInt size;
155:   PetscMPIInt rank;
156:   PetscInt    nloc, neloc, niloc;

158:   MPI_Comm_size(PETSC_COMM_WORLD, &size);
159:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
160:   user->noeqflag = PETSC_FALSE;
161:   PetscOptionsGetBool(NULL, NULL, "-no_eq", &user->noeqflag, NULL);
162:   user->initview = PETSC_FALSE;
163:   PetscOptionsGetBool(NULL, NULL, "-init_view", &user->initview, NULL);

165:   if (!user->noeqflag) {
166:     /* Tell user the correct solution, not an error checking */
167:     PetscPrintf(PETSC_COMM_WORLD, "Solution should be f(1,1)=-2\n");
168:   }

170:   /* create vector x and set initial values */
171:   user->n = 2; /* global length */
172:   nloc    = (size == 1) ? user->n : 1;
173:   VecCreate(PETSC_COMM_WORLD, &user->x);
174:   VecSetSizes(user->x, nloc, user->n);
175:   VecSetFromOptions(user->x);
176:   VecSet(user->x, 0);

178:   /* create and set lower and upper bound vectors */
179:   VecDuplicate(user->x, &user->xl);
180:   VecDuplicate(user->x, &user->xu);
181:   VecSet(user->xl, -1.0);
182:   VecSet(user->xu, 2.0);

184:   /* create scater to zero */
185:   VecScatterCreateToZero(user->x, &user->scat, &user->Xseq);

187:   user->ne = 1;
188:   user->ni = 2;
189:   neloc    = (rank == 0) ? user->ne : 0;
190:   niloc    = (size == 1) ? user->ni : 1;

192:   if (!user->noeqflag) {
193:     VecCreate(PETSC_COMM_WORLD, &user->ce); /* a 1x1 vec for equality constraints */
194:     VecSetSizes(user->ce, neloc, user->ne);
195:     VecSetFromOptions(user->ce);
196:     VecSetUp(user->ce);
197:   }

199:   VecCreate(PETSC_COMM_WORLD, &user->ci); /* a 2x1 vec for inequality constraints */
200:   VecSetSizes(user->ci, niloc, user->ni);
201:   VecSetFromOptions(user->ci);
202:   VecSetUp(user->ci);

204:   /* nexn & nixn matricies for equally and inequalty constraints */
205:   if (!user->noeqflag) {
206:     MatCreate(PETSC_COMM_WORLD, &user->Ae);
207:     MatSetSizes(user->Ae, neloc, nloc, user->ne, user->n);
208:     MatSetFromOptions(user->Ae);
209:     MatSetUp(user->Ae);
210:   }

212:   MatCreate(PETSC_COMM_WORLD, &user->Ai);
213:   MatSetSizes(user->Ai, niloc, nloc, user->ni, user->n);
214:   MatSetFromOptions(user->Ai);
215:   MatSetUp(user->Ai);

217:   MatCreate(PETSC_COMM_WORLD, &user->H);
218:   MatSetSizes(user->H, nloc, nloc, user->n, user->n);
219:   MatSetFromOptions(user->H);
220:   MatSetUp(user->H);
221:   return 0;
222: }

224: PetscErrorCode DestroyProblem(AppCtx *user)
225: {
226:   if (!user->noeqflag) MatDestroy(&user->Ae);
227:   MatDestroy(&user->Ai);
228:   MatDestroy(&user->H);

230:   VecDestroy(&user->x);
231:   if (!user->noeqflag) VecDestroy(&user->ce);
232:   VecDestroy(&user->ci);
233:   VecDestroy(&user->xl);
234:   VecDestroy(&user->xu);
235:   VecDestroy(&user->Xseq);
236:   VecScatterDestroy(&user->scat);
237:   return 0;
238: }

240: /* Evaluate
241:    f(x) = (x0 - 2)^2 + (x1 - 2)^2 - 2*(x0 + x1)
242:    G = grad f = [2*(x0 - 2) - 2;
243:                  2*(x1 - 2) - 2]
244: */
245: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ctx)
246: {
247:   PetscScalar        g;
248:   const PetscScalar *x;
249:   MPI_Comm           comm;
250:   PetscMPIInt        rank;
251:   PetscReal          fin;
252:   AppCtx            *user = (AppCtx *)ctx;
253:   Vec                Xseq = user->Xseq;
254:   VecScatter         scat = user->scat;

256:   PetscObjectGetComm((PetscObject)tao, &comm);
257:   MPI_Comm_rank(comm, &rank);

259:   VecScatterBegin(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD);
260:   VecScatterEnd(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD);

262:   fin = 0.0;
263:   if (rank == 0) {
264:     VecGetArrayRead(Xseq, &x);
265:     fin = (x[0] - 2.0) * (x[0] - 2.0) + (x[1] - 2.0) * (x[1] - 2.0) - 2.0 * (x[0] + x[1]);
266:     g   = 2.0 * (x[0] - 2.0) - 2.0;
267:     VecSetValue(G, 0, g, INSERT_VALUES);
268:     g = 2.0 * (x[1] - 2.0) - 2.0;
269:     VecSetValue(G, 1, g, INSERT_VALUES);
270:     VecRestoreArrayRead(Xseq, &x);
271:   }
272:   MPI_Allreduce(&fin, f, 1, MPIU_REAL, MPIU_SUM, comm);
273:   VecAssemblyBegin(G);
274:   VecAssemblyEnd(G);
275:   return 0;
276: }

278: /* Evaluate
279:    H = fxx + grad (grad g^T*DI) - grad (grad h^T*DE)]
280:      = [ 2*(1+de[0]-di[0]+di[1]), 0;
281:                    0,             2]
282: */
283: PetscErrorCode FormHessian(Tao tao, Vec x, Mat H, Mat Hpre, void *ctx)
284: {
285:   AppCtx            *user = (AppCtx *)ctx;
286:   Vec                DE, DI;
287:   const PetscScalar *de, *di;
288:   PetscInt           zero = 0, one = 1;
289:   PetscScalar        two = 2.0;
290:   PetscScalar        val = 0.0;
291:   Vec                Deseq, Diseq;
292:   VecScatter         Descat, Discat;
293:   PetscMPIInt        rank;
294:   MPI_Comm           comm;

296:   TaoGetDualVariables(tao, &DE, &DI);

298:   PetscObjectGetComm((PetscObject)tao, &comm);
299:   MPI_Comm_rank(comm, &rank);

301:   if (!user->noeqflag) {
302:     VecScatterCreateToZero(DE, &Descat, &Deseq);
303:     VecScatterBegin(Descat, DE, Deseq, INSERT_VALUES, SCATTER_FORWARD);
304:     VecScatterEnd(Descat, DE, Deseq, INSERT_VALUES, SCATTER_FORWARD);
305:   }
306:   VecScatterCreateToZero(DI, &Discat, &Diseq);
307:   VecScatterBegin(Discat, DI, Diseq, INSERT_VALUES, SCATTER_FORWARD);
308:   VecScatterEnd(Discat, DI, Diseq, INSERT_VALUES, SCATTER_FORWARD);

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

314:     if (!user->noeqflag) {
315:       val = 2.0 * (1 + de[0] - di[0] + di[1]);
316:       VecRestoreArrayRead(Deseq, &de);
317:       VecRestoreArrayRead(Diseq, &di);
318:     } else {
319:       val = 2.0 * (1 - di[0] + di[1]);
320:     }
321:     VecRestoreArrayRead(Diseq, &di);
322:     MatSetValues(H, 1, &zero, 1, &zero, &val, INSERT_VALUES);
323:     MatSetValues(H, 1, &one, 1, &one, &two, INSERT_VALUES);
324:   }
325:   MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY);
326:   MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY);
327:   if (!user->noeqflag) {
328:     VecScatterDestroy(&Descat);
329:     VecDestroy(&Deseq);
330:   }
331:   VecScatterDestroy(&Discat);
332:   VecDestroy(&Diseq);
333:   return 0;
334: }

336: /* Evaluate
337:    h = [ x0^2 - x1;
338:          1 -(x0^2 - x1)]
339: */
340: PetscErrorCode FormInequalityConstraints(Tao tao, Vec X, Vec CI, void *ctx)
341: {
342:   const PetscScalar *x;
343:   PetscScalar        ci;
344:   MPI_Comm           comm;
345:   PetscMPIInt        rank;
346:   AppCtx            *user = (AppCtx *)ctx;
347:   Vec                Xseq = user->Xseq;
348:   VecScatter         scat = user->scat;

350:   PetscObjectGetComm((PetscObject)tao, &comm);
351:   MPI_Comm_rank(comm, &rank);

353:   VecScatterBegin(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD);
354:   VecScatterEnd(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD);

356:   if (rank == 0) {
357:     VecGetArrayRead(Xseq, &x);
358:     ci = x[0] * x[0] - x[1];
359:     VecSetValue(CI, 0, ci, INSERT_VALUES);
360:     ci = -x[0] * x[0] + x[1] + 1.0;
361:     VecSetValue(CI, 1, ci, INSERT_VALUES);
362:     VecRestoreArrayRead(Xseq, &x);
363:   }
364:   VecAssemblyBegin(CI);
365:   VecAssemblyEnd(CI);
366:   return 0;
367: }

369: /* Evaluate
370:    g = [ x0^2 + x1 - 2]
371: */
372: PetscErrorCode FormEqualityConstraints(Tao tao, Vec X, Vec CE, void *ctx)
373: {
374:   const PetscScalar *x;
375:   PetscScalar        ce;
376:   MPI_Comm           comm;
377:   PetscMPIInt        rank;
378:   AppCtx            *user = (AppCtx *)ctx;
379:   Vec                Xseq = user->Xseq;
380:   VecScatter         scat = user->scat;

382:   PetscObjectGetComm((PetscObject)tao, &comm);
383:   MPI_Comm_rank(comm, &rank);

385:   VecScatterBegin(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD);
386:   VecScatterEnd(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD);

388:   if (rank == 0) {
389:     VecGetArrayRead(Xseq, &x);
390:     ce = x[0] * x[0] + x[1] - 2.0;
391:     VecSetValue(CE, 0, ce, INSERT_VALUES);
392:     VecRestoreArrayRead(Xseq, &x);
393:   }
394:   VecAssemblyBegin(CE);
395:   VecAssemblyEnd(CE);
396:   return 0;
397: }

399: /*
400:   grad h = [  2*x0, -1;
401:              -2*x0,  1]
402: */
403: PetscErrorCode FormInequalityJacobian(Tao tao, Vec X, Mat JI, Mat JIpre, void *ctx)
404: {
405:   AppCtx            *user = (AppCtx *)ctx;
406:   PetscInt           zero = 0, one = 1, cols[2];
407:   PetscScalar        vals[2];
408:   const PetscScalar *x;
409:   Vec                Xseq = user->Xseq;
410:   VecScatter         scat = user->scat;
411:   MPI_Comm           comm;
412:   PetscMPIInt        rank;

414:   PetscObjectGetComm((PetscObject)tao, &comm);
415:   MPI_Comm_rank(comm, &rank);
416:   VecScatterBegin(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD);
417:   VecScatterEnd(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD);

419:   VecGetArrayRead(Xseq, &x);
420:   if (rank == 0) {
421:     cols[0] = 0;
422:     cols[1] = 1;
423:     vals[0] = 2 * x[0];
424:     vals[1] = -1.0;
425:     MatSetValues(JI, 1, &zero, 2, cols, vals, INSERT_VALUES);
426:     vals[0] = -2 * x[0];
427:     vals[1] = 1.0;
428:     MatSetValues(JI, 1, &one, 2, cols, vals, INSERT_VALUES);
429:   }
430:   VecRestoreArrayRead(Xseq, &x);
431:   MatAssemblyBegin(JI, MAT_FINAL_ASSEMBLY);
432:   MatAssemblyEnd(JI, MAT_FINAL_ASSEMBLY);
433:   return 0;
434: }

436: /*
437:   grad g = [2*x0
438:              1.0 ]
439: */
440: PetscErrorCode FormEqualityJacobian(Tao tao, Vec X, Mat JE, Mat JEpre, void *ctx)
441: {
442:   PetscInt           zero = 0, cols[2];
443:   PetscScalar        vals[2];
444:   const PetscScalar *x;
445:   PetscMPIInt        rank;
446:   MPI_Comm           comm;

448:   PetscObjectGetComm((PetscObject)tao, &comm);
449:   MPI_Comm_rank(comm, &rank);

451:   if (rank == 0) {
452:     VecGetArrayRead(X, &x);
453:     cols[0] = 0;
454:     cols[1] = 1;
455:     vals[0] = 2 * x[0];
456:     vals[1] = 1.0;
457:     MatSetValues(JE, 1, &zero, 2, cols, vals, INSERT_VALUES);
458:     VecRestoreArrayRead(X, &x);
459:   }
460:   MatAssemblyBegin(JE, MAT_FINAL_ASSEMBLY);
461:   MatAssemblyEnd(JE, MAT_FINAL_ASSEMBLY);
462:   return 0;
463: }

465: /*TEST

467:    build:
468:       requires: !complex !defined(PETSC_USE_CXX)

470:    test:
471:       args: -tao_converged_reason -tao_pdipm_kkt_shift_pd
472:       requires: mumps

474:    test:
475:       suffix: 2
476:       nsize: 2
477:       args: -tao_converged_reason -tao_pdipm_kkt_shift_pd
478:       requires: mumps

480:    test:
481:       suffix: 3
482:       args: -tao_converged_reason -no_eq
483:       requires: mumps

485:    test:
486:       suffix: 4
487:       nsize: 2
488:       args: -tao_converged_reason -no_eq
489:       requires: mumps

491:    test:
492:       suffix: 5
493:       args: -tao_cmonitor -tao_type almm
494:       requires: mumps

496:    test:
497:       suffix: 6
498:       args: -tao_cmonitor -tao_type almm -tao_almm_type phr
499:       requires: mumps

501:    test:
502:       suffix: 7
503:       nsize: 2
504:       requires: mumps
505:       args: -tao_cmonitor -tao_type almm

507:    test:
508:       suffix: 8
509:       nsize: 2
510:       requires: cuda mumps
511:       args: -tao_cmonitor -tao_type almm -vec_type cuda -mat_type aijcusparse

513:    test:
514:       suffix: 9
515:       nsize: 2
516:       args: -tao_cmonitor -tao_type almm -no_eq
517:       requires: mumps

519:    test:
520:       suffix: 10
521:       nsize: 2
522:       args: -tao_cmonitor -tao_type almm -tao_almm_type phr -no_eq
523:       requires: mumps

525: TEST*/