Actual source code: ex1.c

petsc-3.14.0 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 application context - contains data needed by the
 23:    application-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;  /* application 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);

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);

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

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

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

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

313:   PetscObjectGetComm((PetscObject)tao,&comm);
314:   MPI_Comm_rank(comm,&rank);

316:   VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
317:   VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);

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

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

344:   PetscObjectGetComm((PetscObject)tao,&comm);
345:   MPI_Comm_rank(comm,&rank);

347:   VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
348:   VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);

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

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

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

379:   VecGetArrayRead(Xseq,&x);
380:   MatGetOwnershipRange(JI,&min,&max);

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

395:   MatAssemblyBegin(JI,MAT_FINAL_ASSEMBLY);
396:   MatAssemblyEnd(JI,MAT_FINAL_ASSEMBLY);
397:   return(0);
398: }

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

410:   PetscObjectGetComm((PetscObject)tao,&comm);
411:   MPI_Comm_rank(comm,&rank);

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


426: /*TEST

428:    build:
429:       requires: !complex !define(PETSC_USE_CXX)

431:    test:
432:       requires: superlu_dist
433:       args: -tao_converged_reason

435:    test:
436:       suffix: 2
437:       requires: superlu_dist
438:       nsize: 2
439:       args: -tao_converged_reason

441:    test:
442:       suffix: 3
443:       requires: superlu_dist
444:       args: -tao_converged_reason -no_eq

446:    test:
447:       suffix: 4
448:       requires: superlu_dist
449:       nsize: 2
450:       args: -tao_converged_reason -no_eq

452: TEST*/