Actual source code: toyf.F90

petsc-3.11.3 2019-06-26
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:       module toymodule
 11: #include "petsc/finclude/petsctao.h"
 12:       use petsctao

 14:       Vec x0,xl,xu
 15:       Vec ce,ci,bl,bu
 16:       Mat Ae,Ai,Hess
 17:       PetscInt n,ne,ni
 18:       end module

 20:       program toyf
 21:       use toymodule
 22:       implicit none

 24:       PetscErrorCode       ierr
 25:       Tao                  tao
 26:       KSP                  ksp
 27:       PC                   pc
 28: !      PetscReal            zero

 30:       external FormFunctionGradient,FormHessian
 31:       external FormInequalityConstraints,FormEqualityConstraints
 32:       external FormInequalityJacobian,FormEqualityJacobian


 35:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
 36:       if (ierr .ne. 0) then
 37:          print*,'Unable to initialize PETSc'
 38:          stop
 39:       endif

 41:       call PetscPrintf(PETSC_COMM_SELF,'Solution should be f(1,1)=-2\n',ierr);CHKERRA(ierr)

 43:       call InitializeProblem(ierr);CHKERRA(ierr)

 45:       call TaoCreate(PETSC_COMM_SELF,tao,ierr);CHKERRA(ierr)
 46:       call TaoSetType(tao,TAOIPM,ierr);CHKERRA(ierr)
 47:       call TaoSetInitialVector(tao,x0,ierr);CHKERRA(ierr)
 48:       call TaoSetVariableBounds(tao,xl,xu,ierr);CHKERRA(ierr)
 49:       call TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,0,ierr);CHKERRA(ierr)
 50:       call TaoSetEqualityConstraintsRoutine(tao,ce,FormEqualityConstraints,0,ierr);CHKERRA(ierr)
 51:       call TaoSetInequalityConstraintsRoutine(tao,ci,FormInequalityConstraints,0,ierr);CHKERRA(ierr)
 52:       call TaoSetJacobianEqualityRoutine(tao,Ae,Ae,FormEqualityJacobian,0,ierr);CHKERRA(ierr)
 53:       call TaoSetJacobianInequalityRoutine(tao,Ai,Ai,FormInequalityJacobian,0,ierr);CHKERRA(ierr)
 54:       call TaoSetHessianRoutine(tao,Hess,Hess,FormHessian,0,ierr);CHKERRA(ierr)

 56:       call TaoSetFromOptions(tao,ierr);CHKERRA(ierr)
 57:       call TaoGetKSP(tao,ksp,ierr);CHKERRA(ierr)
 58:       call KSPGetPC(ksp,pc,ierr);CHKERRA(ierr)
 59:       call PCSetType(pc,PCLU,ierr);CHKERRA(ierr)
 60:       call PCFactorSetMatSolverType(pc,MATSOLVERSUPERLU,ierr);CHKERRA(ierr)
 61:       call KSPSetType(ksp,KSPPREONLY,ierr);CHKERRA(ierr)
 62:       call KSPSetFromOptions(ksp,ierr);CHKERRA(ierr)

 64: !      zero = 0;
 65: !      call TaoSetTolerances(tao,zero,zero,zero,ierr);CHKERRA(ierr)

 67:       ! Solve
 68:       call TaoSolve(tao,ierr);CHKERRA(ierr)

 70:       ! Finalize Memory
 71:       call DestroyProblem(ierr);CHKERRA(ierr)

 73:       call TaoDestroy(tao,ierr);CHKERRA(ierr)
 74:       call PetscFinalize(ierr)

 76:       stop
 77:       end program toyf


 80:       subroutine InitializeProblem(ierr)
 81:       use toymodule
 82:       implicit none
 83:       PetscReal zero,minus1,two,one
 84:       PetscInt done
 85:       PetscErrorCode ierr
 86:       n = 2
 87:       zero = 0
 88:       minus1 = -1
 89:       two= 2
 90:       one = 1
 91:       done = 1

 93:       call VecCreateSeq(PETSC_COMM_SELF,n,x0,ierr);CHKERRQ(ierr)
 94:       call VecDuplicate(x0,xl,ierr);CHKERRQ(ierr)
 95:       call VecDuplicate(x0,xu,ierr);CHKERRQ(ierr)
 96:       call VecSet(x0,zero,ierr);CHKERRQ(ierr)
 97:       call VecSet(xl,minus1,ierr);CHKERRQ(ierr)
 98:       call VecSet(xu,two,ierr);CHKERRQ(ierr)

100:       ne = 1
101:       call VecCreateSeq(PETSC_COMM_SELF,ne,ce,ierr);CHKERRQ(ierr)

103:       ni = 2
104:       call VecCreateSeq(PETSC_COMM_SELF,ni,ci,ierr);CHKERRQ(ierr)

106:       call MatCreateSeqAIJ(PETSC_COMM_SELF,ne,n,n,PETSC_NULL_INTEGER,Ae,ierr);CHKERRQ(ierr)
107:       call MatCreateSeqAIJ(PETSC_COMM_SELF,ni,n,n,PETSC_NULL_INTEGER,Ai,ierr);CHKERRQ(ierr)
108:       call MatSetFromOptions(Ae,ierr);CHKERRQ(ierr)
109:       call MatSetFromOptions(Ai,ierr);CHKERRQ(ierr)


112:       call MatCreateSeqAIJ(PETSC_COMM_SELF,n,n,done,PETSC_NULL_INTEGER,Hess,ierr);CHKERRQ(ierr)
113:       call MatSetFromOptions(Hess,ierr);CHKERRQ(ierr)
114:       0
115:       end subroutine InitializeProblem


118:       subroutine DestroyProblem(ierr)
119:       use toymodule
120:       implicit none

122:       PetscErrorCode ierr

124:       call MatDestroy(Ae,ierr);CHKERRQ(ierr)
125:       call MatDestroy(Ai,ierr);CHKERRQ(ierr)
126:       call MatDestroy(Hess,ierr);CHKERRQ(ierr)

128:       call VecDestroy(x0,ierr);CHKERRQ(ierr)
129:       call VecDestroy(ce,ierr);CHKERRQ(ierr)
130:       call VecDestroy(ci,ierr);CHKERRQ(ierr)
131:       call VecDestroy(xl,ierr);CHKERRQ(ierr)
132:       call VecDestroy(xu,ierr);CHKERRQ(ierr)
133:       0
134:       end subroutine DestroyProblem

136:       subroutine FormFunctionGradient(tao, X, f, G, dummy, ierr)
137:       use toymodule
138:       implicit none

140:       PetscErrorCode ierr
141:       PetscInt dummy
142:       Vec X,G
143:       Tao tao
144:       PetscScalar f
145:       PetscScalar x_v(0:1),g_v(0:1)
146:       PetscOffset x_i,g_i


149:       call VecGetArrayRead(X,x_v,x_i,ierr);CHKERRQ(ierr)
150:       call VecGetArray(G,g_v,g_i,ierr);CHKERRQ(ierr)
151:       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)-2.0*(x_v(x_i)+x_v(x_i+1))
152:       g_v(g_i) = 2.0*(x_v(x_i)-2.0) - 2.0
153:       g_v(g_i+1) = 2.0*(x_v(x_i+1)-2.0) - 2.0
154:       call VecRestoreArrayRead(X,x_v,x_i,ierr);CHKERRQ(ierr)
155:       call VecRestoreArray(G,g_v,g_i,ierr);CHKERRQ(ierr)
156:       0
157:       end subroutine FormFunctionGradient


160:       subroutine FormHessian(tao,X,H,Hpre,dummy,ierr)
161:       use toymodule
162:       implicit none

164:       Tao        tao
165:       Vec              X
166:       Mat              H, Hpre
167:       PetscErrorCode   ierr
168:       PetscInt         dummy

170:       PetscScalar      de_v(0:1),di_v(0:1)
171:       PetscOffset      de_i,di_i
172:       PetscInt         zero(1),ones
173:       PetscInt         one(1)
174:       PetscScalar      two(1)
175:       PetscScalar      val(1)
176:       Vec DE,DI
177:       zero(1) = 0
178:       one(1) = 1
179:       two(1) = 2.0d0
180:       ones = 1


183:       ! fix indices on matsetvalues
184:       call TaoGetDualVariables(tao,DE,DI,ierr);CHKERRQ(ierr)

186:       call VecGetArrayRead(DE,de_v,de_i,ierr);CHKERRQ(ierr)
187:       call VecGetArrayRead(DI,di_v,di_i,ierr);CHKERRQ(ierr)

189:       val(1)=2 * (1 + de_v(de_i) + di_v(di_i) - di_v(di_i+1))

191:       call VecRestoreArrayRead(DE,de_v,de_i,ierr);CHKERRQ(ierr)
192:       call VecRestoreArrayRead(DI,di_v,di_i,ierr);CHKERRQ(ierr)

194:       call MatSetValues(H,ones,zero,ones,zero,val,INSERT_VALUES,ierr);CHKERRQ(ierr)
195:       call MatSetValues(H,ones,one,ones,one,two,INSERT_VALUES,ierr);CHKERRQ(ierr)

197:       call MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr)
198:       call MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr)

200:       0
201:       end subroutine FormHessian

203:       subroutine FormInequalityConstraints(tao,X,C,dummy,ierr)
204:       use toymodule
205:       implicit none

207:       Tao      tao
208:       Vec            X,C
209:       PetscInt       dummy
210:       PetscErrorCode ierr
211:       PetscScalar    x_v(0:1),c_v(0:1)
212:       PetscOffset    x_i,c_i

214:       call VecGetArrayRead(X,x_v,x_i,ierr);CHKERRQ(ierr)
215:       call VecGetArray(C,c_v,c_i,ierr);CHKERRQ(ierr)
216:       c_v(c_i) = x_v(x_i)*x_v(x_i) - x_v(x_i+1)
217:       c_v(c_i+1) = -x_v(x_i)*x_v(x_i) + x_v(x_i+1) + 1
218:       call VecRestoreArrayRead(X,x_v,x_i,ierr);CHKERRQ(ierr)
219:       call VecRestoreArray(C,c_v,c_i,ierr);CHKERRQ(ierr)

221:       0
222:       end subroutine FormInequalityConstraints


225:       subroutine FormEqualityConstraints(tao,X,C,dummy,ierr)
226:       use toymodule
227:       implicit none

229:       Tao      tao
230:       Vec            X,C
231:       PetscInt       dummy
232:       PetscErrorCode ierr
233:       PetscScalar    x_v(0:1),c_v(0:1)
234:       PetscOffset    x_i,c_i
235:       call VecGetArrayRead(X,x_v,x_i,ierr);CHKERRQ(ierr)
236:       call VecGetArray(C,c_v,c_i,ierr);CHKERRQ(ierr)
237:       c_v(c_i) = x_v(x_i)*x_v(x_i) + x_v(x_i+1) - 2
238:       call VecRestoreArrayRead(X,x_v,x_i,ierr);CHKERRQ(ierr)
239:       call VecRestoreArray(C,c_v,c_i,ierr);CHKERRQ(ierr)
240:       0
241:       end subroutine FormEqualityConstraints


244:       subroutine FormInequalityJacobian(tao,X,JI,JIpre,dummy,ierr)
245:       use toymodule
246:       implicit none

248:       Tao       tao
249:       Vec             X
250:       Mat             JI,JIpre
251:       PetscInt        dummy
252:       PetscErrorCode  ierr

254:       PetscInt        rows(2),two
255:       PetscInt        cols(2)
256:       PetscScalar     vals(4),x_v(0:1)
257:       PetscOffset     x_i

259:       two = 2
260:       call VecGetArrayRead(X,x_v,x_i,ierr)
261:       CHKERRQ(ierr)
262:       rows(1)=0
263:       rows(2) = 1
264:       cols(1) = 0
265:       cols(2) = 1
266:       vals(1) = 2*x_v(x_i)
267:       vals(2) = -1
268:       vals(3) = -2*x_v(x_i)
269:       vals(4) = 1

271:       call VecRestoreArrayRead(X,x_v,x_i,ierr);CHKERRQ(ierr)
272:       call MatSetValues(JI,two,rows,two,cols,vals,INSERT_VALUES,ierr);CHKERRQ(ierr)
273:       call MatAssemblyBegin(JI,MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr)
274:       call MatAssemblyEnd(JI,MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr)
275:       0
276:       end subroutine FormInequalityJacobian

278:       subroutine FormEqualityJacobian(tao,X,JE,JEpre,dummy,ierr)
279:       use toymodule
280:       implicit none

282:       Tao       tao
283:       Vec             X
284:       Mat             JE,JEpre
285:       PetscInt        dummy
286:       PetscErrorCode  ierr

288:       PetscInt        rows(2),one,two
289:       PetscScalar     vals(4),x_v(0:1)
290:       PetscOffset     x_i
291:       one = 1
292:       two = 2
293: 
294:       call VecGetArrayRead(X,x_v,x_i,ierr);CHKERRQ(ierr)
295:       rows(1)=0
296:       rows(2) = 1
297:       vals(1) = 2*x_v(x_i)
298:       vals(2) = 1

300:       call VecRestoreArrayRead(X,x_v,x_i,ierr);CHKERRQ(ierr)
301:       call MatSetValues(JE,one,rows,two,rows,vals,INSERT_VALUES,ierr);CHKERRQ(ierr)
302:       call MatAssemblyBegin(JE,MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr)
303:       call MatAssemblyEnd(JE,MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr)
304:       0
305:       end subroutine FormEqualityJacobian

307: !/*TEST
308: !
309: !   build:
310: !      requires: !complex !single
311: !
312: !   test:
313: !      requires: superlu
314: !      args: -tao_monitor -tao_view -tao_gatol 1.e-5
315: !      filter: Error: grep -v IEEE_DENORMAL
316: !
317: !TEST*/