Actual source code: ex1.c

petsc-main 2021-04-20
Report Typos and Errors
  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:   -snes_fd           : snes with finite difference Jacobian (needed for pdipm)\n\
 25:   -snes_compare_explicit : compare user Jacobian with finite difference Jacobian \n\
 26:   -tao_cmonitor      : convergence monitor with constraint norm \n\
 27:   -tao_view_solution : view exact solution at each itteration\n\
 28:   Note: external package MUMPS is required to run pdipm. This is designed for a maximum of 2 processors, the code will error if size > 2.\n";

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

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

 54: PetscErrorCode main(int argc,char **argv)
 55: {
 57:   Tao            tao;
 58:   KSP            ksp;
 59:   PC             pc;
 60:   AppCtx         user;  /* application context */
 61:   Vec            x;
 62:   PetscMPIInt    rank;
 63:   TaoType        type;
 64:   PetscBool      pdipm;

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

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

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

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

100:   TaoSetFromOptions(tao);
101:   TaoGetType(tao,&type);
102:   PetscObjectTypeCompare((PetscObject)tao,TAOPDIPM,&pdipm);
103:   if (pdipm) {
104:     TaoSetHessianRoutine(tao,user.H,user.H,FormHessian,(void*)&user);
105:   }

107:   TaoSolve(tao);
108:   TaoGetSolutionVector(tao,&x);
109:   VecView(x,PETSC_VIEWER_STDOUT_WORLD);

111:   /* Free objects */
112:   DestroyProblem(&user);
113:   TaoDestroy(&tao);
114:   PetscFinalize();
115:   return ierr;
116: }

118: PetscErrorCode InitializeProblem(AppCtx *user)
119: {
121:   PetscMPIInt    size;
122:   PetscMPIInt    rank;
123:   PetscInt       nloc,neloc,niloc;

126:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
127:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
128:   user->noeqflag = PETSC_FALSE;
129:   PetscOptionsGetBool(NULL,NULL,"-no_eq",&user->noeqflag,NULL);

131:   if (!user->noeqflag) {
132:     PetscPrintf(PETSC_COMM_WORLD,"Solution should be f(1,1)=-2\n");
133:   }

135:   /* create vector x and set initial values */
136:   user->n = 2; /* global length */
137:   nloc = (rank==0)?user->n:0;
138:   VecCreate(PETSC_COMM_WORLD,&user->x);
139:   VecSetSizes(user->x,nloc,user->n);
140:   VecSetFromOptions(user->x);
141:   VecSet(user->x,0);

143:   /* create and set lower and upper bound vectors */
144:   VecDuplicate(user->x,&user->xl);
145:   VecDuplicate(user->x,&user->xu);
146:   VecSet(user->xl,-1.0);
147:   VecSet(user->xu,2.0);

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

152:     user->ne = 1;
153:     user->ni = 2;
154:     neloc = (rank==0)?user->ne:0;
155:     niloc = (rank==0)?user->ni:0;

157:   if (!user->noeqflag){
158:     VecCreate(PETSC_COMM_WORLD,&user->ce); /* a 1x1 vec for equality constraints */
159:     VecSetSizes(user->ce,neloc,user->ne);
160:     VecSetFromOptions(user->ce);
161:     VecSetUp(user->ce);
162:   }

164:   VecCreate(PETSC_COMM_WORLD,&user->ci); /* a 2x1 vec for inequality constraints */
165:   VecSetSizes(user->ci,niloc,user->ni);
166:   VecSetFromOptions(user->ci);
167:   VecSetUp(user->ci);

169:   /* nexn & nixn matricies for equaly and inequalty constriants */
170:   if (!user->noeqflag){
171:     MatCreate(PETSC_COMM_WORLD,&user->Ae);
172:     MatSetSizes(user->Ae,neloc,nloc,user->ne,user->n);
173:     MatSetFromOptions(user->Ae);
174:     MatSetUp(user->Ae);
175:   }

177:   MatCreate(PETSC_COMM_WORLD,&user->Ai);
178:   MatSetSizes(user->Ai,niloc,nloc,user->ni,user->n);
179:   MatSetFromOptions(user->Ai);
180:   MatSetUp(user->Ai);

182:   MatCreate(PETSC_COMM_WORLD,&user->H);
183:   MatSetSizes(user->H,nloc,nloc,user->n,user->n);
184:   MatSetFromOptions(user->H);
185:   MatSetUp(user->H);
186:   return(0);
187: }

189: PetscErrorCode DestroyProblem(AppCtx *user)
190: {

194:   if (!user->noeqflag){
195:    MatDestroy(&user->Ae);
196:   }
197:   MatDestroy(&user->Ai);
198:   MatDestroy(&user->H);

200:   VecDestroy(&user->x);
201:   if (!user->noeqflag){
202:     VecDestroy(&user->ce);
203:   }
204:   VecDestroy(&user->ci);
205:   VecDestroy(&user->xl);
206:   VecDestroy(&user->xu);
207:   VecDestroy(&user->Xseq);
208:   VecScatterDestroy(&user->scat);
209:   return(0);
210: }

212: /* Evaluate
213:    f(x) = (x0 - 2)^2 + (x1 - 2)^2 - 2*(x0 + x1)
214:    G = grad f = [2*(x0 - 2) - 2;
215:                  2*(x1 - 2) - 2]
216: */
217: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ctx)
218: {
219:   PetscScalar       g;
220:   const PetscScalar *x;
221:   MPI_Comm          comm;
222:   PetscMPIInt       rank;
223:   PetscErrorCode    ierr;
224:   PetscReal         fin;
225:   AppCtx            *user=(AppCtx*)ctx;
226:   Vec               Xseq=user->Xseq;
227:   VecScatter        scat=user->scat;

230:   PetscObjectGetComm((PetscObject)tao,&comm);
231:   MPI_Comm_rank(comm,&rank);

233:   VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
234:   VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);

236:   fin = 0.0;
237:   if (!rank) {
238:     VecGetArrayRead(Xseq,&x);
239:     fin = (x[0]-2.0)*(x[0]-2.0) + (x[1]-2.0)*(x[1]-2.0) - 2.0*(x[0]+x[1]);
240:     g = 2.0*(x[0]-2.0) - 2.0;
241:     VecSetValue(G,0,g,INSERT_VALUES);
242:     g = 2.0*(x[1]-2.0) - 2.0;
243:     VecSetValue(G,1,g,INSERT_VALUES);
244:     VecRestoreArrayRead(Xseq,&x);
245:   }
246:   MPI_Allreduce(&fin,f,1,MPIU_REAL,MPIU_SUM,comm);
247:   VecAssemblyBegin(G);
248:   VecAssemblyEnd(G);
249:   return(0);
250: }

252: /* Evaluate
253:    H = fxx + grad (grad g^T*DI) - grad (grad h^T*DE)]
254:      = [ 2*(1+de[0]-di[0]+di[1]), 0;
255:                    0,             2]
256: */
257: PetscErrorCode FormHessian(Tao tao, Vec x,Mat H, Mat Hpre, void *ctx)
258: {
259:   AppCtx            *user=(AppCtx*)ctx;
260:   Vec               DE,DI;
261:   const PetscScalar *de,*di;
262:   PetscInt          zero=0,one=1;
263:   PetscScalar       two=2.0;
264:   PetscScalar       val=0.0;
265:   Vec               Deseq,Diseq=user->Xseq;
266:   VecScatter        Descat,Discat=user->scat;
267:   PetscMPIInt       rank;
268:   MPI_Comm          comm;
269:   PetscErrorCode    ierr;

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

274:   PetscObjectGetComm((PetscObject)tao,&comm);
275:   MPI_Comm_rank(comm,&rank);

277:   if (!user->noeqflag){
278:    VecScatterCreateToZero(DE,&Descat,&Deseq);
279:    VecScatterBegin(Descat,DE,Deseq,INSERT_VALUES,SCATTER_FORWARD);
280:    VecScatterEnd(Descat,DE,Deseq,INSERT_VALUES,SCATTER_FORWARD);
281:   }
282:   VecScatterBegin(Discat,DI,Diseq,INSERT_VALUES,SCATTER_FORWARD);
283:   VecScatterEnd(Discat,DI,Diseq,INSERT_VALUES,SCATTER_FORWARD);

285:   if (!rank){
286:     if (!user->noeqflag){
287:       VecGetArrayRead(Deseq,&de);  /* places equality constraint dual into array */
288:     }
289:     VecGetArrayRead(Diseq,&di);  /* places inequality constraint dual into array */

291:     if (!user->noeqflag){
292:       val = 2.0 * (1 + de[0] - di[0] + di[1]);
293:       VecRestoreArrayRead(Deseq,&de);
294:       VecRestoreArrayRead(Diseq,&di);
295:     } else {
296:       val = 2.0 * (1 - di[0] + di[1]);
297:     }
298:     VecRestoreArrayRead(Diseq,&di);
299:     MatSetValues(H,1,&zero,1,&zero,&val,INSERT_VALUES);
300:     MatSetValues(H,1,&one,1,&one,&two,INSERT_VALUES);
301:   }
302:   MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY);
303:   MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY);
304:   if (!user->noeqflag){
305:     VecScatterDestroy(&Descat);
306:     VecDestroy(&Deseq);
307:   }
308:   return(0);
309: }

311: /* Evaluate
312:    h = [ x0^2 - x1;
313:          1 -(x0^2 - x1)]
314: */
315: PetscErrorCode FormInequalityConstraints(Tao tao,Vec X,Vec CI,void *ctx)
316: {
317:   const PetscScalar *x;
318:   PetscScalar       ci;
319:   PetscErrorCode    ierr;
320:   MPI_Comm          comm;
321:   PetscMPIInt       rank;
322:   AppCtx            *user=(AppCtx*)ctx;
323:   Vec               Xseq=user->Xseq;
324:   VecScatter        scat=user->scat;

327:   PetscObjectGetComm((PetscObject)tao,&comm);
328:   MPI_Comm_rank(comm,&rank);

330:   VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
331:   VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);

333:   if (!rank) {
334:     VecGetArrayRead(Xseq,&x);
335:     ci = x[0]*x[0] - x[1];
336:     VecSetValue(CI,0,ci,INSERT_VALUES);
337:     ci = -x[0]*x[0] + x[1] + 1.0;
338:     VecSetValue(CI,1,ci,INSERT_VALUES);
339:     VecRestoreArrayRead(Xseq,&x);
340:   }
341:   VecAssemblyBegin(CI);
342:   VecAssemblyEnd(CI);
343:   return(0);
344: }

346: /* Evaluate
347:    g = [ x0^2 + x1 - 2]
348: */
349: PetscErrorCode FormEqualityConstraints(Tao tao,Vec X,Vec CE,void *ctx)
350: {
351:   const PetscScalar *x;
352:   PetscScalar       ce;
353:   PetscErrorCode    ierr;
354:   MPI_Comm          comm;
355:   PetscMPIInt       rank;
356:   AppCtx            *user=(AppCtx*)ctx;
357:   Vec               Xseq=user->Xseq;
358:   VecScatter        scat=user->scat;

361:   PetscObjectGetComm((PetscObject)tao,&comm);
362:   MPI_Comm_rank(comm,&rank);

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

367:   if (!rank) {
368:     VecGetArrayRead(Xseq,&x);
369:     ce = x[0]*x[0] + x[1] - 2.0;
370:     VecSetValue(CE,0,ce,INSERT_VALUES);
371:     VecRestoreArrayRead(Xseq,&x);
372:   }
373:   VecAssemblyBegin(CE);
374:   VecAssemblyEnd(CE);
375:   return(0);
376: }

378: /*
379:   grad h = [  2*x0, -1;
380:              -2*x0,  1]
381: */
382: PetscErrorCode FormInequalityJacobian(Tao tao, Vec X, Mat JI, Mat JIpre,  void *ctx)
383: {
384:   AppCtx            *user=(AppCtx*)ctx;
385:   PetscInt          cols[2],min,max,i;
386:   PetscScalar       vals[2];
387:   const PetscScalar *x;
388:   PetscErrorCode    ierr;
389:   Vec               Xseq=user->Xseq;
390:   VecScatter        scat=user->scat;
391:   MPI_Comm          comm;
392:   PetscMPIInt       rank;

395:   PetscObjectGetComm((PetscObject)tao,&comm);
396:   MPI_Comm_rank(comm,&rank);
397:   VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
398:   VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);

400:   VecGetArrayRead(Xseq,&x);
401:   MatGetOwnershipRange(JI,&min,&max);

403:   cols[0] = 0; cols[1] = 1;
404:   for (i=min;i<max;i++) {
405:     if (i==0){
406:       vals[0] = 2*x[0]; vals[1] = -1.0;
407:       MatSetValues(JI,1,&i,2,cols,vals,INSERT_VALUES);
408:     }
409:     if (i==1) {
410:       vals[0] = -2*x[0]; vals[1] = 1.0;
411:       MatSetValues(JI,1,&i,2,cols,vals,INSERT_VALUES);
412:     }
413:   }
414:   VecRestoreArrayRead(Xseq,&x);

416:   MatAssemblyBegin(JI,MAT_FINAL_ASSEMBLY);
417:   MatAssemblyEnd(JI,MAT_FINAL_ASSEMBLY);
418:   return(0);
419: }

421: /*
422:   grad g = [2*x0
423:              1.0 ]
424: */
425: PetscErrorCode FormEqualityJacobian(Tao tao,Vec X,Mat JE,Mat JEpre,void *ctx)
426: {
427:   PetscInt          rows[2];
428:   PetscScalar       vals[2];
429:   const PetscScalar *x;
430:   PetscMPIInt       rank;
431:   MPI_Comm          comm;
432:   PetscErrorCode    ierr;

435:   PetscObjectGetComm((PetscObject)tao,&comm);
436:   MPI_Comm_rank(comm,&rank);

438:   if (!rank) {
439:     VecGetArrayRead(X,&x);
440:     rows[0] = 0;       rows[1] = 1;
441:     vals[0] = 2*x[0];  vals[1] = 1.0;
442:     MatSetValues(JE,1,rows,2,rows,vals,INSERT_VALUES);
443:     VecRestoreArrayRead(X,&x);
444:   }
445:   MatAssemblyBegin(JE,MAT_FINAL_ASSEMBLY);
446:   MatAssemblyEnd(JE,MAT_FINAL_ASSEMBLY);
447:   return(0);
448: }


451: /*TEST

453:    build:
454:       requires: !complex !define(PETSC_USE_CXX) mumps

456:    test:
457:       args: -tao_converged_reason -tao_pdipm_kkt_shift_pd

459:    test:
460:       suffix: 2
461:       nsize: 2
462:       args: -tao_converged_reason -tao_pdipm_kkt_shift_pd

464:    test:
465:       suffix: 3
466:       args: -tao_converged_reason -no_eq

468:    test:
469:       suffix: 4
470:       nsize: 2
471:       args: -tao_converged_reason -no_eq

473:    test:
474:       suffix: 5
475:       args: -tao_cmonitor -tao_type almm

477:    test:
478:       suffix: 6
479:       args: -tao_cmonitor -tao_type almm -tao_almm_type phr

481:    test:
482:       suffix: 7
483:       nsize: 2
484:       args: -tao_cmonitor -tao_type almm

486:    test:
487:       suffix: 8
488:       nsize: 2
489:       requires: cuda
490:       args: -tao_cmonitor -tao_type almm -vec_type cuda -mat_type aijcusparse

492:    test:
493:       suffix: 9
494:       nsize: 2
495:       args: -tao_cmonitor -tao_type almm -no_eq

497:    test:
498:       suffix: 10
499:       nsize: 2
500:       args: -tao_cmonitor -tao_type almm -tao_almm_type phr -no_eq

502: TEST*/