Actual source code: ex1.c

petsc-3.13.6 2020-09-29
Report Typos and Errors
  1: /* Program usage: mpiexec -n 2 ./ex1 [-help] [all TAO options] */

  3: /* ----------------------------------------------------------------------
  4: min f = (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: ---------------------------------------------------------------------- */

 10:  #include <petsctao.h>

 12: static  char help[]= "Solves constrained optimiztion problem using pdipm.\n\
 13: Input parameters include:\n\
 14:   -tao_type pdipm    : sets Tao solver\n\
 15:   -no_eq             : removes the equaility constraints from the problem\n\
 16:   -snes_fd           : snes with finite difference Jacobian (needed for pdipm)\n\
 17:   -tao_cmonitor      : convergence monitor with constraint norm \n\
 18:   -tao_view_solution : view exact solution at each itteration\n\
 19:   Note: external package superlu_dist is requried to run either for ipm or pdipm. Additionally This is designed for a maximum of 2 processors, the code will error if size > 2.\n\n";

 21: /*
 22:    User-defined Section 1.5 Writing Application Codes with PETSc context - contains data needed by the
 23:    Section 1.5 Writing Application Codes with PETSc-provided call-back routines, FormFunction(),
 24:    FormGradient(), and FormHessian().
 25: */

 27: /*
 28:    x,d in R^n
 29:    f in R
 30:    bin in R^mi
 31:    beq in R^me
 32:    Aeq in R^(me x n)
 33:    Ain in R^(mi x n)
 34:    H in R^(n x n)
 35:    min f=(1/2)*x'*H*x + d'*x
 36:    s.t.  Aeq*x == beq
 37:          Ain*x >= bin
 38: */
 39: typedef struct {
 40:   PetscInt n;  /* Global length of x */
 41:   PetscInt ne; /* Global number of equality constraints */
 42:   PetscInt ni; /* Global number of inequality constraints */
 43:   PetscBool noeqflag;
 44:   Vec      x,xl,xu;
 45:   Vec      ce,ci,bl,bu,Xseq;
 46:   Mat      Ae,Ai,H;
 47:   VecScatter scat;
 48: } AppCtx;

 50: /* -------- User-defined Routines --------- */
 51: PetscErrorCode InitializeProblem(AppCtx *);
 52: PetscErrorCode DestroyProblem(AppCtx *);
 53: PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal *,Vec,void *);
 54: PetscErrorCode FormHessian(Tao,Vec,Mat,Mat, void*);
 55: PetscErrorCode FormInequalityConstraints(Tao,Vec,Vec,void*);
 56: PetscErrorCode FormEqualityConstraints(Tao,Vec,Vec,void*);
 57: PetscErrorCode FormInequalityJacobian(Tao,Vec,Mat,Mat, void*);
 58: PetscErrorCode FormEqualityJacobian(Tao,Vec,Mat,Mat, void*);

 60: PetscErrorCode main(int argc,char **argv)
 61: {
 63:   Tao            tao;
 64:   KSP            ksp;
 65:   PC             pc;
 66:   AppCtx         user;  /* Section 1.5 Writing Application Codes with PETSc context */
 67:   Vec            x;
 68:   PetscMPIInt    rank;

 70:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
 71:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 72:   if (rank>1){
 73:     SETERRQ(PETSC_COMM_SELF,1,"More than 2 processors detected. Example written to use max of 2 processors.\n");
 74:   }

 76:   PetscPrintf(PETSC_COMM_WORLD,"---- Constrained Problem -----\n");
 77:   InitializeProblem(&user); /* sets up problem, function below */
 78:   TaoCreate(PETSC_COMM_WORLD,&tao);
 79:   TaoSetType(tao,TAOPDIPM);
 80:   TaoSetInitialVector(tao,user.x); /* gets solution vector from problem */
 81:   TaoSetVariableBounds(tao,user.xl,user.xu); /* sets lower upper bounds from given solution */
 82:   TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void*)&user);

 84:   if (!user.noeqflag){
 85:     TaoSetEqualityConstraintsRoutine(tao,user.ce,FormEqualityConstraints,(void*)&user);
 86:   }
 87:   TaoSetInequalityConstraintsRoutine(tao,user.ci,FormInequalityConstraints,(void*)&user);

 89:   if (!user.noeqflag){
 90:     TaoSetJacobianEqualityRoutine(tao,user.Ae,user.Ae,FormEqualityJacobian,(void*)&user); /* equality jacobian */
 91:   }
 92:   TaoSetJacobianInequalityRoutine(tao,user.Ai,user.Ai,FormInequalityJacobian,(void*)&user); /* inequality jacobian */
 93:   TaoSetHessianRoutine(tao,user.H,user.H,FormHessian,(void*)&user);
 94:   TaoSetTolerances(tao,1.e-6,1.e-6,1.e-6);
 95:   TaoSetConstraintTolerances(tao,1.e-6,1.e-6);

 97:   TaoGetKSP(tao,&ksp);
 98:   KSPGetPC(ksp,&pc);
 99:   PCSetType(pc,PCLU);
100:   /*
101:       This algorithm produces matrices with zeros along the diagonal therefore we use
102:     SuperLU_DIST which provides shift to the diagonal
103:   */
104:   PCFactorSetMatSolverType(pc,MATSOLVERSUPERLU_DIST);  /* requires superlu_dist to solve pdipm */
105:   KSPSetType(ksp,KSPPREONLY);
106:   KSPSetFromOptions(ksp);

108:   TaoSetFromOptions(tao);

110:   TaoSolve(tao);
111:   TaoGetSolutionVector(tao,&x);
112:   VecView(x,PETSC_VIEWER_STDOUT_WORLD);

114:   /* Free objects */
115:   DestroyProblem(&user);
116:   TaoDestroy(&tao);
117:   PetscFinalize();
118:   return ierr;
119: }

121: PetscErrorCode InitializeProblem(AppCtx *user)
122: {
124:   PetscMPIInt    size;
125:   PetscMPIInt    rank;
126:   PetscInt       nloc,neloc,niloc;

129:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
130:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
131:   user->noeqflag = PETSC_FALSE;
132:   PetscOptionsGetBool(NULL,NULL,"-no_eq",&user->noeqflag,NULL);
133:   if (!user->noeqflag) {
134:     PetscPrintf(PETSC_COMM_WORLD,"Solution should be f(1,1)=-2\n");
135:   }

137:   /* create vector x and set initial values */
138:   user->n = 2; /* global length */
139:   nloc = (rank==0)?user->n:0;
140:   VecCreateMPI(PETSC_COMM_WORLD,nloc,user->n,&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:     neloc = (rank==0)?user->ne:0;

155:   if (!user->noeqflag){
156:     VecCreateMPI(PETSC_COMM_WORLD,neloc,user->ne,&user->ce); /* a 1x1 vec for equality constraints */
157:   }
158:   user->ni = 2;
159:   niloc = (rank==0)?user->ni:0;
160:   VecCreateMPI(PETSC_COMM_WORLD,niloc,user->ni,&user->ci); /* a 2x1 vec for inequality constraints */

162:   /* nexn & nixn matricies for equaly and inequalty constriants */
163:   if (!user->noeqflag){
164:     MatCreate(PETSC_COMM_WORLD,&user->Ae);
165:     MatSetSizes(user->Ae,neloc,nloc,user->ne,user->n);
166:     MatSetUp(user->Ae);
167:     MatSetFromOptions(user->Ae);
168:   }

170:   MatCreate(PETSC_COMM_WORLD,&user->Ai);
171:   MatCreate(PETSC_COMM_WORLD,&user->H);
172:  
173:   MatSetSizes(user->Ai,niloc,nloc,user->ni,user->n);
174:   MatSetSizes(user->H,nloc,nloc,user->n,user->n);

176:   MatSetUp(user->Ai);
177:   MatSetUp(user->H);

179:   MatSetFromOptions(user->Ai);
180:   MatSetFromOptions(user->H);
181:   return(0);
182: }

184: PetscErrorCode DestroyProblem(AppCtx *user)
185: {

189:   if (!user->noeqflag){
190:    MatDestroy(&user->Ae);
191:   }
192:   MatDestroy(&user->Ai);
193:   MatDestroy(&user->H);

195:   VecDestroy(&user->x);
196:   if (!user->noeqflag){
197:     VecDestroy(&user->ce);
198:   }
199:   VecDestroy(&user->ci);
200:   VecDestroy(&user->xl);
201:   VecDestroy(&user->xu);
202:   VecDestroy(&user->Xseq);
203:   VecScatterDestroy(&user->scat);
204:   return(0);
205: }

207: /*
208:   f(X) = (x0 - 2)^2 + (x1 - 2)^2 - 2*(x0 + x1)
209:   G(X) = fx = [2.0*(x[0]-2.0) - 2.0; 2.0*(x[1]-2.0) - 2.0]
210: */
211: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ctx)
212: {
213:   PetscScalar       g;
214:   const PetscScalar *x;
215:   MPI_Comm          comm;
216:   PetscMPIInt       rank;
217:   PetscErrorCode    ierr;
218:   PetscReal         fin;
219:   AppCtx            *user=(AppCtx*)ctx;
220:   Vec               Xseq=user->Xseq;
221:   VecScatter        scat=user->scat;

224:   PetscObjectGetComm((PetscObject)tao,&comm);
225:   MPI_Comm_rank(comm,&rank);

227:   VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
228:   VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);

230:   fin = 0.0;
231:   if (!rank) {
232:     VecGetArrayRead(Xseq,&x);
233:     fin = (x[0]-2.0)*(x[0]-2.0) + (x[1]-2.0)*(x[1]-2.0) - 2.0*(x[0]+x[1]);
234:     g = 2.0*(x[0]-2.0) - 2.0;
235:     VecSetValue(G,0,g,INSERT_VALUES);
236:     g = 2.0*(x[1]-2.0) - 2.0;
237:     VecSetValue(G,1,g,INSERT_VALUES);
238:     VecRestoreArrayRead(Xseq,&x);
239:   }
240:   MPI_Allreduce(&fin,f,1,MPIU_REAL,MPIU_SUM,comm);
241:   VecAssemblyBegin(G);
242:   VecAssemblyEnd(G);
243:   return(0);
244: }

246: PetscErrorCode FormHessian(Tao tao, Vec x,Mat H, Mat Hpre, void *ctx)
247: {
248:   AppCtx            *user=(AppCtx*)ctx;
249:   Vec               DE,DI;
250:   const PetscScalar *de, *di;
251:   PetscInt          zero=0,one=1;
252:   PetscScalar       two=2.0;
253:   PetscScalar       val=0.0;
254:   Vec               Deseq,Diseq=user->Xseq;
255:   VecScatter        Descat,Discat=user->scat;
256:   PetscMPIInt       rank;
257:   MPI_Comm          comm;
258:   PetscErrorCode    ierr;

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

263:   PetscObjectGetComm((PetscObject)tao,&comm);
264:   MPI_Comm_rank(comm,&rank);

266:   if (!user->noeqflag){
267:    VecScatterCreateToZero(DE,&Descat,&Deseq);
268:    VecScatterBegin(Descat,DE,Deseq,INSERT_VALUES,SCATTER_FORWARD);
269:    VecScatterEnd(Descat,DE,Deseq,INSERT_VALUES,SCATTER_FORWARD);
270:   }

272:   VecScatterBegin(Discat,DI,Diseq,INSERT_VALUES,SCATTER_FORWARD);
273:   VecScatterEnd(Discat,DI,Diseq,INSERT_VALUES,SCATTER_FORWARD);
274:   

276:   if (!rank){
277:     if (!user->noeqflag){
278:       VecGetArrayRead(Deseq,&de);  /* places equality constraint dual into array */
279:     }

281:     VecGetArrayRead(Diseq,&di);  /* places inequality constraint dual into array */
282:     
283:     if (!user->noeqflag){
284:       val = 2.0 * (1 + de[0] + di[0] - di[1]);
285:       VecRestoreArrayRead(Deseq,&de);
286:     }else{
287:       val = 2.0 * (1 + di[0] - di[1]);
288:     }
289:     VecRestoreArrayRead(Diseq,&di);
290:     MatSetValues(H,1,&zero,1,&zero,&val,INSERT_VALUES);
291:     MatSetValues(H,1,&one,1,&one,&two,INSERT_VALUES);
292:   }

294:   MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY);
295:   MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY);
296:   if (!user->noeqflag){
297:     VecScatterDestroy(&Descat);
298:     VecDestroy(&Deseq);
299:   }
300:   return(0);
301: }

303: PetscErrorCode FormInequalityConstraints(Tao tao,Vec X,Vec CI,void *ctx)
304: {
305:   const PetscScalar *x;
306:   PetscScalar       ci;
307:   PetscErrorCode    ierr;
308:   MPI_Comm          comm;
309:   PetscMPIInt       rank;
310:   AppCtx            *user=(AppCtx*)ctx;
311:   Vec               Xseq=user->Xseq;
312:   VecScatter        scat=user->scat;

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

318:   VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
319:   VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);

321:   if (!rank) {
322:     VecGetArrayRead(Xseq,&x);
323:     ci = x[0]*x[0] - x[1];
324:     VecSetValue(CI,0,ci,INSERT_VALUES);
325:     ci = -x[0]*x[0] + x[1] + 1.0;
326:     VecSetValue(CI,1,ci,INSERT_VALUES);
327:     VecRestoreArrayRead(Xseq,&x);
328:   }
329:   VecAssemblyBegin(CI);
330:   VecAssemblyEnd(CI);
331:   return(0);
332: }

334: PetscErrorCode FormEqualityConstraints(Tao tao,Vec X,Vec CE,void *ctx)
335: {
336:   const PetscScalar *x;
337:   PetscScalar       ce;
338:   PetscErrorCode    ierr;
339:   MPI_Comm          comm;
340:   PetscMPIInt       rank;
341:   AppCtx            *user=(AppCtx*)ctx;
342:   Vec               Xseq=user->Xseq;
343:   VecScatter        scat=user->scat;

346:   PetscObjectGetComm((PetscObject)tao,&comm);
347:   MPI_Comm_rank(comm,&rank);

349:   VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
350:   VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);

352:   if (!rank) {
353:     VecGetArrayRead(Xseq,&x);
354:     ce = x[0]*x[0] + x[1] - 2.0;
355:     VecSetValue(CE,0,ce,INSERT_VALUES);
356:     VecRestoreArrayRead(Xseq,&x);
357:   }
358:   VecAssemblyBegin(CE);
359:   VecAssemblyEnd(CE);
360:   return(0);
361: }

363: PetscErrorCode FormInequalityJacobian(Tao tao, Vec X, Mat JI, Mat JIpre,  void *ctx)
364: {
365:   AppCtx            *user=(AppCtx*)ctx;
366:   PetscInt          cols[2],min,max,i;
367:   PetscScalar       vals[2];
368:   const PetscScalar *x;
369:   PetscErrorCode    ierr;
370:   Vec               Xseq=user->Xseq;
371:   VecScatter        scat=user->scat;
372:   MPI_Comm          comm;
373:   PetscMPIInt       rank;

376:   PetscObjectGetComm((PetscObject)tao,&comm);
377:   MPI_Comm_rank(comm,&rank);
378:   VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
379:   VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);

381:   VecGetArrayRead(Xseq,&x);
382:   MatGetOwnershipRange(JI,&min,&max);

384:   cols[0] = 0; cols[1] = 1;
385:   for (i=min;i<max;i++) {
386:     if (i==0){
387:       vals[0] = +2*x[0]; vals[1] = -1.0;
388:       MatSetValues(JI,1,&i,2,cols,vals,INSERT_VALUES);
389:     }
390:     if (i==1) {
391:       vals[0] = -2*x[0]; vals[1] = +1.0;
392:       MatSetValues(JI,1,&i,2,cols,vals,INSERT_VALUES);
393:     }
394:   }
395:   VecRestoreArrayRead(Xseq,&x);

397:   MatAssemblyBegin(JI,MAT_FINAL_ASSEMBLY);
398:   MatAssemblyEnd(JI,MAT_FINAL_ASSEMBLY);
399:   return(0);
400: }

402: PetscErrorCode FormEqualityJacobian(Tao tao,Vec X,Mat JE,Mat JEpre,void *ctx)
403: {
404:   PetscInt          rows[2];
405:   PetscScalar       vals[2];
406:   const PetscScalar *x;
407:   PetscMPIInt       rank;
408:   MPI_Comm          comm;
409:   PetscErrorCode    ierr;

412:   PetscObjectGetComm((PetscObject)tao,&comm);
413:   MPI_Comm_rank(comm,&rank);

415:   if (!rank) {
416:     VecGetArrayRead(X,&x);
417:     rows[0] = 0;       rows[1] = 1;
418:     vals[0] = 2*x[0];  vals[1] = 1.0;
419:     MatSetValues(JE,1,rows,2,rows,vals,INSERT_VALUES);
420:     VecRestoreArrayRead(X,&x);
421:   }
422:   MatAssemblyBegin(JE,MAT_FINAL_ASSEMBLY);
423:   MatAssemblyEnd(JE,MAT_FINAL_ASSEMBLY);
424:   return(0);
425: }


428: /*TEST

430:    build:
431:       requires: !complex !define(PETSC_USE_CXX)

433:    test:
434:       requires: superlu_dist
435:       args: -tao_converged_reason

437:    test:
438:       suffix: 2
439:       requires: superlu_dist
440:       nsize: 2
441:       args: -tao_converged_reason

443:    test:
444:       suffix: 3
445:       requires: superlu_dist
446:       args: -tao_converged_reason -no_eq

448:    test:
449:       suffix: 4
450:       requires: superlu_dist
451:       nsize: 2
452:       args: -tao_converged_reason -no_eq

454: TEST*/