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