xref: /petsc/src/tao/leastsquares/tutorials/chwirut1f.F90 (revision 9b88ac225e01f016352a5f4cd90e158abe5f5675) 
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!
10
11!
12! ----------------------------------------------------------------------
13!
14#include <petsc/finclude/petsctao.h>
15module chwirut1fmodule
16  use petsctao
17  implicit none
18
19  PetscReal t(0:213)
20  PetscReal y(0:213)
21  PetscInt m, n
22
23contains
24
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)
36
37    Tao ta
38    Vec x, f
39    PetscErrorCode ierr
40    PetscInt dummy
41
42    PetscInt i
43    PetscScalar, pointer, dimension(:)  :: x_v, f_v
44
45    ierr = 0
46
47!     Get pointers to vector data
48    PetscCall(VecGetArray(x, x_v, ierr))
49    PetscCall(VecGetArray(f, f_v, ierr))
50
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
55
56!     Restore vectors
57    PetscCall(VecRestoreArray(X, x_v, ierr))
58    PetscCall(VecRestoreArray(F, f_v, ierr))
59
60  end
61
62  subroutine FormStartingPoint(x)
63    Vec x
64    PetscScalar, pointer, dimension(:)  :: x_v
65    PetscErrorCode ierr
66
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
73
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
291
292  end
293end module chwirut1fmodule
294
295program main
296  use chwirut1fmodule
297  implicit none
298! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
299!                   Variable declarations
300! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
301!
302!  See additional variable declarations in the file chwirut1f.h
303
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
316
317!  Initialize TAO and PETSc
318  PetscCallA(PetscInitialize(ierr))
319
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')
323
324!  Initialize problem parameters
325  m = 214
326  n = 3
327
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))
331
332!  The TAO code begins here
333
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
338
339  PetscCallA(TaoSetResidualRoutine(ta, f, FormFunction, 0, ierr))
340
341!  Optional: Set initial guess
342  call InitializeData()
343  call FormStartingPoint(x)
344  PetscCallA(TaoSetSolution(ta, x, ierr))
345
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
359
360!  Free TAO data structures
361  PetscCallA(TaoDestroy(ta, ierr))
362
363!  Free PETSc data structures
364  PetscCallA(VecDestroy(x, ierr))
365  PetscCallA(VecDestroy(f, ierr))
366
367  PetscCallA(PetscFinalize(ierr))
368
369end
370
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*/
381