Actual source code: ex1.c

petsc-3.15.0 2021-04-05
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 mumps is requried to run either for 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;


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

 61: PetscErrorCode main(int argc,char **argv)
 62: {
 64:   Tao            tao;
 65:   KSP            ksp;
 66:   PC             pc;
 67:   AppCtx         user;  /* application context */
 68:   Vec            x;
 69:   PetscMPIInt    rank;
 70:   TaoType        type;
 71:   PetscBool      pdipm;

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

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

 87:   if (!user.noeqflag){
 88:     TaoSetEqualityConstraintsRoutine(tao,user.ce,FormEqualityConstraints,(void*)&user);
 89:   }
 90:   TaoSetInequalityConstraintsRoutine(tao,user.ci,FormInequalityConstraints,(void*)&user);

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

 99:   TaoGetKSP(tao,&ksp);
100:   KSPGetPC(ksp,&pc);
101:   PCSetType(pc,PCCHOLESKY);
102:   /*
103:       This algorithm produces matrices with zeros along the diagonal therefore we use
104:     MUMPS which provides solver for indefinite matrices
105:   */
106:   PCFactorSetMatSolverType(pc,MATSOLVERMUMPS);  /* requires mumps to solve pdipm */
107:   KSPSetType(ksp,KSPPREONLY);
108:   KSPSetFromOptions(ksp);

110:   TaoSetFromOptions(tao);
111:   TaoGetType(tao,&type);
112:   PetscObjectTypeCompare((PetscObject)tao,TAOPDIPM,&pdipm);
113:   if (pdipm) {
114:     TaoSetHessianRoutine(tao,user.H,user.H,FormHessian,(void*)&user);
115:   }

117:   TaoSolve(tao);
118:   TaoGetSolutionVector(tao,&x);
119:   VecView(x,PETSC_VIEWER_STDOUT_WORLD);

121:   /* Free objects */
122:   DestroyProblem(&user);
123:   TaoDestroy(&tao);
124:   PetscFinalize();
125:   return ierr;
126: }

128: PetscErrorCode InitializeProblem(AppCtx *user)
129: {
131:   PetscMPIInt    size;
132:   PetscMPIInt    rank;
133:   PetscInt       nloc,neloc,niloc;

136:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
137:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
138:   user->noeqflag = PETSC_FALSE;
139:   PetscOptionsGetBool(NULL,NULL,"-no_eq",&user->noeqflag,NULL);
140:   if (!user->noeqflag) {
141:     PetscPrintf(PETSC_COMM_WORLD,"Solution should be f(1,1)=-2\n");
142:   }

144:   /* create vector x and set initial values */
145:   user->n = 2; /* global length */
146:   nloc = (rank==0)?user->n:0;
147:   VecCreate(PETSC_COMM_WORLD,&user->x);
148:   VecSetSizes(user->x,nloc,user->n);
149:   VecSetFromOptions(user->x);
150:   VecSet(user->x,0);

152:   /* create and set lower and upper bound vectors */
153:   VecDuplicate(user->x,&user->xl);
154:   VecDuplicate(user->x,&user->xu);
155:   VecSet(user->xl,-1.0);
156:   VecSet(user->xu,2.0);

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

161:     user->ne = 1;
162:     neloc = (rank==0)?user->ne:0;

164:   if (!user->noeqflag){
165:     VecCreate(PETSC_COMM_WORLD,&user->ce); /* a 1x1 vec for equality constraints */
166:     VecSetSizes(user->ce,neloc,user->ne);
167:     VecSetFromOptions(user->ce);
168:     VecSetUp(user->ce);
169:   }
170:   user->ni = 2;
171:   niloc = (rank==0)?user->ni:0;
172:   VecCreate(PETSC_COMM_WORLD,&user->ci); /* a 2x1 vec for inequality constraints */
173:   VecSetSizes(user->ci,niloc,user->ni);
174:   VecSetFromOptions(user->ci);
175:   VecSetUp(user->ci);

177:   /* nexn & nixn matricies for equaly and inequalty constriants */
178:   if (!user->noeqflag){
179:     MatCreate(PETSC_COMM_WORLD,&user->Ae);
180:     MatSetSizes(user->Ae,neloc,nloc,user->ne,user->n);
181:     MatSetFromOptions(user->Ae);
182:     MatSetUp(user->Ae);
183:   }

185:   MatCreate(PETSC_COMM_WORLD,&user->Ai);
186:   MatCreate(PETSC_COMM_WORLD,&user->H);

188:   MatSetSizes(user->Ai,niloc,nloc,user->ni,user->n);
189:   MatSetSizes(user->H,nloc,nloc,user->n,user->n);

191:   MatSetFromOptions(user->Ai);
192:   MatSetFromOptions(user->H);

194:   MatSetUp(user->Ai);
195:   MatSetUp(user->H);
196:   return(0);
197: }

199: PetscErrorCode DestroyProblem(AppCtx *user)
200: {

204:   if (!user->noeqflag){
205:    MatDestroy(&user->Ae);
206:   }
207:   MatDestroy(&user->Ai);
208:   MatDestroy(&user->H);

210:   VecDestroy(&user->x);
211:   if (!user->noeqflag){
212:     VecDestroy(&user->ce);
213:   }
214:   VecDestroy(&user->ci);
215:   VecDestroy(&user->xl);
216:   VecDestroy(&user->xu);
217:   VecDestroy(&user->Xseq);
218:   VecScatterDestroy(&user->scat);
219:   return(0);
220: }

222: /*
223:   f(X) = (x0 - 2)^2 + (x1 - 2)^2 - 2*(x0 + x1)
224:   G(X) = fx = [2.0*(x[0]-2.0) - 2.0; 2.0*(x[1]-2.0) - 2.0]
225: */
226: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ctx)
227: {
228:   PetscScalar       g;
229:   const PetscScalar *x;
230:   MPI_Comm          comm;
231:   PetscMPIInt       rank;
232:   PetscErrorCode    ierr;
233:   PetscReal         fin;
234:   AppCtx            *user=(AppCtx*)ctx;
235:   Vec               Xseq=user->Xseq;
236:   VecScatter        scat=user->scat;

239:   PetscObjectGetComm((PetscObject)tao,&comm);
240:   MPI_Comm_rank(comm,&rank);

242:   VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
243:   VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);

245:   fin = 0.0;
246:   if (!rank) {
247:     VecGetArrayRead(Xseq,&x);
248:     fin = (x[0]-2.0)*(x[0]-2.0) + (x[1]-2.0)*(x[1]-2.0) - 2.0*(x[0]+x[1]);
249:     g = 2.0*(x[0]-2.0) - 2.0;
250:     VecSetValue(G,0,g,INSERT_VALUES);
251:     g = 2.0*(x[1]-2.0) - 2.0;
252:     VecSetValue(G,1,g,INSERT_VALUES);
253:     VecRestoreArrayRead(Xseq,&x);
254:   }
255:   MPI_Allreduce(&fin,f,1,MPIU_REAL,MPIU_SUM,comm);
256:   VecAssemblyBegin(G);
257:   VecAssemblyEnd(G);
258:   return(0);
259: }

261: PetscErrorCode FormHessian(Tao tao, Vec x,Mat H, Mat Hpre, void *ctx)
262: {
263:   AppCtx            *user=(AppCtx*)ctx;
264:   Vec               DE,DI;
265:   const PetscScalar *de, *di;
266:   PetscInt          zero=0,one=1;
267:   PetscScalar       two=2.0;
268:   PetscScalar       val=0.0;
269:   Vec               Deseq,Diseq=user->Xseq;
270:   VecScatter        Descat,Discat=user->scat;
271:   PetscMPIInt       rank;
272:   MPI_Comm          comm;
273:   PetscErrorCode    ierr;

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

278:   PetscObjectGetComm((PetscObject)tao,&comm);
279:   MPI_Comm_rank(comm,&rank);

281:   if (!user->noeqflag){
282:    VecScatterCreateToZero(DE,&Descat,&Deseq);
283:    VecScatterBegin(Descat,DE,Deseq,INSERT_VALUES,SCATTER_FORWARD);
284:    VecScatterEnd(Descat,DE,Deseq,INSERT_VALUES,SCATTER_FORWARD);
285:   }

287:   VecScatterBegin(Discat,DI,Diseq,INSERT_VALUES,SCATTER_FORWARD);
288:   VecScatterEnd(Discat,DI,Diseq,INSERT_VALUES,SCATTER_FORWARD);

290:   if (!rank){
291:     if (!user->noeqflag){
292:       VecGetArrayRead(Deseq,&de);  /* places equality constraint dual into array */
293:     }

295:     VecGetArrayRead(Diseq,&di);  /* places inequality constraint dual into array */
296:     if (!user->noeqflag){
297:       val = 2.0 * (1 + de[0] + di[0] - di[1]);
298:       VecRestoreArrayRead(Deseq,&de);
299:     }else{
300:       val = 2.0 * (1 + di[0] - di[1]);
301:     }
302:     VecRestoreArrayRead(Diseq,&di);
303:     MatSetValues(H,1,&zero,1,&zero,&val,INSERT_VALUES);
304:     MatSetValues(H,1,&one,1,&one,&two,INSERT_VALUES);
305:   }

307:   MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY);
308:   MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY);
309:   if (!user->noeqflag){
310:     VecScatterDestroy(&Descat);
311:     VecDestroy(&Deseq);
312:   }
313:   return(0);
314: }

316: PetscErrorCode FormInequalityConstraints(Tao tao,Vec X,Vec CI,void *ctx)
317: {
318:   const PetscScalar *x;
319:   PetscScalar       ci;
320:   PetscErrorCode    ierr;
321:   MPI_Comm          comm;
322:   PetscMPIInt       rank;
323:   AppCtx            *user=(AppCtx*)ctx;
324:   Vec               Xseq=user->Xseq;
325:   VecScatter        scat=user->scat;

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

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

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

347: PetscErrorCode FormEqualityConstraints(Tao tao,Vec X,Vec CE,void *ctx)
348: {
349:   const PetscScalar *x;
350:   PetscScalar       ce;
351:   PetscErrorCode    ierr;
352:   MPI_Comm          comm;
353:   PetscMPIInt       rank;
354:   AppCtx            *user=(AppCtx*)ctx;
355:   Vec               Xseq=user->Xseq;
356:   VecScatter        scat=user->scat;

359:   PetscObjectGetComm((PetscObject)tao,&comm);
360:   MPI_Comm_rank(comm,&rank);

362:   VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
363:   VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);

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

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

389:   PetscObjectGetComm((PetscObject)tao,&comm);
390:   MPI_Comm_rank(comm,&rank);
391:   VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
392:   VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);

394:   VecGetArrayRead(Xseq,&x);
395:   MatGetOwnershipRange(JI,&min,&max);

397:   cols[0] = 0; cols[1] = 1;
398:   for (i=min;i<max;i++) {
399:     if (i==0){
400:       vals[0] = +2*x[0]; vals[1] = -1.0;
401:       MatSetValues(JI,1,&i,2,cols,vals,INSERT_VALUES);
402:     }
403:     if (i==1) {
404:       vals[0] = -2*x[0]; vals[1] = +1.0;
405:       MatSetValues(JI,1,&i,2,cols,vals,INSERT_VALUES);
406:     }
407:   }
408:   VecRestoreArrayRead(Xseq,&x);

410:   MatAssemblyBegin(JI,MAT_FINAL_ASSEMBLY);
411:   MatAssemblyEnd(JI,MAT_FINAL_ASSEMBLY);
412:   return(0);
413: }

415: PetscErrorCode FormEqualityJacobian(Tao tao,Vec X,Mat JE,Mat JEpre,void *ctx)
416: {
417:   PetscInt          rows[2];
418:   PetscScalar       vals[2];
419:   const PetscScalar *x;
420:   PetscMPIInt       rank;
421:   MPI_Comm          comm;
422:   PetscErrorCode    ierr;

425:   PetscObjectGetComm((PetscObject)tao,&comm);
426:   MPI_Comm_rank(comm,&rank);

428:   if (!rank) {
429:     VecGetArrayRead(X,&x);
430:     rows[0] = 0;       rows[1] = 1;
431:     vals[0] = 2*x[0];  vals[1] = 1.0;
432:     MatSetValues(JE,1,rows,2,rows,vals,INSERT_VALUES);
433:     VecRestoreArrayRead(X,&x);
434:   }
435:   MatAssemblyBegin(JE,MAT_FINAL_ASSEMBLY);
436:   MatAssemblyEnd(JE,MAT_FINAL_ASSEMBLY);
437:   return(0);
438: }


441: /*TEST

443:    build:
444:       requires: !complex !define(PETSC_USE_CXX) mumps

446:    test:
447:       args: -tao_converged_reason -tao_pdipm_kkt_shift_pd

449:    test:
450:       suffix: 2
451:       nsize: 2
452:       args: -tao_converged_reason -tao_pdipm_kkt_shift_pd

454:    test:
455:       suffix: 3
456:       args: -tao_converged_reason -no_eq

458:    test:
459:       suffix: 4
460:       nsize: 2
461:       args: -tao_converged_reason -no_eq

463:    test:
464:       suffix: 5
465:       args: -tao_cmonitor -tao_type almm

467:    test:
468:       suffix: 6
469:       args: -tao_cmonitor -tao_type almm -tao_almm_type phr

471:    test:
472:       suffix: 7
473:       nsize: 2
474:       args: -tao_cmonitor -tao_type almm

476:    test:
477:       suffix: 8
478:       nsize: 2
479:       requires: cuda
480:       args: -tao_cmonitor -tao_type almm -vec_type cuda -mat_type aijcusparse

482:    test:
483:       suffix: 9
484:       nsize: 2
485:       args: -tao_cmonitor -tao_type almm -no_eq

487:    test:
488:       suffix: 10
489:       nsize: 2
490:       args: -tao_cmonitor -tao_type almm -tao_almm_type phr -no_eq

492: TEST*/