Actual source code: toyf.F
petsc-3.8.4 2018-03-24
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: CHKERRA(ierr)
36: call PetscPrintf(PETSC_COMM_SELF,'Solution should be f(1,1)=-2\n',&
37: & ierr)
38: CHKERRA(ierr)
40: call InitializeProblem(ierr)
41: CHKERRA(ierr)
43: call TaoCreate(PETSC_COMM_SELF,tao,ierr)
44: CHKERRA(ierr)
46: call TaoSetType(tao,TAOIPM,ierr)
47: CHKERRA(ierr)
49: call TaoSetInitialVector(tao,x0,ierr)
50: CHKERRA(ierr)
52: call TaoSetVariableBounds(tao,xl,xu,ierr)
53: CHKERRA(ierr)
55: call TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient, &
56: & 0,ierr)
57: CHKERRA(ierr)
59: call TaoSetEqualityConstraintsRoutine(tao,ce, &
60: & FormEqualityConstraints,0,ierr)
61: CHKERRA(ierr)
63: call TaoSetInequalityConstraintsRoutine(tao,ci, &
64: & FormInequalityConstraints,0,ierr)
65: CHKERRA(ierr)
67: call TaoSetJacobianEqualityRoutine(tao,Ae,Ae,FormEqualityJacobian, &
68: & 0,ierr)
69: CHKERRA(ierr)
71: call TaoSetJacobianInequalityRoutine(tao,Ai,Ai, &
72: & FormInequalityJacobian,0,ierr)
73: CHKERRA(ierr)
75: call TaoSetHessianRoutine(tao,Hess,Hess,FormHessian, &
76: & 0,ierr)
77: CHKERRA(ierr)
79: call TaoSetTolerances(tao,0,0,0,ierr)
80: CHKERRA(ierr)
82: call TaoSetFromOptions(tao,ierr)
83: CHKERRA(ierr)
85: call TaoGetKSP(tao,ksp,ierr)
86: CHKERRA(ierr)
88: call KSPGetPC(ksp,pc,ierr)
89: CHKERRA(ierr)
91: ! call PCFactorSetMatSolverPackage(pc,MATSOLVERSUPERLU)
92: ! CHKERRA(ierr)
94: call PetscOptionsSetValue(PETSC_NULL_OPTIONS, &
95: & '-pc_factor_mat_solver_package','superlu',ierr)
96: CHKERRA(ierr)
98: call PCSetType(pc,PCLU,ierr)
99: CHKERRA(ierr)
101: call KSPSetType(ksp,KSPPREONLY,ierr)
102: CHKERRA(ierr)
104: call KSPSetFromOptions(ksp,ierr)
105: CHKERRA(ierr)
107: call TaoSetTolerances(tao,0.0d0,0.0d0,0.0d0,ierr)
108: CHKERRA(ierr)
110: ! Solve
111: call TaoSolve(tao,ierr)
112: CHKERRA(ierr)
114: ! Finalize Memory
115: call DestroyProblem(ierr)
116: CHKERRA(ierr)
118: call TaoDestroy(tao,ierr)
119: CHKERRA(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