Actual source code: toyf.F
petsc-3.5.4 2015-05-23
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