Actual source code: toy.c

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

  3: /* ----------------------------------------------------------------------
  4: min f=(x1-x2)^2 + (x2-2)^2 -2*x1-2*x2
  5: s.t.     x1^2 + x2 = 2
  6:       0 <= x1^2 - x2 <= 1
  7:       -1 <= x1,x2 <= 2
  8: ---------------------------------------------------------------------- */

 10:  #include <petsctao.h>

 12: static  char help[]="";

 14: /*
 15:    User-defined Section 1.5 Writing Application Codes with PETSc context - contains data needed by the
 16:    Section 1.5 Writing Application Codes with PETSc-provided call-back routines, FormFunction(),
 17:    FormGradient(), and FormHessian().
 18: */

 20: /*
 21:    x,d in R^n
 22:    f in R
 23:    bin in R^mi
 24:    beq in R^me
 25:    Aeq in R^(me x n)
 26:    Ain in R^(mi x n)
 27:    H in R^(n x n)
 28:    min f=(1/2)*x'*H*x + d'*x
 29:    s.t.  Aeq*x == beq
 30:          Ain*x >= bin
 31: */
 32: typedef struct {
 33:   PetscInt n; /* Length x */
 34:   PetscInt ne; /* number of equality constraints */
 35:   PetscInt ni; /* number of inequality constraints */
 36:   Vec      x,xl,xu;
 37:   Vec      ce,ci,bl,bu;
 38:   Mat      Ae,Ai,H;
 39: } AppCtx;

 41: /* -------- User-defined Routines --------- */

 43: PetscErrorCode InitializeProblem(AppCtx *);
 44: PetscErrorCode DestroyProblem(AppCtx *);
 45: PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal *,Vec,void *);
 46: PetscErrorCode FormHessian(Tao,Vec,Mat,Mat, void*);
 47: PetscErrorCode FormInequalityConstraints(Tao,Vec,Vec,void*);
 48: PetscErrorCode FormEqualityConstraints(Tao,Vec,Vec,void*);
 49: PetscErrorCode FormInequalityJacobian(Tao,Vec,Mat,Mat, void*);
 50: PetscErrorCode FormEqualityJacobian(Tao,Vec,Mat,Mat, void*);



 54: PetscErrorCode main(int argc,char **argv)
 55: {
 56:   PetscErrorCode     ierr;                /* used to check for functions returning nonzeros */
 57:   Tao                tao;
 58:   KSP                ksp;
 59:   PC                 pc;
 60:   AppCtx             user;                /* Section 1.5 Writing Application Codes with PETSc context */

 62:   PetscInitialize(&argc,&argv,(char *)0,help);if (ierr) return ierr;
 63:   PetscPrintf(PETSC_COMM_WORLD,"\n---- TOY Problem -----\n");
 64:   PetscPrintf(PETSC_COMM_WORLD,"Solution should be f(1,1)=-2\n");
 65:   InitializeProblem(&user);
 66:   TaoCreate(PETSC_COMM_WORLD,&tao);
 67:   TaoSetType(tao,TAOIPM);
 68:   TaoSetInitialVector(tao,user.x);
 69:   TaoSetVariableBounds(tao,user.xl,user.xu);
 70:   TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void*)&user);

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

 75:   TaoSetJacobianEqualityRoutine(tao,user.Ae,user.Ae,FormEqualityJacobian,(void*)&user);
 76:   TaoSetJacobianInequalityRoutine(tao,user.Ai,user.Ai,FormInequalityJacobian,(void*)&user);
 77:   TaoSetHessianRoutine(tao,user.H,user.H,FormHessian,(void*)&user);
 78:   /* TaoSetTolerances(tao,0,0,0); */

 80:   TaoGetKSP(tao,&ksp);
 81:   KSPGetPC(ksp,&pc);
 82:   PCSetType(pc,PCLU);
 83:   /*
 84:       This algorithm produces matrices with zeros along the diagonal therefore we need to use
 85:     SuperLU which does partial pivoting
 86:   */
 87:   PCFactorSetMatSolverType(pc,MATSOLVERSUPERLU);
 88:   KSPSetType(ksp,KSPPREONLY);
 89:   KSPSetFromOptions(ksp);

 91:   TaoSetFromOptions(tao);
 92:   /* TaoSetTolerances(tao,0,0,0); */
 93:   TaoSolve(tao);

 95:   DestroyProblem(&user);
 96:   TaoDestroy(&tao);
 97:   PetscFinalize();
 98:   return ierr;
 99: }

101: PetscErrorCode InitializeProblem(AppCtx *user)
102: {

106:   user->n = 2;
107:   VecCreateSeq(PETSC_COMM_SELF,user->n,&user->x);
108:   VecDuplicate(user->x,&user->xl);
109:   VecDuplicate(user->x,&user->xu);
110:   VecSet(user->x,0.0);
111:   VecSet(user->xl,-1.0);
112:   VecSet(user->xu,2.0);

114:   user->ne = 1;
115:   VecCreateSeq(PETSC_COMM_SELF,user->ne,&user->ce);

117:   user->ni = 2;
118:   VecCreateSeq(PETSC_COMM_SELF,user->ni,&user->ci);

120:   MatCreateSeqAIJ(PETSC_COMM_SELF,user->ne,user->n,user->n,NULL,&user->Ae);
121:   MatCreateSeqAIJ(PETSC_COMM_SELF,user->ni,user->n,user->n,NULL,&user->Ai);
122:   MatSetFromOptions(user->Ae);
123:   MatSetFromOptions(user->Ai);


126:   MatCreateSeqAIJ(PETSC_COMM_SELF,user->n,user->n,1,NULL,&user->H);
127:   MatSetFromOptions(user->H);

129:   return(0);
130: }

132: PetscErrorCode DestroyProblem(AppCtx *user)
133: {

137:   MatDestroy(&user->Ae);
138:   MatDestroy(&user->Ai);
139:   MatDestroy(&user->H);

141:   VecDestroy(&user->x);
142:   VecDestroy(&user->ce);
143:   VecDestroy(&user->ci);
144:   VecDestroy(&user->xl);
145:   VecDestroy(&user->xu);
146:   return(0);
147: }

149: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ctx)
150: {
151:   PetscScalar       *g;
152:   const PetscScalar *x;
153:   PetscErrorCode    ierr;

156:   VecGetArrayRead(X,&x);
157:   VecGetArray(G,&g);
158:   *f = (x[0]-2.0)*(x[0]-2.0) + (x[1]-2.0)*(x[1]-2.0) - 2.0*(x[0]+x[1]);
159:   g[0] = 2.0*(x[0]-2.0) - 2.0;
160:   g[1] = 2.0*(x[1]-2.0) - 2.0;
161:   VecRestoreArrayRead(X,&x);
162:   VecRestoreArray(G,&g);
163:   return(0);
164: }

166: PetscErrorCode FormHessian(Tao tao, Vec x, Mat H, Mat Hpre, void *ctx)
167: {
168:   Vec               DE,DI;
169:   const PetscScalar *de, *di;
170:   PetscInt          zero=0,one=1;
171:   PetscScalar       two=2.0;
172:   PetscScalar       val;
173:   PetscErrorCode    ierr;

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

178:   VecGetArrayRead(DE,&de);
179:   VecGetArrayRead(DI,&di);
180:   val=2.0 * (1 + de[0] + di[0] - di[1]);
181:   VecRestoreArrayRead(DE,&de);
182:   VecRestoreArrayRead(DI,&di);

184:   MatSetValues(H,1,&zero,1,&zero,&val,INSERT_VALUES);
185:   MatSetValues(H,1,&one,1,&one,&two,INSERT_VALUES);

187:   MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY);
188:   MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY);
189:   return(0);
190: }

192: PetscErrorCode FormInequalityConstraints(Tao tao, Vec X, Vec CI, void *ctx)
193: {
194:   const PetscScalar *x;
195:   PetscScalar       *c;
196:   PetscErrorCode    ierr;

199:   VecGetArrayRead(X,&x);
200:   VecGetArray(CI,&c);
201:   c[0] = x[0]*x[0] - x[1];
202:   c[1] = -x[0]*x[0] + x[1] + 1.0;
203:   VecRestoreArrayRead(X,&x);
204:   VecRestoreArray(CI,&c);
205:   return(0);
206: }

208: PetscErrorCode FormEqualityConstraints(Tao tao, Vec X, Vec CE,void *ctx)
209: {
210:   PetscScalar    *x,*c;

214:   VecGetArray(X,&x);
215:   VecGetArray(CE,&c);
216:   c[0] = x[0]*x[0] + x[1] - 2.0;
217:   VecRestoreArray(X,&x);
218:   VecRestoreArray(CE,&c);
219:   return(0);
220: }

222: PetscErrorCode FormInequalityJacobian(Tao tao, Vec X, Mat JI, Mat JIpre,  void *ctx)
223: {
224:   PetscInt          rows[2];
225:   PetscInt          cols[2];
226:   PetscScalar       vals[4];
227:   const PetscScalar *x;
228:   PetscErrorCode    ierr;

231:   VecGetArrayRead(X,&x);
232:   rows[0] = 0;       rows[1] = 1;
233:   cols[0] = 0;       cols[1] = 1;
234:   vals[0] = +2*x[0]; vals[1] = -1.0;
235:   vals[2] = -2*x[0]; vals[3] = +1.0;
236:   VecRestoreArrayRead(X,&x);
237:   MatSetValues(JI,2,rows,2,cols,vals,INSERT_VALUES);
238:   MatAssemblyBegin(JI,MAT_FINAL_ASSEMBLY);
239:   MatAssemblyEnd(JI,MAT_FINAL_ASSEMBLY);

241:   return(0);
242: }

244: PetscErrorCode FormEqualityJacobian(Tao tao, Vec X, Mat JE, Mat JEpre, void *ctx)
245: {
246:   PetscInt          rows[2];
247:   PetscScalar       vals[2];
248:   const PetscScalar *x;
249:   PetscErrorCode    ierr;

252:   VecGetArrayRead(X,&x);
253:   rows[0] = 0;       rows[1] = 1;
254:   vals[0] = 2*x[0];  vals[1] = 1.0;
255:   VecRestoreArrayRead(X,&x);
256:   MatSetValues(JE,1,rows,2,rows,vals,INSERT_VALUES);
257:   MatAssemblyBegin(JE,MAT_FINAL_ASSEMBLY);
258:   MatAssemblyEnd(JE,MAT_FINAL_ASSEMBLY);
259:   return(0);
260: }


263: /*TEST

265:    build:
266:       requires: !complex !define(PETSC_USE_CXX)

268:    test:
269:       requires: superlu
270:       args: -tao_smonitor -tao_view -tao_gatol 1.e-5

272: TEST*/