Actual source code: toyf.F

petsc-3.5.2 2014-09-08
Report Typos and Errors
  1: ! Program usage: mpirun -np 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,1.0d-12,0,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('-pc_factor_mat_solver_package',         &
 90:      &     '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,1.0d-12,0.0d0,0.0d0,0.0d0,0.0d0,ierr)
103:       CHKERRQ(ierr)

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


110:       ! Analyze solution
111:       call TaoGetConvergedReason(tao,reason,ierr)
112:       CHKERRQ(ierr)

114: !      if (reason .lt. 0) then
115: !         call PetscPrintf(MPI_COMM_WORLD, "TAO failed to converge.\n"!,  &
116: !     &         ierr)
117: !      else
118: !         call PetscPrintf(MPI_COMM_WORLD,                            !   &
119: !     &         "Optimization terminated with status %D.\n", reason,ierr)
120: !      end if

122:       ! Finalize Memory
123:       call DestroyProblem(ierr)
124:       CHKERRQ(ierr)

126:       call TaoDestroy(tao,ierr)
127:       CHKERRQ(ierr)

129:       call PetscFinalize(ierr)

131:       stop
132:       end program toyf


135:       subroutine InitializeProblem(ierr)
136:       implicit none
137:  #include toyf.h
138:       PetscReal zero,minus1,two
139:       PetscErrorCode ierr
140:       n = 2
141:       zero =0.0d0
142:       minus1=-1.0d0
143:       two=2.0d0

145:       call VecCreateSeq(PETSC_COMM_SELF,n,x0,ierr)
146:       CHKERRQ(ierr)
147:       call VecDuplicate(x0,xl,ierr)
148:       CHKERRQ(ierr)
149:       call VecDuplicate(x0,xu,ierr)
150:       CHKERRQ(ierr)
151:       call VecSet(x0,zero,ierr)
152:       CHKERRQ(ierr)
153:       call VecSet(xl,minus1,ierr)
154:       CHKERRQ(ierr)
155:       call VecSet(xu,two,ierr)
156:       CHKERRQ(ierr)

158:       ne = 1
159:       call VecCreateSeq(PETSC_COMM_SELF,ne,ce,ierr)
160:       CHKERRQ(ierr)

162:       ni = 2
163:       call VecCreateSeq(PETSC_COMM_SELF,ni,ci,ierr)
164:       CHKERRQ(ierr)

166:       call MatCreateSeqAIJ(PETSC_COMM_SELF,ne,n,n,PETSC_NULL_INTEGER,Ae,&
167:      &     ierr)
168:       CHKERRQ(ierr)
169:       call MatCreateSeqAIJ(PETSC_COMM_SELF,ni,n,n,PETSC_NULL_INTEGER,Ai,&
170:      &     ierr)
171:       CHKERRQ(ierr)
172:       call MatSetFromOptions(Ae,ierr)
173:       CHKERRQ(ierr)
174:       call MatSetFromOptions(Ai,ierr)
175:       CHKERRQ(ierr)


178:       call MatCreateSeqAIJ(PETSC_COMM_SELF,n,n,1,PETSC_NULL_INTEGER,Hess&
179:      &     ,ierr)
180:       CHKERRQ(ierr)
181:       call MatSetFromOptions(Hess,ierr)
182:       CHKERRQ(ierr)
183:       0
184:       end subroutine InitializeProblem


187:       subroutine DestroyProblem(ierr)
188:       implicit none
189:  #include toyf.h

191:       PetscErrorCode ierr

193:       call MatDestroy(Ae,ierr)
194:       CHKERRQ(ierr)
195:       call MatDestroy(Ai,ierr)
196:       CHKERRQ(ierr)
197:       call MatDestroy(Hess,ierr)
198:       CHKERRQ(ierr)

200:       call VecDestroy(x0,ierr)
201:       CHKERRQ(ierr)
202:       call VecDestroy(ce,ierr)
203:       CHKERRQ(ierr)
204:       call VecDestroy(ci,ierr)
205:       CHKERRQ(ierr)
206:       call VecDestroy(xl,ierr)
207:       CHKERRQ(ierr)
208:       call VecDestroy(xu,ierr)
209:       CHKERRQ(ierr)
210:       0
211:       end subroutine DestroyProblem

213:       subroutine FormFunctionGradient(tao, X, f, G, dummy, ierr)
214:       implicit none
215:  #include toyf.h

217:       PetscErrorCode ierr
218:       PetscInt dummy
219:       Vec X,G
220:       Tao tao
221:       PetscScalar f
222:       PetscScalar x_v(0:1),g_v(0:1)
223:       PetscOffset x_i,g_i


226:       call VecGetArray(X,x_v,x_i,ierr)
227:       CHKERRQ(ierr)
228:       call VecGetArray(G,g_v,g_i,ierr)
229:       CHKERRQ(ierr)
230:       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)  &
231:      &       - 2.0*(x_v(x_i)+x_v(x_i+1))
232:       g_v(g_i) = 2.0*(x_v(x_i)-2.0) - 2.0
233:       g_v(g_i+1) = 2.0*(x_v(x_i+1)-2.0) - 2.0
234:       call VecRestoreArray(X,x_v,x_i,ierr)
235:       CHKERRQ(ierr)
236:       call VecRestoreArray(G,g_v,g_i,ierr)
237:       CHKERRQ(ierr)
238:       0
239:       end subroutine FormFunctionGradient


242:       subroutine FormHessian(tao,X,H,Hpre,dummy,ierr)
243:       implicit none
244:  #include toyf.h

246:       Tao        tao
247:       Vec              X
248:       Mat              H, Hpre
249:       PetscErrorCode   ierr
250:       PetscInt         dummy

252:       PetscScalar      de_v(0:1),di_v(0:1)
253:       PetscOffset      de_i,di_i
254:       PetscInt         zero(1)
255:       PetscInt         one(1)
256:       PetscScalar      two(1)
257:       PetscScalar      val(1)
258:       Vec DE,DI
259:       zero(1) = 0
260:       one(1) = 1
261:       two(1) = 2.0d0


264:       ! fix indices on matsetvalues
265:       call TaoGetDualVariables(tao,DE,DI,ierr)
266:       CHKERRQ(ierr)

268:       call VecGetArray(DE,de_v,de_i,ierr)
269:       CHKERRQ(ierr)
270:       call VecGetArray(DI,di_v,di_i,ierr)
271:       CHKERRQ(ierr)

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

275:       call VecRestoreArray(DE,de_v,de_i,ierr)
276:       CHKERRQ(ierr)
277:       call VecRestoreArray(DI,di_v,di_i,ierr)
278:       CHKERRQ(ierr)

280:       call MatSetValues(H,1,zero,1,zero,val,INSERT_VALUES,ierr)
281:       CHKERRQ(ierr)
282:       call MatSetValues(H,1,one,1,one,two,INSERT_VALUES,ierr)
283:       CHKERRQ(ierr)

285:       call MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY,ierr)
286:       CHKERRQ(ierr)
287:       call MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY,ierr)
288:       CHKERRQ(ierr)

290:       flag = SAME_NONZERO_PATTERN
291:       0
292:       end subroutine FormHessian

294:       subroutine FormInequalityConstraints(tao,X,C,dummy,ierr)
295:       implicit none
296:  #include toyf.h
297:       Tao      tao
298:       Vec            X,C
299:       PetscInt       dummy
300:       PetscErrorCode ierr
301:       PetscScalar    x_v(0:1),c_v(0:1)
302:       PetscOffset    x_i,c_i

304:       call VecGetArray(X,x_v,x_i,ierr)
305:       CHKERRQ(ierr)
306:       call VecGetArray(C,c_v,c_i,ierr)
307:       CHKERRQ(ierr)
308:       c_v(c_i) = x_v(x_i)*x_v(x_i) - x_v(x_i+1)
309:       c_v(c_i+1) = -x_v(x_i)*x_v(x_i) + x_v(x_i+1) + 1.0d0
310:       call VecRestoreArray(X,x_v,x_i,ierr)
311:       CHKERRQ(ierr)
312:       call VecRestoreArray(C,c_v,c_i,ierr)
313:       CHKERRQ(ierr)

315:       0
316:       end subroutine FormInequalityConstraints


319:       subroutine FormEqualityConstraints(tao,X,C,dummy,ierr)
320:       implicit none
321:  #include toyf.h
322:       Tao      tao
323:       Vec            X,C
324:       PetscInt       dummy
325:       PetscErrorCode ierr
326:       PetscScalar    x_v(0:1),c_v(0:1)
327:       PetscOffset    x_i,c_i
328:       call VecGetArray(X,x_v,x_i,ierr)
329:       CHKERRQ(ierr)
330:       call VecGetArray(C,c_v,c_i,ierr)
331:       CHKERRQ(ierr)
332:       c_v(c_i) = x_v(x_i)*x_v(x_i) + x_v(x_i+1) - 2.0d0
333:       call VecRestoreArray(X,x_v,x_i,ierr)
334:       CHKERRQ(ierr)
335:       call VecRestoreArray(C,c_v,c_i,ierr)
336:       CHKERRQ(ierr)
337:       0
338:       end subroutine FormEqualityConstraints


341:       subroutine FormInequalityJacobian(tao,X,JI,JIpre,dummy,ierr)
342:       implicit none
343:  #include toyf.h

345:       Tao       tao
346:       Vec             X
347:       Mat             JI,JIpre
348:       PetscInt        dummy
349:       PetscErrorCode  ierr

351:       PetscInt        rows(2)
352:       PetscInt        cols(2)
353:       PetscScalar     vals(4),x_v(0:1)
354:       PetscOffset     x_i

356:       call VecGetArray(X,x_v,x_i,ierr)
357:       CHKERRQ(ierr)
358:       rows(1)=0
359:       rows(2) = 1
360:       cols(1) = 0
361:       cols(2) = 1
362:       vals(1) = 2.0*x_v(x_i)
363:       vals(2) = -1.0d0
364:       vals(3) = -2.0*x_v(x_i)
365:       vals(4) = 1.0d0

367:       call VecRestoreArray(X,x_v,x_i,ierr)
368:       CHKERRQ(ierr)
369:       call MatSetValues(JI,2,rows,2,cols,vals,INSERT_VALUES,ierr)
370:       CHKERRQ(ierr)
371:       call MatAssemblyBegin(JI,MAT_FINAL_ASSEMBLY,ierr)
372:       CHKERRQ(ierr)
373:       call MatAssemblyEnd(JI,MAT_FINAL_ASSEMBLY,ierr)
374:       CHKERRQ(ierr)
375:       0
376:       end subroutine FormInequalityJacobian

378:       subroutine FormEqualityJacobian(tao,X,JE,JEpre,dummy,ierr)
379:       implicit none
380:  #include toyf.h

382:       Tao       tao
383:       Vec             X
384:       Mat             JE,JEpre
385:       PetscInt        dummy
386:       PetscErrorCode  ierr

388:       PetscInt        rows(2)
389:       PetscScalar     vals(4),x_v(0:1)
390:       PetscOffset     x_i

392:       call VecGetArray(X,x_v,x_i,ierr)
393:       CHKERRQ(ierr)
394:       rows(1)=0
395:       rows(2) = 1
396:       vals(1) = 2.0*x_v(x_i)
397:       vals(2) = 1.0d0

399:       call VecRestoreArray(X,x_v,x_i,ierr)
400:       CHKERRQ(ierr)
401:       call MatSetValues(JE,1,rows,2,rows,vals,INSERT_VALUES,ierr)
402:       CHKERRQ(ierr)
403:       call MatAssemblyBegin(JE,MAT_FINAL_ASSEMBLY,ierr)
404:       CHKERRQ(ierr)
405:       call MatAssemblyEnd(JE,MAT_FINAL_ASSEMBLY,ierr)
406:       CHKERRQ(ierr)
407:       0
408:       end subroutine FormEqualityJacobian