Actual source code: chwirut1f.F90
1: ! Program usage: mpiexec -n 1 chwirut1f [-help] [all TAO options]
2: !
3: ! Description: This example demonstrates use of the TAO package to solve a
4: ! nonlinear least-squares problem on a single processor. We minimize the
5: ! Chwirut function:
6: ! sum_{i=0}^{n/2-1} ( alpha*(x_{2i+1}-x_{2i}^2)^2 + (1-x_{2i})^2)
7: !
8: ! The C version of this code is test_chwirut1.c
9: !
11: !
12: ! ----------------------------------------------------------------------
13: !
14: module chwirut1fmodule
15: use petsctao
16: #include <petsc/finclude/petsctao.h>
17: PetscReal t(0:213)
18: PetscReal y(0:213)
19: PetscInt m,n
20: end module chwirut1fmodule
22: program main
23: use chwirut1fmodule
24: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
25: ! Variable declarations
26: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
27: !
28: ! See additional variable declarations in the file chwirut1f.h
30: PetscErrorCode ierr ! used to check for functions returning nonzeros
31: Vec x ! solution vector
32: Vec f ! vector of functions
33: Tao tao ! Tao context
34: PetscInt nhist
35: PetscMPIInt size,rank ! number of processes running
36: PetscReal hist(100) ! objective value history
37: PetscReal resid(100)! residual history
38: PetscReal cnorm(100)! cnorm history
39: PetscInt lits(100) ! #ksp history
40: PetscInt oh
41: TaoConvergedReason reason
43: ! Note: Any user-defined Fortran routines (such as FormGradient)
44: ! MUST be declared as external.
46: external FormFunction
48: ! Initialize TAO and PETSc
49: PetscCallA(PetscInitialize(ierr))
51: PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD,size,ierr))
52: PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr))
53: PetscCheckA(size .eq. 1,PETSC_COMM_SELF,PETSC_ERR_WRONG_MPI_SIZE,'This is a uniprocessor example only')
55: ! Initialize problem parameters
56: m = 214
57: n = 3
59: ! Allocate vectors for the solution and gradient
60: PetscCallA(VecCreateSeq(PETSC_COMM_SELF,n,x,ierr))
61: PetscCallA(VecCreateSeq(PETSC_COMM_SELF,m,f,ierr))
63: ! The TAO code begins here
65: ! Create TAO solver
66: PetscCallA(TaoCreate(PETSC_COMM_SELF,tao,ierr))
67: PetscCallA(TaoSetType(tao,TAOPOUNDERS,ierr))
68: ! Set routines for function, gradient, and hessian evaluation
70: PetscCallA(TaoSetResidualRoutine(tao,f,FormFunction,0,ierr))
72: ! Optional: Set initial guess
73: call InitializeData()
74: call FormStartingPoint(x)
75: PetscCallA(TaoSetSolution(tao, x, ierr))
77: ! Check for TAO command line options
78: PetscCallA(TaoSetFromOptions(tao,ierr))
79: oh = 100
80: PetscCallA(TaoSetConvergenceHistory(tao,hist,resid,cnorm,lits,oh,PETSC_TRUE,ierr))
81: ! SOLVE THE APPLICATION
82: PetscCallA(TaoSolve(tao,ierr))
83: PetscCallA(TaoGetConvergenceHistory(tao,nhist,ierr))
84: PetscCallA(TaoGetConvergedReason(tao, reason, ierr))
85: if (reason .le. 0) then
86: print *,'Tao failed.'
87: print *,'Try a different TAO method, adjust some parameters,'
88: print *,'or check the function evaluation routines.'
89: endif
91: ! Free TAO data structures
92: PetscCallA(TaoDestroy(tao,ierr))
94: ! Free PETSc data structures
95: PetscCallA(VecDestroy(x,ierr))
96: PetscCallA(VecDestroy(f,ierr))
98: PetscCallA(PetscFinalize(ierr))
100: end
102: ! --------------------------------------------------------------------
103: ! FormFunction - Evaluates the function f(X) and gradient G(X)
104: !
105: ! Input Parameters:
106: ! tao - the Tao context
107: ! X - input vector
108: ! dummy - not used
109: !
110: ! Output Parameters:
111: ! f - function vector
113: subroutine FormFunction(tao, x, f, dummy, ierr)
114: use chwirut1fmodule
116: Tao tao
117: Vec x,f
118: PetscErrorCode ierr
119: PetscInt dummy
121: PetscInt i
122: PetscScalar, pointer, dimension(:) :: x_v,f_v
124: ierr = 0
126: ! Get pointers to vector data
127: PetscCall(VecGetArrayF90(x,x_v,ierr))
128: PetscCall(VecGetArrayF90(f,f_v,ierr))
130: ! Compute F(X)
131: do i=0,m-1
132: f_v(i+1) = y(i) - exp(-x_v(1)*t(i))/(x_v(2) + x_v(3)*t(i))
133: enddo
135: ! Restore vectors
136: PetscCall(VecRestoreArrayF90(X,x_v,ierr))
137: PetscCall(VecRestoreArrayF90(F,f_v,ierr))
139: end
141: subroutine FormStartingPoint(x)
142: use chwirut1fmodule
144: Vec x
145: PetscScalar, pointer, dimension(:) :: x_v
146: PetscErrorCode ierr
148: PetscCall(VecGetArrayF90(x,x_v,ierr))
149: x_v(1) = 0.15
150: x_v(2) = 0.008
151: x_v(3) = 0.01
152: PetscCall(VecRestoreArrayF90(x,x_v,ierr))
153: end
155: subroutine InitializeData()
156: use chwirut1fmodule
158: integer i
159: i=0
160: y(i) = 92.9000; t(i) = 0.5000; i=i+1
161: y(i) = 78.7000; t(i) = 0.6250; i=i+1
162: y(i) = 64.2000; t(i) = 0.7500; i=i+1
163: y(i) = 64.9000; t(i) = 0.8750; i=i+1
164: y(i) = 57.1000; t(i) = 1.0000; i=i+1
165: y(i) = 43.3000; t(i) = 1.2500; i=i+1
166: y(i) = 31.1000; t(i) = 1.7500; i=i+1
167: y(i) = 23.6000; t(i) = 2.2500; i=i+1
168: y(i) = 31.0500; t(i) = 1.7500; i=i+1
169: y(i) = 23.7750; t(i) = 2.2500; i=i+1
170: y(i) = 17.7375; t(i) = 2.7500; i=i+1
171: y(i) = 13.8000; t(i) = 3.2500; i=i+1
172: y(i) = 11.5875; t(i) = 3.7500; i=i+1
173: y(i) = 9.4125; t(i) = 4.2500; i=i+1
174: y(i) = 7.7250; t(i) = 4.7500; i=i+1
175: y(i) = 7.3500; t(i) = 5.2500; i=i+1
176: y(i) = 8.0250; t(i) = 5.7500; i=i+1
177: y(i) = 90.6000; t(i) = 0.5000; i=i+1
178: y(i) = 76.9000; t(i) = 0.6250; i=i+1
179: y(i) = 71.6000; t(i) = 0.7500; i=i+1
180: y(i) = 63.6000; t(i) = 0.8750; i=i+1
181: y(i) = 54.0000; t(i) = 1.0000; i=i+1
182: y(i) = 39.2000; t(i) = 1.2500; i=i+1
183: y(i) = 29.3000; t(i) = 1.7500; i=i+1
184: y(i) = 21.4000; t(i) = 2.2500; i=i+1
185: y(i) = 29.1750; t(i) = 1.7500; i=i+1
186: y(i) = 22.1250; t(i) = 2.2500; i=i+1
187: y(i) = 17.5125; t(i) = 2.7500; i=i+1
188: y(i) = 14.2500; t(i) = 3.2500; i=i+1
189: y(i) = 9.4500; t(i) = 3.7500; i=i+1
190: y(i) = 9.1500; t(i) = 4.2500; i=i+1
191: y(i) = 7.9125; t(i) = 4.7500; i=i+1
192: y(i) = 8.4750; t(i) = 5.2500; i=i+1
193: y(i) = 6.1125; t(i) = 5.7500; i=i+1
194: y(i) = 80.0000; t(i) = 0.5000; i=i+1
195: y(i) = 79.0000; t(i) = 0.6250; i=i+1
196: y(i) = 63.8000; t(i) = 0.7500; i=i+1
197: y(i) = 57.2000; t(i) = 0.8750; i=i+1
198: y(i) = 53.2000; t(i) = 1.0000; i=i+1
199: y(i) = 42.5000; t(i) = 1.2500; i=i+1
200: y(i) = 26.8000; t(i) = 1.7500; i=i+1
201: y(i) = 20.4000; t(i) = 2.2500; i=i+1
202: y(i) = 26.8500; t(i) = 1.7500; i=i+1
203: y(i) = 21.0000; t(i) = 2.2500; i=i+1
204: y(i) = 16.4625; t(i) = 2.7500; i=i+1
205: y(i) = 12.5250; t(i) = 3.2500; i=i+1
206: y(i) = 10.5375; t(i) = 3.7500; i=i+1
207: y(i) = 8.5875; t(i) = 4.2500; i=i+1
208: y(i) = 7.1250; t(i) = 4.7500; i=i+1
209: y(i) = 6.1125; t(i) = 5.2500; i=i+1
210: y(i) = 5.9625; t(i) = 5.7500; i=i+1
211: y(i) = 74.1000; t(i) = 0.5000; i=i+1
212: y(i) = 67.3000; t(i) = 0.6250; i=i+1
213: y(i) = 60.8000; t(i) = 0.7500; i=i+1
214: y(i) = 55.5000; t(i) = 0.8750; i=i+1
215: y(i) = 50.3000; t(i) = 1.0000; i=i+1
216: y(i) = 41.0000; t(i) = 1.2500; i=i+1
217: y(i) = 29.4000; t(i) = 1.7500; i=i+1
218: y(i) = 20.4000; t(i) = 2.2500; i=i+1
219: y(i) = 29.3625; t(i) = 1.7500; i=i+1
220: y(i) = 21.1500; t(i) = 2.2500; i=i+1
221: y(i) = 16.7625; t(i) = 2.7500; i=i+1
222: y(i) = 13.2000; t(i) = 3.2500; i=i+1
223: y(i) = 10.8750; t(i) = 3.7500; i=i+1
224: y(i) = 8.1750; t(i) = 4.2500; i=i+1
225: y(i) = 7.3500; t(i) = 4.7500; i=i+1
226: y(i) = 5.9625; t(i) = 5.2500; i=i+1
227: y(i) = 5.6250; t(i) = 5.7500; i=i+1
228: y(i) = 81.5000; t(i) = .5000; i=i+1
229: y(i) = 62.4000; t(i) = .7500; i=i+1
230: y(i) = 32.5000; t(i) = 1.5000; i=i+1
231: y(i) = 12.4100; t(i) = 3.0000; i=i+1
232: y(i) = 13.1200; t(i) = 3.0000; i=i+1
233: y(i) = 15.5600; t(i) = 3.0000; i=i+1
234: y(i) = 5.6300; t(i) = 6.0000; i=i+1
235: y(i) = 78.0000; t(i) = .5000; i=i+1
236: y(i) = 59.9000; t(i) = .7500; i=i+1
237: y(i) = 33.2000; t(i) = 1.5000; i=i+1
238: y(i) = 13.8400; t(i) = 3.0000; i=i+1
239: y(i) = 12.7500; t(i) = 3.0000; i=i+1
240: y(i) = 14.6200; t(i) = 3.0000; i=i+1
241: y(i) = 3.9400; t(i) = 6.0000; i=i+1
242: y(i) = 76.8000; t(i) = .5000; i=i+1
243: y(i) = 61.0000; t(i) = .7500; i=i+1
244: y(i) = 32.9000; t(i) = 1.5000; i=i+1
245: y(i) = 13.8700; t(i) = 3.0000; i=i+1
246: y(i) = 11.8100; t(i) = 3.0000; i=i+1
247: y(i) = 13.3100; t(i) = 3.0000; i=i+1
248: y(i) = 5.4400; t(i) = 6.0000; i=i+1
249: y(i) = 78.0000; t(i) = .5000; i=i+1
250: y(i) = 63.5000; t(i) = .7500; i=i+1
251: y(i) = 33.8000; t(i) = 1.5000; i=i+1
252: y(i) = 12.5600; t(i) = 3.0000; i=i+1
253: y(i) = 5.6300; t(i) = 6.0000; i=i+1
254: y(i) = 12.7500; t(i) = 3.0000; i=i+1
255: y(i) = 13.1200; t(i) = 3.0000; i=i+1
256: y(i) = 5.4400; t(i) = 6.0000; i=i+1
257: y(i) = 76.8000; t(i) = .5000; i=i+1
258: y(i) = 60.0000; t(i) = .7500; i=i+1
259: y(i) = 47.8000; t(i) = 1.0000; i=i+1
260: y(i) = 32.0000; t(i) = 1.5000; i=i+1
261: y(i) = 22.2000; t(i) = 2.0000; i=i+1
262: y(i) = 22.5700; t(i) = 2.0000; i=i+1
263: y(i) = 18.8200; t(i) = 2.5000; i=i+1
264: y(i) = 13.9500; t(i) = 3.0000; i=i+1
265: y(i) = 11.2500; t(i) = 4.0000; i=i+1
266: y(i) = 9.0000; t(i) = 5.0000; i=i+1
267: y(i) = 6.6700; t(i) = 6.0000; i=i+1
268: y(i) = 75.8000; t(i) = .5000; i=i+1
269: y(i) = 62.0000; t(i) = .7500; i=i+1
270: y(i) = 48.8000; t(i) = 1.0000; i=i+1
271: y(i) = 35.2000; t(i) = 1.5000; i=i+1
272: y(i) = 20.0000; t(i) = 2.0000; i=i+1
273: y(i) = 20.3200; t(i) = 2.0000; i=i+1
274: y(i) = 19.3100; t(i) = 2.5000; i=i+1
275: y(i) = 12.7500; t(i) = 3.0000; i=i+1
276: y(i) = 10.4200; t(i) = 4.0000; i=i+1
277: y(i) = 7.3100; t(i) = 5.0000; i=i+1
278: y(i) = 7.4200; t(i) = 6.0000; i=i+1
279: y(i) = 70.5000; t(i) = .5000; i=i+1
280: y(i) = 59.5000; t(i) = .7500; i=i+1
281: y(i) = 48.5000; t(i) = 1.0000; i=i+1
282: y(i) = 35.8000; t(i) = 1.5000; i=i+1
283: y(i) = 21.0000; t(i) = 2.0000; i=i+1
284: y(i) = 21.6700; t(i) = 2.0000; i=i+1
285: y(i) = 21.0000; t(i) = 2.5000; i=i+1
286: y(i) = 15.6400; t(i) = 3.0000; i=i+1
287: y(i) = 8.1700; t(i) = 4.0000; i=i+1
288: y(i) = 8.5500; t(i) = 5.0000; i=i+1
289: y(i) = 10.1200; t(i) = 6.0000; i=i+1
290: y(i) = 78.0000; t(i) = .5000; i=i+1
291: y(i) = 66.0000; t(i) = .6250; i=i+1
292: y(i) = 62.0000; t(i) = .7500; i=i+1
293: y(i) = 58.0000; t(i) = .8750; i=i+1
294: y(i) = 47.7000; t(i) = 1.0000; i=i+1
295: y(i) = 37.8000; t(i) = 1.2500; i=i+1
296: y(i) = 20.2000; t(i) = 2.2500; i=i+1
297: y(i) = 21.0700; t(i) = 2.2500; i=i+1
298: y(i) = 13.8700; t(i) = 2.7500; i=i+1
299: y(i) = 9.6700; t(i) = 3.2500; i=i+1
300: y(i) = 7.7600; t(i) = 3.7500; i=i+1
301: y(i) = 5.4400; t(i) = 4.2500; i=i+1
302: y(i) = 4.8700; t(i) = 4.7500; i=i+1
303: y(i) = 4.0100; t(i) = 5.2500; i=i+1
304: y(i) = 3.7500; t(i) = 5.7500; i=i+1
305: y(i) = 24.1900; t(i) = 3.0000; i=i+1
306: y(i) = 25.7600; t(i) = 3.0000; i=i+1
307: y(i) = 18.0700; t(i) = 3.0000; i=i+1
308: y(i) = 11.8100; t(i) = 3.0000; i=i+1
309: y(i) = 12.0700; t(i) = 3.0000; i=i+1
310: y(i) = 16.1200; t(i) = 3.0000; i=i+1
311: y(i) = 70.8000; t(i) = .5000; i=i+1
312: y(i) = 54.7000; t(i) = .7500; i=i+1
313: y(i) = 48.0000; t(i) = 1.0000; i=i+1
314: y(i) = 39.8000; t(i) = 1.5000; i=i+1
315: y(i) = 29.8000; t(i) = 2.0000; i=i+1
316: y(i) = 23.7000; t(i) = 2.5000; i=i+1
317: y(i) = 29.6200; t(i) = 2.0000; i=i+1
318: y(i) = 23.8100; t(i) = 2.5000; i=i+1
319: y(i) = 17.7000; t(i) = 3.0000; i=i+1
320: y(i) = 11.5500; t(i) = 4.0000; i=i+1
321: y(i) = 12.0700; t(i) = 5.0000; i=i+1
322: y(i) = 8.7400; t(i) = 6.0000; i=i+1
323: y(i) = 80.7000; t(i) = .5000; i=i+1
324: y(i) = 61.3000; t(i) = .7500; i=i+1
325: y(i) = 47.5000; t(i) = 1.0000; i=i+1
326: y(i) = 29.0000; t(i) = 1.5000; i=i+1
327: y(i) = 24.0000; t(i) = 2.0000; i=i+1
328: y(i) = 17.7000; t(i) = 2.5000; i=i+1
329: y(i) = 24.5600; t(i) = 2.0000; i=i+1
330: y(i) = 18.6700; t(i) = 2.5000; i=i+1
331: y(i) = 16.2400; t(i) = 3.0000; i=i+1
332: y(i) = 8.7400; t(i) = 4.0000; i=i+1
333: y(i) = 7.8700; t(i) = 5.0000; i=i+1
334: y(i) = 8.5100; t(i) = 6.0000; i=i+1
335: y(i) = 66.7000; t(i) = .5000; i=i+1
336: y(i) = 59.2000; t(i) = .7500; i=i+1
337: y(i) = 40.8000; t(i) = 1.0000; i=i+1
338: y(i) = 30.7000; t(i) = 1.5000; i=i+1
339: y(i) = 25.7000; t(i) = 2.0000; i=i+1
340: y(i) = 16.3000; t(i) = 2.5000; i=i+1
341: y(i) = 25.9900; t(i) = 2.0000; i=i+1
342: y(i) = 16.9500; t(i) = 2.5000; i=i+1
343: y(i) = 13.3500; t(i) = 3.0000; i=i+1
344: y(i) = 8.6200; t(i) = 4.0000; i=i+1
345: y(i) = 7.2000; t(i) = 5.0000; i=i+1
346: y(i) = 6.6400; t(i) = 6.0000; i=i+1
347: y(i) = 13.6900; t(i) = 3.0000; i=i+1
348: y(i) = 81.0000; t(i) = .5000; i=i+1
349: y(i) = 64.5000; t(i) = .7500; i=i+1
350: y(i) = 35.5000; t(i) = 1.5000; i=i+1
351: y(i) = 13.3100; t(i) = 3.0000; i=i+1
352: y(i) = 4.8700; t(i) = 6.0000; i=i+1
353: y(i) = 12.9400; t(i) = 3.0000; i=i+1
354: y(i) = 5.0600; t(i) = 6.0000; i=i+1
355: y(i) = 15.1900; t(i) = 3.0000; i=i+1
356: y(i) = 14.6200; t(i) = 3.0000; i=i+1
357: y(i) = 15.6400; t(i) = 3.0000; i=i+1
358: y(i) = 25.5000; t(i) = 1.7500; i=i+1
359: y(i) = 25.9500; t(i) = 1.7500; i=i+1
360: y(i) = 81.7000; t(i) = .5000; i=i+1
361: y(i) = 61.6000; t(i) = .7500; i=i+1
362: y(i) = 29.8000; t(i) = 1.7500; i=i+1
363: y(i) = 29.8100; t(i) = 1.7500; i=i+1
364: y(i) = 17.1700; t(i) = 2.7500; i=i+1
365: y(i) = 10.3900; t(i) = 3.7500; i=i+1
366: y(i) = 28.4000; t(i) = 1.7500; i=i+1
367: y(i) = 28.6900; t(i) = 1.7500; i=i+1
368: y(i) = 81.3000; t(i) = .5000; i=i+1
369: y(i) = 60.9000; t(i) = .7500; i=i+1
370: y(i) = 16.6500; t(i) = 2.7500; i=i+1
371: y(i) = 10.0500; t(i) = 3.7500; i=i+1
372: y(i) = 28.9000; t(i) = 1.7500; i=i+1
373: y(i) = 28.9500; t(i) = 1.7500; i=i+1
375: end
377: !/*TEST
378: !
379: ! build:
380: ! requires: !complex
381: !
382: ! test:
383: ! args: -tao_monitor_short -tao_max_it 100 -tao_type pounders -tao_gatol 1.e-5
384: ! requires: !single
385: !
386: !TEST*/