Actual source code: toyf.F90
petsc-3.14.6 2021-03-30
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
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*/