Actual source code: toyf.F

petsc-3.7.3 2016-08-01
Report Typos and Errors
  1: ! Program usage: mpiexec -n 1 toyf[-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:       program toyf
 11:       implicit none
 12:  #include toyf.h

 14:       PetscErrorCode       ierr
 15:       Tao                  tao
 16:       TaoConvergedReason   reason
 17:       KSP                  ksp
 18:       PC                   pc
 19:       external FormFunctionGradient,FormHessian
 20:       external FormInequalityConstraints,FormEqualityConstraints
 21:       external FormInequalityJacobian,FormEqualityJacobian


 24:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)

 26:       call PetscPrintf(PETSC_COMM_SELF,                                 &
 27:      &           '\n---- TOY Problem -----\n',                          &
 28:      &           ierr)
 29:       CHKERRQ(ierr)

 31:       call PetscPrintf(PETSC_COMM_SELF,'Solution should be f(1,1)=-2\n',&
 32:      &     ierr)
 33:       CHKERRQ(ierr)

 35:       call InitializeProblem(ierr)
 36:       CHKERRQ(ierr)

 38:       call TaoCreate(PETSC_COMM_SELF,tao,ierr)
 39:       CHKERRQ(ierr)

 41:       call TaoSetType(tao,TAOIPM,ierr)
 42:       CHKERRQ(ierr)

 44:       call TaoSetInitialVector(tao,x0,ierr)
 45:       CHKERRQ(ierr)

 47:       call TaoSetVariableBounds(tao,xl,xu,ierr)
 48:       CHKERRQ(ierr)

 50:       call TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,   &
 51:      &     PETSC_NULL_OBJECT,ierr)
 52:       CHKERRQ(ierr)

 54:       call TaoSetEqualityConstraintsRoutine(tao,ce,                      &
 55:      &     FormEqualityConstraints,PETSC_NULL_OBJECT,ierr)
 56:       CHKERRQ(ierr)

 58:       call TaoSetInequalityConstraintsRoutine(tao,ci,                      &
 59:      &     FormInequalityConstraints,PETSC_NULL_OBJECT,ierr)
 60:       CHKERRQ(ierr)

 62:       call TaoSetJacobianEqualityRoutine(tao,Ae,Ae,FormEqualityJacobian, &
 63:      &      PETSC_NULL_OBJECT,ierr)
 64:       CHKERRQ(ierr)

 66:       call TaoSetJacobianInequalityRoutine(tao,Ai,Ai,                    &
 67:      &     FormInequalityJacobian,PETSC_NULL_OBJECT,ierr)
 68:       CHKERRQ(ierr)

 70:       call TaoSetHessianRoutine(tao,Hess,Hess,FormHessian,               &
 71:      &     PETSC_NULL_OBJECT,ierr)
 72:       CHKERRQ(ierr)

 74:       call TaoSetTolerances(tao,0,0,0,ierr)
 75:       CHKERRQ(ierr)

 77:       call TaoSetFromOptions(tao,ierr)
 78:       CHKERRQ(ierr)

 80:       call TaoGetKSP(tao,ksp,ierr)
 81:       CHKERRQ(ierr)

 83:       call KSPGetPC(ksp,pc,ierr)
 84:       CHKERRQ(ierr)

 86: !      call PCFactorSetMatSolverPackage(pc,MATSOLVERSUPERLU)
 87: !      CHKERRQ(ierr)

 89:       call PetscOptionsSetValue(PETSC_NULL_OBJECT,                      &
 90:      &              '-pc_factor_mat_solver_package','superlu',ierr)
 91:       CHKERRQ(ierr)

 93:       call PCSetType(pc,PCLU,ierr)
 94:       CHKERRQ(ierr)

 96:       call KSPSetType(ksp,KSPPREONLY,ierr)
 97:       CHKERRQ(ierr)

 99:       call KSPSetFromOptions(ksp,ierr)
100:       CHKERRQ(ierr)

102:       call TaoSetTolerances(tao,0.0d0,0.0d0,0.0d0,ierr)
103:       CHKERRQ(ierr)

105:       ! Solve
106:       call TaoSolve(tao,ierr)
107:       CHKERRQ(ierr)

109:       ! Finalize Memory
110:       call DestroyProblem(ierr)
111:       CHKERRQ(ierr)

113:       call TaoDestroy(tao,ierr)
114:       CHKERRQ(ierr)

116:       call PetscFinalize(ierr)

118:       stop
119:       end program toyf


122:       subroutine InitializeProblem(ierr)
123:       implicit none
124:  #include toyf.h
125:       PetscReal zero,minus1,two
126:       PetscErrorCode ierr
127:       n = 2
128:       zero =0.0d0
129:       minus1=-1.0d0
130:       two=2.0d0

132:       call VecCreateSeq(PETSC_COMM_SELF,n,x0,ierr)
133:       CHKERRQ(ierr)
134:       call VecDuplicate(x0,xl,ierr)
135:       CHKERRQ(ierr)
136:       call VecDuplicate(x0,xu,ierr)
137:       CHKERRQ(ierr)
138:       call VecSet(x0,zero,ierr)
139:       CHKERRQ(ierr)
140:       call VecSet(xl,minus1,ierr)
141:       CHKERRQ(ierr)
142:       call VecSet(xu,two,ierr)
143:       CHKERRQ(ierr)

145:       ne = 1
146:       call VecCreateSeq(PETSC_COMM_SELF,ne,ce,ierr)
147:       CHKERRQ(ierr)

149:       ni = 2
150:       call VecCreateSeq(PETSC_COMM_SELF,ni,ci,ierr)
151:       CHKERRQ(ierr)

153:       call MatCreateSeqAIJ(PETSC_COMM_SELF,ne,n,n,PETSC_NULL_INTEGER,Ae,&
154:      &     ierr)
155:       CHKERRQ(ierr)
156:       call MatCreateSeqAIJ(PETSC_COMM_SELF,ni,n,n,PETSC_NULL_INTEGER,Ai,&
157:      &     ierr)
158:       CHKERRQ(ierr)
159:       call MatSetFromOptions(Ae,ierr)
160:       CHKERRQ(ierr)
161:       call MatSetFromOptions(Ai,ierr)
162:       CHKERRQ(ierr)


165:       call MatCreateSeqAIJ(PETSC_COMM_SELF,n,n,1,PETSC_NULL_INTEGER,Hess&
166:      &     ,ierr)
167:       CHKERRQ(ierr)
168:       call MatSetFromOptions(Hess,ierr)
169:       CHKERRQ(ierr)
170:       0
171:       end subroutine InitializeProblem


174:       subroutine DestroyProblem(ierr)
175:       implicit none
176:  #include toyf.h

178:       PetscErrorCode ierr

180:       call MatDestroy(Ae,ierr)
181:       CHKERRQ(ierr)
182:       call MatDestroy(Ai,ierr)
183:       CHKERRQ(ierr)
184:       call MatDestroy(Hess,ierr)
185:       CHKERRQ(ierr)

187:       call VecDestroy(x0,ierr)
188:       CHKERRQ(ierr)
189:       call VecDestroy(ce,ierr)
190:       CHKERRQ(ierr)
191:       call VecDestroy(ci,ierr)
192:       CHKERRQ(ierr)
193:       call VecDestroy(xl,ierr)
194:       CHKERRQ(ierr)
195:       call VecDestroy(xu,ierr)
196:       CHKERRQ(ierr)
197:       0
198:       end subroutine DestroyProblem

200:       subroutine FormFunctionGradient(tao, X, f, G, dummy, ierr)
201:       implicit none
202:  #include toyf.h

204:       PetscErrorCode ierr
205:       PetscInt dummy
206:       Vec X,G
207:       Tao tao
208:       PetscScalar f
209:       PetscScalar x_v(0:1),g_v(0:1)
210:       PetscOffset x_i,g_i


213:       call VecGetArray(X,x_v,x_i,ierr)
214:       CHKERRQ(ierr)
215:       call VecGetArray(G,g_v,g_i,ierr)
216:       CHKERRQ(ierr)
217:       f=(x_v(x_i)-2.0)*(x_v(x_i)-2.0)+(x_v(x_i+1)-2.0)*(x_v(x_i+1)-2.0)  &
218:      &       - 2.0*(x_v(x_i)+x_v(x_i+1))
219:       g_v(g_i) = 2.0*(x_v(x_i)-2.0) - 2.0
220:       g_v(g_i+1) = 2.0*(x_v(x_i+1)-2.0) - 2.0
221:       call VecRestoreArray(X,x_v,x_i,ierr)
222:       CHKERRQ(ierr)
223:       call VecRestoreArray(G,g_v,g_i,ierr)
224:       CHKERRQ(ierr)
225:       0
226:       end subroutine FormFunctionGradient


229:       subroutine FormHessian(tao,X,H,Hpre,dummy,ierr)
230:       implicit none
231:  #include toyf.h

233:       Tao        tao
234:       Vec              X
235:       Mat              H, Hpre
236:       PetscErrorCode   ierr
237:       PetscInt         dummy

239:       PetscScalar      de_v(0:1),di_v(0:1)
240:       PetscOffset      de_i,di_i
241:       PetscInt         zero(1)
242:       PetscInt         one(1)
243:       PetscScalar      two(1)
244:       PetscScalar      val(1)
245:       Vec DE,DI
246:       zero(1) = 0
247:       one(1) = 1
248:       two(1) = 2.0d0


251:       ! fix indices on matsetvalues
252:       call TaoGetDualVariables(tao,DE,DI,ierr)
253:       CHKERRQ(ierr)

255:       call VecGetArray(DE,de_v,de_i,ierr)
256:       CHKERRQ(ierr)
257:       call VecGetArray(DI,di_v,di_i,ierr)
258:       CHKERRQ(ierr)

260:       val(1)=2.0d0 * (1.0d0 + de_v(de_i) + di_v(di_i) - di_v(di_i+1))

262:       call VecRestoreArray(DE,de_v,de_i,ierr)
263:       CHKERRQ(ierr)
264:       call VecRestoreArray(DI,di_v,di_i,ierr)
265:       CHKERRQ(ierr)

267:       call MatSetValues(H,1,zero,1,zero,val,INSERT_VALUES,ierr)
268:       CHKERRQ(ierr)
269:       call MatSetValues(H,1,one,1,one,two,INSERT_VALUES,ierr)
270:       CHKERRQ(ierr)

272:       call MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY,ierr)
273:       CHKERRQ(ierr)
274:       call MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY,ierr)
275:       CHKERRQ(ierr)

277:       flag = SAME_NONZERO_PATTERN
278:       0
279:       end subroutine FormHessian

281:       subroutine FormInequalityConstraints(tao,X,C,dummy,ierr)
282:       implicit none
283:  #include toyf.h
284:       Tao      tao
285:       Vec            X,C
286:       PetscInt       dummy
287:       PetscErrorCode ierr
288:       PetscScalar    x_v(0:1),c_v(0:1)
289:       PetscOffset    x_i,c_i

291:       call VecGetArray(X,x_v,x_i,ierr)
292:       CHKERRQ(ierr)
293:       call VecGetArray(C,c_v,c_i,ierr)
294:       CHKERRQ(ierr)
295:       c_v(c_i) = x_v(x_i)*x_v(x_i) - x_v(x_i+1)
296:       c_v(c_i+1) = -x_v(x_i)*x_v(x_i) + x_v(x_i+1) + 1.0d0
297:       call VecRestoreArray(X,x_v,x_i,ierr)
298:       CHKERRQ(ierr)
299:       call VecRestoreArray(C,c_v,c_i,ierr)
300:       CHKERRQ(ierr)

302:       0
303:       end subroutine FormInequalityConstraints


306:       subroutine FormEqualityConstraints(tao,X,C,dummy,ierr)
307:       implicit none
308:  #include toyf.h
309:       Tao      tao
310:       Vec            X,C
311:       PetscInt       dummy
312:       PetscErrorCode ierr
313:       PetscScalar    x_v(0:1),c_v(0:1)
314:       PetscOffset    x_i,c_i
315:       call VecGetArray(X,x_v,x_i,ierr)
316:       CHKERRQ(ierr)
317:       call VecGetArray(C,c_v,c_i,ierr)
318:       CHKERRQ(ierr)
319:       c_v(c_i) = x_v(x_i)*x_v(x_i) + x_v(x_i+1) - 2.0d0
320:       call VecRestoreArray(X,x_v,x_i,ierr)
321:       CHKERRQ(ierr)
322:       call VecRestoreArray(C,c_v,c_i,ierr)
323:       CHKERRQ(ierr)
324:       0
325:       end subroutine FormEqualityConstraints


328:       subroutine FormInequalityJacobian(tao,X,JI,JIpre,dummy,ierr)
329:       implicit none
330:  #include toyf.h

332:       Tao       tao
333:       Vec             X
334:       Mat             JI,JIpre
335:       PetscInt        dummy
336:       PetscErrorCode  ierr

338:       PetscInt        rows(2)
339:       PetscInt        cols(2)
340:       PetscScalar     vals(4),x_v(0:1)
341:       PetscOffset     x_i

343:       call VecGetArray(X,x_v,x_i,ierr)
344:       CHKERRQ(ierr)
345:       rows(1)=0
346:       rows(2) = 1
347:       cols(1) = 0
348:       cols(2) = 1
349:       vals(1) = 2.0*x_v(x_i)
350:       vals(2) = -1.0d0
351:       vals(3) = -2.0*x_v(x_i)
352:       vals(4) = 1.0d0

354:       call VecRestoreArray(X,x_v,x_i,ierr)
355:       CHKERRQ(ierr)
356:       call MatSetValues(JI,2,rows,2,cols,vals,INSERT_VALUES,ierr)
357:       CHKERRQ(ierr)
358:       call MatAssemblyBegin(JI,MAT_FINAL_ASSEMBLY,ierr)
359:       CHKERRQ(ierr)
360:       call MatAssemblyEnd(JI,MAT_FINAL_ASSEMBLY,ierr)
361:       CHKERRQ(ierr)
362:       0
363:       end subroutine FormInequalityJacobian

365:       subroutine FormEqualityJacobian(tao,X,JE,JEpre,dummy,ierr)
366:       implicit none
367:  #include toyf.h

369:       Tao       tao
370:       Vec             X
371:       Mat             JE,JEpre
372:       PetscInt        dummy
373:       PetscErrorCode  ierr

375:       PetscInt        rows(2)
376:       PetscScalar     vals(4),x_v(0:1)
377:       PetscOffset     x_i

379:       call VecGetArray(X,x_v,x_i,ierr)
380:       CHKERRQ(ierr)
381:       rows(1)=0
382:       rows(2) = 1
383:       vals(1) = 2.0*x_v(x_i)
384:       vals(2) = 1.0d0

386:       call VecRestoreArray(X,x_v,x_i,ierr)
387:       CHKERRQ(ierr)
388:       call MatSetValues(JE,1,rows,2,rows,vals,INSERT_VALUES,ierr)
389:       CHKERRQ(ierr)
390:       call MatAssemblyBegin(JE,MAT_FINAL_ASSEMBLY,ierr)
391:       CHKERRQ(ierr)
392:       call MatAssemblyEnd(JE,MAT_FINAL_ASSEMBLY,ierr)
393:       CHKERRQ(ierr)
394:       0
395:       end subroutine FormEqualityJacobian