Actual source code: toyf.F

petsc-master 2017-03-25
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:  #include <petsc/finclude/petsctao.h>
 12:       use petsctao
 13:       implicit none
 14: #include "toyf.h"

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


 25:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
 26:       if (ierr .ne. 0) then
 27:          print*,'Unable to initialize PETSc'
 28:          stop
 29:       endif

 31:       call PetscPrintf(PETSC_COMM_SELF,                                 &
 32:      &           '\n---- TOY Problem -----\n',                          &
 33:      &           ierr)
 34:       CHKERRQ(ierr)

 36:       call PetscPrintf(PETSC_COMM_SELF,'Solution should be f(1,1)=-2\n',&
 37:      &     ierr)
 38:       CHKERRQ(ierr)

 40:       call InitializeProblem(ierr)
 41:       CHKERRQ(ierr)

 43:       call TaoCreate(PETSC_COMM_SELF,tao,ierr)
 44:       CHKERRQ(ierr)

 46:       call TaoSetType(tao,TAOIPM,ierr)
 47:       CHKERRQ(ierr)

 49:       call TaoSetInitialVector(tao,x0,ierr)
 50:       CHKERRQ(ierr)

 52:       call TaoSetVariableBounds(tao,xl,xu,ierr)
 53:       CHKERRQ(ierr)

 55:       call TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,   &
 56:      &     0,ierr)
 57:       CHKERRQ(ierr)

 59:       call TaoSetEqualityConstraintsRoutine(tao,ce,                      &
 60:      &     FormEqualityConstraints,0,ierr)
 61:       CHKERRQ(ierr)

 63:       call TaoSetInequalityConstraintsRoutine(tao,ci,                      &
 64:      &     FormInequalityConstraints,0,ierr)
 65:       CHKERRQ(ierr)

 67:       call TaoSetJacobianEqualityRoutine(tao,Ae,Ae,FormEqualityJacobian, &
 68:      &      0,ierr)
 69:       CHKERRQ(ierr)

 71:       call TaoSetJacobianInequalityRoutine(tao,Ai,Ai,                    &
 72:      &     FormInequalityJacobian,0,ierr)
 73:       CHKERRQ(ierr)

 75:       call TaoSetHessianRoutine(tao,Hess,Hess,FormHessian,               &
 76:      &     0,ierr)
 77:       CHKERRQ(ierr)

 79:       call TaoSetTolerances(tao,0,0,0,ierr)
 80:       CHKERRQ(ierr)

 82:       call TaoSetFromOptions(tao,ierr)
 83:       CHKERRQ(ierr)

 85:       call TaoGetKSP(tao,ksp,ierr)
 86:       CHKERRQ(ierr)

 88:       call KSPGetPC(ksp,pc,ierr)
 89:       CHKERRQ(ierr)

 91: !      call PCFactorSetMatSolverPackage(pc,MATSOLVERSUPERLU)
 92: !      CHKERRQ(ierr)

 94:       call PetscOptionsSetValue(PETSC_NULL_OPTIONS,                      &
 95:      &              '-pc_factor_mat_solver_package','superlu',ierr)
 96:       CHKERRQ(ierr)

 98:       call PCSetType(pc,PCLU,ierr)
 99:       CHKERRQ(ierr)

101:       call KSPSetType(ksp,KSPPREONLY,ierr)
102:       CHKERRQ(ierr)

104:       call KSPSetFromOptions(ksp,ierr)
105:       CHKERRQ(ierr)

107:       call TaoSetTolerances(tao,0.0d0,0.0d0,0.0d0,ierr)
108:       CHKERRQ(ierr)

110:       ! Solve
111:       call TaoSolve(tao,ierr)
112:       CHKERRQ(ierr)

114:       ! Finalize Memory
115:       call DestroyProblem(ierr)
116:       CHKERRQ(ierr)

118:       call TaoDestroy(tao,ierr)
119:       CHKERRQ(ierr)

121:       call PetscFinalize(ierr)

123:       stop
124:       end program toyf


127:       subroutine InitializeProblem(ierr)
128:       use petsctao
129:       implicit none
130: #include "toyf.h"
131:       PetscReal zero,minus1,two
132:       PetscErrorCode ierr
133:       n = 2
134:       zero =0.0d0
135:       minus1=-1.0d0
136:       two=2.0d0

138:       call VecCreateSeq(PETSC_COMM_SELF,n,x0,ierr)
139:       CHKERRQ(ierr)
140:       call VecDuplicate(x0,xl,ierr)
141:       CHKERRQ(ierr)
142:       call VecDuplicate(x0,xu,ierr)
143:       CHKERRQ(ierr)
144:       call VecSet(x0,zero,ierr)
145:       CHKERRQ(ierr)
146:       call VecSet(xl,minus1,ierr)
147:       CHKERRQ(ierr)
148:       call VecSet(xu,two,ierr)
149:       CHKERRQ(ierr)

151:       ne = 1
152:       call VecCreateSeq(PETSC_COMM_SELF,ne,ce,ierr)
153:       CHKERRQ(ierr)

155:       ni = 2
156:       call VecCreateSeq(PETSC_COMM_SELF,ni,ci,ierr)
157:       CHKERRQ(ierr)

159:       call MatCreateSeqAIJ(PETSC_COMM_SELF,ne,n,n,PETSC_NULL_INTEGER,Ae,&
160:      &     ierr)
161:       CHKERRQ(ierr)
162:       call MatCreateSeqAIJ(PETSC_COMM_SELF,ni,n,n,PETSC_NULL_INTEGER,Ai,&
163:      &     ierr)
164:       CHKERRQ(ierr)
165:       call MatSetFromOptions(Ae,ierr)
166:       CHKERRQ(ierr)
167:       call MatSetFromOptions(Ai,ierr)
168:       CHKERRQ(ierr)


171:       call MatCreateSeqAIJ(PETSC_COMM_SELF,n,n,1,PETSC_NULL_INTEGER,Hess&
172:      &     ,ierr)
173:       CHKERRQ(ierr)
174:       call MatSetFromOptions(Hess,ierr)
175:       CHKERRQ(ierr)
176:       0
177:       end subroutine InitializeProblem


180:       subroutine DestroyProblem(ierr)
181:       use petsctao
182:       implicit none
183: #include "toyf.h"

185:       PetscErrorCode ierr

187:       call MatDestroy(Ae,ierr)
188:       CHKERRQ(ierr)
189:       call MatDestroy(Ai,ierr)
190:       CHKERRQ(ierr)
191:       call MatDestroy(Hess,ierr)
192:       CHKERRQ(ierr)

194:       call VecDestroy(x0,ierr)
195:       CHKERRQ(ierr)
196:       call VecDestroy(ce,ierr)
197:       CHKERRQ(ierr)
198:       call VecDestroy(ci,ierr)
199:       CHKERRQ(ierr)
200:       call VecDestroy(xl,ierr)
201:       CHKERRQ(ierr)
202:       call VecDestroy(xu,ierr)
203:       CHKERRQ(ierr)
204:       0
205:       end subroutine DestroyProblem

207:       subroutine FormFunctionGradient(tao, X, f, G, dummy, ierr)
208:       use petsctao
209:       implicit none
210: #include "toyf.h"

212:       PetscErrorCode ierr
213:       PetscInt dummy
214:       Vec X,G
215:       Tao tao
216:       PetscScalar f
217:       PetscScalar x_v(0:1),g_v(0:1)
218:       PetscOffset x_i,g_i


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


237:       subroutine FormHessian(tao,X,H,Hpre,dummy,ierr)
238:       use petsctao
239:       implicit none
240: #include "toyf.h"

242:       Tao        tao
243:       Vec              X
244:       Mat              H, Hpre
245:       PetscErrorCode   ierr
246:       PetscInt         dummy

248:       PetscScalar      de_v(0:1),di_v(0:1)
249:       PetscOffset      de_i,di_i
250:       PetscInt         zero(1)
251:       PetscInt         one(1)
252:       PetscScalar      two(1)
253:       PetscScalar      val(1)
254:       Vec DE,DI
255:       zero(1) = 0
256:       one(1) = 1
257:       two(1) = 2.0d0


260:       ! fix indices on matsetvalues
261:       call TaoGetDualVariables(tao,DE,DI,ierr)
262:       CHKERRQ(ierr)

264:       call VecGetArray(DE,de_v,de_i,ierr)
265:       CHKERRQ(ierr)
266:       call VecGetArray(DI,di_v,di_i,ierr)
267:       CHKERRQ(ierr)

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

271:       call VecRestoreArray(DE,de_v,de_i,ierr)
272:       CHKERRQ(ierr)
273:       call VecRestoreArray(DI,di_v,di_i,ierr)
274:       CHKERRQ(ierr)

276:       call MatSetValues(H,1,zero,1,zero,val,INSERT_VALUES,ierr)
277:       CHKERRQ(ierr)
278:       call MatSetValues(H,1,one,1,one,two,INSERT_VALUES,ierr)
279:       CHKERRQ(ierr)

281:       call MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY,ierr)
282:       CHKERRQ(ierr)
283:       call MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY,ierr)
284:       CHKERRQ(ierr)

286:       0
287:       end subroutine FormHessian

289:       subroutine FormInequalityConstraints(tao,X,C,dummy,ierr)
290:       use petsctao
291:       implicit none
292: #include "toyf.h"
293:       Tao      tao
294:       Vec            X,C
295:       PetscInt       dummy
296:       PetscErrorCode ierr
297:       PetscScalar    x_v(0:1),c_v(0:1)
298:       PetscOffset    x_i,c_i

300:       call VecGetArray(X,x_v,x_i,ierr)
301:       CHKERRQ(ierr)
302:       call VecGetArray(C,c_v,c_i,ierr)
303:       CHKERRQ(ierr)
304:       c_v(c_i) = x_v(x_i)*x_v(x_i) - x_v(x_i+1)
305:       c_v(c_i+1) = -x_v(x_i)*x_v(x_i) + x_v(x_i+1) + 1.0d0
306:       call VecRestoreArray(X,x_v,x_i,ierr)
307:       CHKERRQ(ierr)
308:       call VecRestoreArray(C,c_v,c_i,ierr)
309:       CHKERRQ(ierr)

311:       0
312:       end subroutine FormInequalityConstraints


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


338:       subroutine FormInequalityJacobian(tao,X,JI,JIpre,dummy,ierr)
339:       use petsctao
340:       implicit none
341: #include "toyf.h"

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

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

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

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

376:       subroutine FormEqualityJacobian(tao,X,JE,JEpre,dummy,ierr)
377:       use petsctao
378:       implicit none
379: #include "toyf.h"

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

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

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

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