Actual source code: toyf.F
petsc-3.7.5 2017-01-01
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