Actual source code: ex1.c

petsc-master 2020-08-10
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:   -snes_fd           : snes with finite difference Jacobian (needed for pdipm)\n\
 16:   -tao_cmonitor      : convergence monitor with constraint norm \n\
 17:   -tao_view_solution : view exact solution at each itteration\n\
 18:   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";

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

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

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

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

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

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

 83:   TaoSetEqualityConstraintsRoutine(tao,user.ce,FormEqualityConstraints,(void*)&user);
 84:   TaoSetInequalityConstraintsRoutine(tao,user.ci,FormInequalityConstraints,(void*)&user);

 86:   TaoSetJacobianEqualityRoutine(tao,user.Ae,user.Ae,FormEqualityJacobian,(void*)&user); /* equality jacobian */
 87:   TaoSetJacobianInequalityRoutine(tao,user.Ai,user.Ai,FormInequalityJacobian,(void*)&user); /* inequality jacobian */
 88:   TaoSetHessianRoutine(tao,user.H,user.H,FormHessian,(void*)&user);
 89:   TaoSetTolerances(tao,1.e-6,1.e-6,1.e-6);
 90:   TaoSetConstraintTolerances(tao,1.e-6,1.e-6);

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

103:   TaoSetFromOptions(tao);

105:   TaoSolve(tao);
106:   TaoGetSolutionVector(tao,&x);
107:   VecView(x,PETSC_VIEWER_STDOUT_WORLD);

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

116: PetscErrorCode InitializeProblem(AppCtx *user)
117: {
119:   PetscMPIInt    size;
120:   PetscMPIInt    rank;
121:   PetscInt       nloc,neloc,niloc;

124:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
125:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);

127:   /* create vector x and set initial values */
128:   user->n = 2; /* global length */
129:   nloc = (rank==0)?user->n:0;
130:   VecCreateMPI(PETSC_COMM_WORLD,nloc,user->n,&user->x);
131:   VecSet(user->x,0);

133:   /* create and set lower and upper bound vectors */
134:   VecDuplicate(user->x,&user->xl);
135:   VecDuplicate(user->x,&user->xu);
136:   VecSet(user->xl,-1.0);
137:   VecSet(user->xu,2.0);

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

142:   user->ne = 1;
143:   neloc = (rank==0)?user->ne:0;
144:   VecCreateMPI(PETSC_COMM_WORLD,neloc,user->ne,&user->ce); /* a 1x1 vec for equality constraints */

146:   user->ni = 2;
147:   niloc = (rank==0)?user->ni:0;
148:   VecCreateMPI(PETSC_COMM_WORLD,niloc,user->ni,&user->ci); /* a 2x1 vec for inequality constraints */

150:   /* nexn & nixn matricies for equaly and inequalty constriants */
151:   MatCreate(PETSC_COMM_WORLD,&user->Ae);
152:   MatCreate(PETSC_COMM_WORLD,&user->Ai);
153:   MatCreate(PETSC_COMM_WORLD,&user->H);

155:   MatSetSizes(user->Ae,neloc,nloc,user->ne,user->n);
156:   MatSetSizes(user->Ai,niloc,nloc,user->ni,user->n);
157:   MatSetSizes(user->H,nloc,nloc,user->n,user->n);

159:   MatSetUp(user->Ae);
160:   MatSetUp(user->Ai);
161:   MatSetUp(user->H);

163:   MatSetFromOptions(user->Ae);
164:   MatSetFromOptions(user->Ai);
165:   MatSetFromOptions(user->H);
166:   return(0);
167: }

169: PetscErrorCode DestroyProblem(AppCtx *user)
170: {

174:   MatDestroy(&user->Ae);
175:   MatDestroy(&user->Ai);
176:   MatDestroy(&user->H);

178:   VecDestroy(&user->x);
179:   VecDestroy(&user->ce);
180:   VecDestroy(&user->ci);
181:   VecDestroy(&user->xl);
182:   VecDestroy(&user->xu);
183:   VecDestroy(&user->Xseq);
184:   VecScatterDestroy(&user->scat);
185:   return(0);
186: }

188: /*
189:   f(X) = (x0 - 2)^2 + (x1 - 2)^2 - 2*(x0 + x1)
190:   G(X) = fx = [2.0*(x[0]-2.0) - 2.0; 2.0*(x[1]-2.0) - 2.0]
191: */
192: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ctx)
193: {
194:   PetscScalar       g;
195:   const PetscScalar *x;
196:   MPI_Comm          comm;
197:   PetscMPIInt       rank;
198:   PetscErrorCode    ierr;
199:   PetscReal         fin;
200:   AppCtx            *user=(AppCtx*)ctx;
201:   Vec               Xseq=user->Xseq;
202:   VecScatter        scat=user->scat;

205:   PetscObjectGetComm((PetscObject)tao,&comm);
206:   MPI_Comm_rank(comm,&rank);

208:   VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
209:   VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);

211:   fin = 0.0;
212:   if (!rank) {
213:     VecGetArrayRead(Xseq,&x);
214:     fin = (x[0]-2.0)*(x[0]-2.0) + (x[1]-2.0)*(x[1]-2.0) - 2.0*(x[0]+x[1]);
215:     g = 2.0*(x[0]-2.0) - 2.0;
216:     VecSetValue(G,0,g,INSERT_VALUES);
217:     g = 2.0*(x[1]-2.0) - 2.0;
218:     VecSetValue(G,1,g,INSERT_VALUES);
219:     VecRestoreArrayRead(Xseq,&x);
220:   }
221:   MPI_Allreduce(&fin,f,1,MPIU_REAL,MPIU_SUM,comm);
222:   VecAssemblyBegin(G);
223:   VecAssemblyEnd(G);
224:   return(0);
225: }

227: PetscErrorCode FormHessian(Tao tao, Vec x,Mat H, Mat Hpre, void *ctx)
228: {
229:   AppCtx            *user=(AppCtx*)ctx;
230:   Vec               DE,DI;
231:   const PetscScalar *de, *di;
232:   PetscInt          zero=0,one=1;
233:   PetscScalar       two=2.0;
234:   PetscScalar       val=0.0;
235:   Vec               Deseq,Diseq=user->Xseq;
236:   VecScatter        Descat,Discat=user->scat;
237:   PetscMPIInt       rank;
238:   MPI_Comm          comm;
239:   PetscErrorCode    ierr;

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

244:   PetscObjectGetComm((PetscObject)tao,&comm);
245:   MPI_Comm_rank(comm,&rank);

247:   VecScatterCreateToZero(DE,&Descat,&Deseq);

249:   VecScatterBegin(Descat,DE,Deseq,INSERT_VALUES,SCATTER_FORWARD);
250:   VecScatterBegin(Discat,DI,Diseq,INSERT_VALUES,SCATTER_FORWARD);
251:   VecScatterEnd(Descat,DE,Deseq,INSERT_VALUES,SCATTER_FORWARD);
252:   VecScatterEnd(Discat,DI,Diseq,INSERT_VALUES,SCATTER_FORWARD);

254:   if (!rank){
255:     VecGetArrayRead(Deseq,&de);  /* places equality constraint dual into array */
256:     VecGetArrayRead(Diseq,&di);  /* places inequality constraint dual into array */
257:     val = 2.0 * (1 + de[0] + di[0] - di[1]);
258:     VecRestoreArrayRead(Deseq,&de);
259:     VecRestoreArrayRead(Diseq,&di);
260:     MatSetValues(H,1,&zero,1,&zero,&val,INSERT_VALUES);
261:     MatSetValues(H,1,&one,1,&one,&two,INSERT_VALUES);
262:   }

264:   MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY);
265:   MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY);

267:   VecScatterDestroy(&Descat);
268:   VecDestroy(&Deseq);
269:   return(0);
270: }

272: PetscErrorCode FormInequalityConstraints(Tao tao,Vec X,Vec CI,void *ctx)
273: {
274:   const PetscScalar *x;
275:   PetscScalar       ci;
276:   PetscErrorCode    ierr;
277:   MPI_Comm          comm;
278:   PetscMPIInt       rank;
279:   AppCtx            *user=(AppCtx*)ctx;
280:   Vec               Xseq=user->Xseq;
281:   VecScatter        scat=user->scat;

284:   PetscObjectGetComm((PetscObject)tao,&comm);
285:   MPI_Comm_rank(comm,&rank);

287:   VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
288:   VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);

290:   if (!rank) {
291:     VecGetArrayRead(Xseq,&x);
292:     ci = x[0]*x[0] - x[1];
293:     VecSetValue(CI,0,ci,INSERT_VALUES);
294:     ci = -x[0]*x[0] + x[1] + 1.0;
295:     VecSetValue(CI,1,ci,INSERT_VALUES);
296:     VecRestoreArrayRead(Xseq,&x);
297:   }
298:   VecAssemblyBegin(CI);
299:   VecAssemblyEnd(CI);
300:   return(0);
301: }

303: PetscErrorCode FormEqualityConstraints(Tao tao,Vec X,Vec CE,void *ctx)
304: {
305:   const PetscScalar *x;
306:   PetscScalar       ce;
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:     ce = x[0]*x[0] + x[1] - 2.0;
324:     VecSetValue(CE,0,ce,INSERT_VALUES);
325:     VecRestoreArrayRead(Xseq,&x);
326:   }
327:   VecAssemblyBegin(CE);
328:   VecAssemblyEnd(CE);
329:   return(0);
330: }

332: PetscErrorCode FormInequalityJacobian(Tao tao, Vec X, Mat JI, Mat JIpre,  void *ctx)
333: {
334:   AppCtx            *user=(AppCtx*)ctx;
335:   PetscInt          cols[2],min,max,i;
336:   PetscScalar       vals[2];
337:   const PetscScalar *x;
338:   PetscErrorCode    ierr;
339:   Vec               Xseq=user->Xseq;
340:   VecScatter        scat=user->scat;
341:   MPI_Comm          comm;
342:   PetscMPIInt       rank;

345:   PetscObjectGetComm((PetscObject)tao,&comm);
346:   MPI_Comm_rank(comm,&rank);
347:   VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
348:   VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);

350:   VecGetArrayRead(Xseq,&x);
351:   MatGetOwnershipRange(JI,&min,&max);

353:   cols[0] = 0; cols[1] = 1;
354:   for (i=min;i<max;i++) {
355:     if (i==0){
356:       vals[0] = +2*x[0]; vals[1] = -1.0;
357:       MatSetValues(JI,1,&i,2,cols,vals,INSERT_VALUES);
358:     }
359:     if (i==1) {
360:       vals[0] = -2*x[0]; vals[1] = +1.0;
361:       MatSetValues(JI,1,&i,2,cols,vals,INSERT_VALUES);
362:     }
363:   }
364:   VecRestoreArrayRead(Xseq,&x);

366:   MatAssemblyBegin(JI,MAT_FINAL_ASSEMBLY);
367:   MatAssemblyEnd(JI,MAT_FINAL_ASSEMBLY);
368:   return(0);
369: }

371: PetscErrorCode FormEqualityJacobian(Tao tao,Vec X,Mat JE,Mat JEpre,void *ctx)
372: {
373:   PetscInt          rows[2];
374:   PetscScalar       vals[2];
375:   const PetscScalar *x;
376:   PetscMPIInt       rank;
377:   MPI_Comm          comm;
378:   PetscErrorCode    ierr;

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

384:   if (!rank) {
385:     VecGetArrayRead(X,&x);
386:     rows[0] = 0;       rows[1] = 1;
387:     vals[0] = 2*x[0];  vals[1] = 1.0;
388:     MatSetValues(JE,1,rows,2,rows,vals,INSERT_VALUES);
389:     VecRestoreArrayRead(X,&x);
390:   }
391:   MatAssemblyBegin(JE,MAT_FINAL_ASSEMBLY);
392:   MatAssemblyEnd(JE,MAT_FINAL_ASSEMBLY);
393:   return(0);
394: }


397: /*TEST

399:    build:
400:       requires: !complex !define(PETSC_USE_CXX)

402:    test:
403:       requires: superlu_dist
404:       args: -tao_converged_reason

406:    test:
407:       suffix: 2
408:       requires: superlu_dist
409:       nsize: 2
410:       args: -tao_converged_reason

412: TEST*/