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: {
 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;

 68:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
 69:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 70:   if (size>2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_WRONG_MPI_SIZE,"More than 2 processors detected. Example written to use max of 2 processors.\n");

 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:   TaoSetInitialVector(tao,user.x); /* gets solution vector from problem */
 77:   TaoSetVariableBounds(tao,user.xl,user.xu); /* sets lower upper bounds from given solution */
 78:   TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void*)&user);

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

 91:   TaoGetKSP(tao,&ksp);
 92:   KSPGetPC(ksp,&pc);
 93:   PCSetType(pc,PCCHOLESKY);
 94:   /*
 95:       This algorithm produces matrices with zeros along the diagonal therefore we use
 96:     MUMPS which provides solver for indefinite matrices
 97:   */
 98: #if defined(PETSC_HAVE_MUMPS)
 99:   PCFactorSetMatSolverType(pc,MATSOLVERMUMPS);  /* requires mumps to solve pdipm */
100: #else
101:   if (size > 1) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Requires an external package that supports parallel PCCHOLESKY, e.g., MUMPS.");
102: #endif
103:   KSPSetType(ksp,KSPPREONLY);
104:   KSPSetFromOptions(ksp);

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

145:   TaoSolve(tao);
146:   TaoGetSolutionVector(tao,&x);
147:   VecView(x,PETSC_VIEWER_STDOUT_WORLD);

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

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

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

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

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

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

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

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

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

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

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

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

223:   MatCreate(PETSC_COMM_WORLD,&user->H);
224:   MatSetSizes(user->H,nloc,nloc,user->n,user->n);
225:   MatSetFromOptions(user->H);
226:   MatSetUp(user->H);
227:   return(0);
228: }

230: PetscErrorCode DestroyProblem(AppCtx *user)
231: {

235:   if (!user->noeqflag) {
236:     MatDestroy(&user->Ae);
237:   }
238:   MatDestroy(&user->Ai);
239:   MatDestroy(&user->H);

241:   VecDestroy(&user->x);
242:   if (!user->noeqflag) {
243:     VecDestroy(&user->ce);
244:   }
245:   VecDestroy(&user->ci);
246:   VecDestroy(&user->xl);
247:   VecDestroy(&user->xu);
248:   VecDestroy(&user->Xseq);
249:   VecScatterDestroy(&user->scat);
250:   return(0);
251: }

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

271:   PetscObjectGetComm((PetscObject)tao,&comm);
272:   MPI_Comm_rank(comm,&rank);

274:   VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
275:   VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);

277:   fin = 0.0;
278:   if (rank == 0) {
279:     VecGetArrayRead(Xseq,&x);
280:     fin = (x[0]-2.0)*(x[0]-2.0) + (x[1]-2.0)*(x[1]-2.0) - 2.0*(x[0]+x[1]);
281:     g = 2.0*(x[0]-2.0) - 2.0;
282:     VecSetValue(G,0,g,INSERT_VALUES);
283:     g = 2.0*(x[1]-2.0) - 2.0;
284:     VecSetValue(G,1,g,INSERT_VALUES);
285:     VecRestoreArrayRead(Xseq,&x);
286:   }
287:   MPI_Allreduce(&fin,f,1,MPIU_REAL,MPIU_SUM,comm);
288:   VecAssemblyBegin(G);
289:   VecAssemblyEnd(G);
290:   return(0);
291: }

293: /* Evaluate
294:    H = fxx + grad (grad g^T*DI) - grad (grad h^T*DE)]
295:      = [ 2*(1+de[0]-di[0]+di[1]), 0;
296:                    0,             2]
297: */
298: PetscErrorCode FormHessian(Tao tao, Vec x,Mat H, Mat Hpre, void *ctx)
299: {
300:   AppCtx            *user=(AppCtx*)ctx;
301:   Vec               DE,DI;
302:   const PetscScalar *de,*di;
303:   PetscInt          zero=0,one=1;
304:   PetscScalar       two=2.0;
305:   PetscScalar       val=0.0;
306:   Vec               Deseq,Diseq;
307:   VecScatter        Descat,Discat;
308:   PetscMPIInt       rank;
309:   MPI_Comm          comm;
310:   PetscErrorCode    ierr;

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

315:   PetscObjectGetComm((PetscObject)tao,&comm);
316:   MPI_Comm_rank(comm,&rank);

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

327:   if (rank == 0) {
328:     if (!user->noeqflag) {
329:       VecGetArrayRead(Deseq,&de);  /* places equality constraint dual into array */
330:     }
331:     VecGetArrayRead(Diseq,&di);  /* places inequality constraint dual into array */

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

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

371:   PetscObjectGetComm((PetscObject)tao,&comm);
372:   MPI_Comm_rank(comm,&rank);

374:   VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
375:   VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);

377:   if (rank == 0) {
378:     VecGetArrayRead(Xseq,&x);
379:     ci = x[0]*x[0] - x[1];
380:     VecSetValue(CI,0,ci,INSERT_VALUES);
381:     ci = -x[0]*x[0] + x[1] + 1.0;
382:     VecSetValue(CI,1,ci,INSERT_VALUES);
383:     VecRestoreArrayRead(Xseq,&x);
384:   }
385:   VecAssemblyBegin(CI);
386:   VecAssemblyEnd(CI);
387:   return(0);
388: }

390: /* Evaluate
391:    g = [ x0^2 + x1 - 2]
392: */
393: PetscErrorCode FormEqualityConstraints(Tao tao,Vec X,Vec CE,void *ctx)
394: {
395:   const PetscScalar *x;
396:   PetscScalar       ce;
397:   PetscErrorCode    ierr;
398:   MPI_Comm          comm;
399:   PetscMPIInt       rank;
400:   AppCtx            *user=(AppCtx*)ctx;
401:   Vec               Xseq=user->Xseq;
402:   VecScatter        scat=user->scat;

405:   PetscObjectGetComm((PetscObject)tao,&comm);
406:   MPI_Comm_rank(comm,&rank);

408:   VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
409:   VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);

411:   if (rank == 0) {
412:     VecGetArrayRead(Xseq,&x);
413:     ce = x[0]*x[0] + x[1] - 2.0;
414:     VecSetValue(CE,0,ce,INSERT_VALUES);
415:     VecRestoreArrayRead(Xseq,&x);
416:   }
417:   VecAssemblyBegin(CE);
418:   VecAssemblyEnd(CE);
419:   return(0);
420: }

422: /*
423:   grad h = [  2*x0, -1;
424:              -2*x0,  1]
425: */
426: PetscErrorCode FormInequalityJacobian(Tao tao, Vec X, Mat JI, Mat JIpre,  void *ctx)
427: {
428:   AppCtx            *user=(AppCtx*)ctx;
429:   PetscInt          zero=0,one=1,cols[2];
430:   PetscScalar       vals[2];
431:   const PetscScalar *x;
432:   PetscErrorCode    ierr;
433:   Vec               Xseq=user->Xseq;
434:   VecScatter        scat=user->scat;
435:   MPI_Comm          comm;
436:   PetscMPIInt       rank;

439:   PetscObjectGetComm((PetscObject)tao,&comm);
440:   MPI_Comm_rank(comm,&rank);
441:   VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
442:   VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);

444:   VecGetArrayRead(Xseq,&x);
445:   if (!rank) {
446:     cols[0] = 0; cols[1] = 1;
447:     vals[0] = 2*x[0]; vals[1] = -1.0;
448:     MatSetValues(JI,1,&zero,2,cols,vals,INSERT_VALUES);
449:     vals[0] = -2*x[0]; vals[1] = 1.0;
450:     MatSetValues(JI,1,&one,2,cols,vals,INSERT_VALUES);
451:   }
452:   VecRestoreArrayRead(Xseq,&x);
453:   MatAssemblyBegin(JI,MAT_FINAL_ASSEMBLY);
454:   MatAssemblyEnd(JI,MAT_FINAL_ASSEMBLY);
455:   return(0);
456: }

458: /*
459:   grad g = [2*x0
460:              1.0 ]
461: */
462: PetscErrorCode FormEqualityJacobian(Tao tao,Vec X,Mat JE,Mat JEpre,void *ctx)
463: {
464:   PetscInt          zero=0,cols[2];
465:   PetscScalar       vals[2];
466:   const PetscScalar *x;
467:   PetscMPIInt       rank;
468:   MPI_Comm          comm;
469:   PetscErrorCode    ierr;

472:   PetscObjectGetComm((PetscObject)tao,&comm);
473:   MPI_Comm_rank(comm,&rank);

475:   if (rank == 0) {
476:     VecGetArrayRead(X,&x);
477:     cols[0] = 0;       cols[1] = 1;
478:     vals[0] = 2*x[0];  vals[1] = 1.0;
479:     MatSetValues(JE,1,&zero,2,cols,vals,INSERT_VALUES);
480:     VecRestoreArrayRead(X,&x);
481:   }
482:   MatAssemblyBegin(JE,MAT_FINAL_ASSEMBLY);
483:   MatAssemblyEnd(JE,MAT_FINAL_ASSEMBLY);
484:   return(0);
485: }

487: /*TEST

489:    build:
490:       requires: !complex !defined(PETSC_USE_CXX)

492:    test:
493:       args: -tao_converged_reason -tao_pdipm_kkt_shift_pd
494:       requires: mumps

496:    test:
497:       suffix: 2
498:       nsize: 2
499:       args: -tao_converged_reason -tao_pdipm_kkt_shift_pd
500:       requires: mumps

502:    test:
503:       suffix: 3
504:       args: -tao_converged_reason -no_eq
505:       requires: mumps

507:    test:
508:       suffix: 4
509:       nsize: 2
510:       args: -tao_converged_reason -no_eq
511:       requires: mumps

513:    test:
514:       suffix: 5
515:       args: -tao_cmonitor -tao_type almm
516:       requires: mumps

518:    test:
519:       suffix: 6
520:       args: -tao_cmonitor -tao_type almm -tao_almm_type phr
521:       requires: mumps

523:    test:
524:       suffix: 7
525:       nsize: 2
526:       requires: mumps
527:       args: -tao_cmonitor -tao_type almm

529:    test:
530:       suffix: 8
531:       nsize: 2
532:       requires: cuda mumps
533:       args: -tao_cmonitor -tao_type almm -vec_type cuda -mat_type aijcusparse

535:    test:
536:       suffix: 9
537:       nsize: 2
538:       args: -tao_cmonitor -tao_type almm -no_eq
539:       requires: mumps

541:    test:
542:       suffix: 10
543:       nsize: 2
544:       args: -tao_cmonitor -tao_type almm -tao_almm_type phr -no_eq
545:       requires: mumps

547: TEST*/