V3FIT
nelmin.f
1  subroutine nelmin ( fn, n, start, xmin, ynewlo, reqmin, step,
2  & konvge, kcount, icount, numres, ifault )
3 
4 !*****************************************************************************80
5 !
6 !! NELMIN minimizes a function using the Nelder-Mead algorithm.
7 !
8 ! Discussion:
9 !
10 ! This routine seeks the minimum value of a user-specified function.
11 !
12 ! Simplex function minimisation procedure due to Nelder and Mead (1965),
13 ! as implemented by O'Neill(1971, Appl.Statist. 20, 338-45), with
14 ! subsequent comments by Chambers+Ertel(1974, 23, 250-1), Benyon(1976,
15 ! 25, 97) and Hill(1978, 27, 380-2)
16 !
17 ! The function to be minimized must be defined by a function of
18 ! the form
19 !
20 ! function fn ( x )
21 ! real ( kind = 8 ) fn
22 ! real ( kind = 8 ) x(*)
23 !
24 ! and the name of this subroutine must be declared EXTERNAL in the
25 ! calling routine and passed as the argument FN.
26 !
27 ! This routine does not include a termination test using the
28 ! fitting of a quadratic surface.
29 !
30 ! Licensing:
31 !
32 ! This code is distributed under the GNU LGPL license.
33 !
34 ! Modified:
35 !
36 ! 27 February 2008
37 !
38 ! Author:
39 !
40 ! Original FORTRAN77 version by R ONeill.
41 ! FORTRAN90 version by John Burkardt.
42 !
43 ! Reference:
44 !
45 ! John Nelder, Roger Mead,
46 ! A simplex method for function minimization,
47 ! Computer Journal,
48 ! Volume 7, 1965, pages 308-313.
49 !
50 ! R ONeill,
51 ! Algorithm AS 47:
52 ! Function Minimization Using a Simplex Procedure,
53 ! Applied Statistics,
54 ! Volume 20, Number 3, 1971, pages 338-345.
55 !
56 ! Parameters:
57 !
58 ! Input, external FN, the name of the function which evaluates
59 ! the function to be minimized.
60 !
61 ! Input, integer ( kind = 4 ) N, the number of variables.
62 ! 0 < N is required.
63 !
64 ! Input/output, real ( kind = 8 ) START(N). On input, a starting point
65 ! for the iteration. On output, this data may have been overwritten.
66 !
67 ! Output, real ( kind = 8 ) XMIN(N), the coordinates of the point which
68 ! is estimated to minimize the function.
69 !
70 ! Output, real ( kind = 8 ) YNEWLO, the minimum value of the function.
71 !
72 ! Input, real ( kind = 8 ) REQMIN, the terminating limit for the variance
73 ! of the function values. 0 < REQMIN is required.
74 !
75 ! Input, real ( kind = 8 ) STEP(N), determines the size and shape of the
76 ! initial simplex. The relative magnitudes of its elements should reflect
77 ! the units of the variables.
78 !
79 ! Input, integer ( kind = 4 ) KONVGE, the convergence check is carried out
80 ! every KONVGE iterations. 0 < KONVGE is required.
81 !
82 ! Input, integer ( kind = 4 ) KCOUNT, the maximum number of function
83 ! evaluations.
84 !
85 ! Output, integer ( kind = 4 ) ICOUNT, the number of function evaluations
86 ! used.
87 !
88 ! Output, integer ( kind = 4 ) NUMRES, the number of restarts.
89 !
90 ! Output, integer ( kind = 4 ) IFAULT, error indicator.
91 ! 0, no errors detected.
92 ! 1, REQMIN, N, or KONVGE has an illegal value.
93 ! 2, iteration terminated because KCOUNT was exceeded without convergence.
94 !
95  implicit none
96 
97  integer ( kind = 4 ) n
98 
99  real ( kind = 8 ), parameter :: ccoeff = 0.5d+00
100  real ( kind = 8 ) del
101  real ( kind = 8 ), parameter :: ecoeff = 2.0d+00
102  real ( kind = 8 ), parameter :: eps = 0.001d+00
103  real ( kind = 8 ), external :: fn
104  integer ( kind = 4 ) i
105  integer ( kind = 4 ) icount
106  integer ( kind = 4 ) ifault
107  integer ( kind = 4 ) ihi
108  integer ( kind = 4 ) ilo
109  integer ( kind = 4 ) j
110  integer ( kind = 4 ) jcount
111  integer ( kind = 4 ) kcount
112  integer ( kind = 4 ) konvge
113  integer ( kind = 4 ) l
114  integer ( kind = 4 ) numres
115  real ( kind = 8 ) p(n,n+1)
116  real ( kind = 8 ) p2star(n)
117  real ( kind = 8 ) pbar(n)
118  real ( kind = 8 ) pstar(n)
119  real ( kind = 8 ), parameter :: rcoeff = 1.0d+00
120  real ( kind = 8 ) reqmin
121  real ( kind = 8 ) rq
122  real ( kind = 8 ) start(n)
123  real ( kind = 8 ) step(n)
124  real ( kind = 8 ) x
125  real ( kind = 8 ) xmin(n)
126  real ( kind = 8 ) y(n+1)
127  real ( kind = 8 ) y2star
128  real ( kind = 8 ) ylo
129  real ( kind = 8 ) ynewlo
130  real ( kind = 8 ) ystar
131  real ( kind = 8 ) z
132 !
133 ! Check the input parameters.
134 !
135  if ( reqmin <= 0.0d+00 ) then
136  ifault = 1
137  return
138  end if
139 
140  if ( n < 1 ) then
141  ifault = 1
142  return
143  end if
144 
145  if ( konvge < 1 ) then
146  ifault = 1
147  return
148  end if
149 !
150 ! Initialization.
151 !
152  icount = 0
153  numres = 0
154  jcount = konvge
155  del = 1.0d+00
156  rq = reqmin * real( n, kind = 8 )
157 !
158 ! Initial or restarted loop.
159 !
160  do
161 
162  p(1:n,n+1) = start(1:n)
163  y(n+1) = fn( start )
164  icount = icount + 1
165 !
166 ! Define the initial simplex.
167 !
168  do j = 1, n
169  x = start(j)
170  start(j) = start(j) + step(j) * del
171  p(1:n,j) = start(1:n)
172  y(j) = fn( start )
173  icount = icount + 1
174  start(j) = x
175  end do
176 !
177 ! Find highest and lowest Y values. YNEWLO = Y(IHI) indicates
178 ! the vertex of the simplex to be replaced.
179 !
180  ilo = minloc( y(1:n+1), 1 )
181  ylo = y(ilo)
182 !
183 ! Inner loop.
184 !
185  do while ( icount < kcount )
186 !
187 ! YNEWLO is, of course, the HIGHEST value???
188 !
189  ihi = maxloc( y(1:n+1), 1 )
190  ynewlo = y(ihi)
191 !
192 ! Calculate PBAR, the centroid of the simplex vertices
193 ! excepting the vertex with Y value YNEWLO.
194 !
195  do i = 1, n
196  pbar(i) = (sum( p(i,1:n+1) ) - p(i,ihi)) / real(n, kind = 8)
197  end do
198 !
199 ! Reflection through the centroid.
200 !
201  pstar(1:n) = pbar(1:n) + rcoeff * ( pbar(1:n) - p(1:n,ihi) )
202  ystar = fn( pstar )
203  icount = icount + 1
204 !
205 ! Successful reflection, so extension.
206 !
207  if ( ystar < ylo ) then
208 
209  p2star(1:n) = pbar(1:n) + ecoeff * ( pstar(1:n) - pbar(1:n) )
210  y2star = fn( p2star )
211  icount = icount + 1
212 !
213 ! Retain extension or contraction.
214 !
215  if ( ystar < y2star ) then
216  p(1:n,ihi) = pstar(1:n)
217  y(ihi) = ystar
218  else
219  p(1:n,ihi) = p2star(1:n)
220  y(ihi) = y2star
221  end if
222 !
223 ! No extension.
224 !
225  else
226 
227  l = 0
228  do i = 1, n + 1
229  if ( ystar < y(i) ) then
230  l = l + 1
231  end if
232  end do
233 
234  if ( 1 < l ) then
235 
236  p(1:n,ihi) = pstar(1:n)
237  y(ihi) = ystar
238 !
239 ! Contraction on the Y(IHI) side of the centroid.
240 !
241  else if ( l == 0 ) then
242 
243  p2star(1:n) = pbar(1:n) + ccoeff * ( p(1:n,ihi) - pbar(1:n) )
244  y2star = fn( p2star )
245  icount = icount + 1
246 !
247 ! Contract the whole simplex.
248 !
249  if ( y(ihi) < y2star ) then
250 
251  do j = 1, n + 1
252  p(1:n,j) = ( p(1:n,j) + p(1:n,ilo) ) * 0.5d+00
253  xmin(1:n) = p(1:n,j)
254  y(j) = fn( xmin )
255  icount = icount + 1
256  end do
257 
258  ilo = minloc( y(1:n+1), 1 )
259  ylo = y(ilo)
260 
261  cycle
262 !
263 ! Retain contraction.
264 !
265  else
266  p(1:n,ihi) = p2star(1:n)
267  y(ihi) = y2star
268  end if
269 !
270 ! Contraction on the reflection side of the centroid.
271 !
272  else if ( l == 1 ) then
273 
274  p2star(1:n) = pbar(1:n) + ccoeff * ( pstar(1:n) - pbar(1:n) )
275  y2star = fn( p2star )
276  icount = icount + 1
277 !
278 ! Retain reflection?
279 !
280  if ( y2star <= ystar ) then
281  p(1:n,ihi) = p2star(1:n)
282  y(ihi) = y2star
283  else
284  p(1:n,ihi) = pstar(1:n)
285  y(ihi) = ystar
286  end if
287 
288  end if
289 
290  end if
291 !
292 ! Check if YLO improved.
293 !
294  if ( y(ihi) < ylo ) then
295  ylo = y(ihi)
296  ilo = ihi
297  end if
298 
299  jcount = jcount - 1
300 
301  if ( 0 < jcount ) then
302  cycle
303  end if
304 !
305 ! Check to see if minimum reached.
306 !
307  if ( icount <= kcount ) then
308 
309  jcount = konvge
310 
311  x = sum( y(1:n+1) ) / real( n + 1, kind = 8 )
312  z = sum( ( y(1:n+1) - x )**2 )
313 
314  if ( z <= rq ) then
315  exit
316  end if
317 
318  end if
319 
320  end do
321 !
322 ! Factorial tests to check that YNEWLO is a local minimum.
323 !
324  xmin(1:n) = p(1:n,ilo)
325  ynewlo = y(ilo)
326 
327  if ( kcount < icount ) then
328  ifault = 2
329  exit
330  end if
331 
332  ifault = 0
333 
334  do i = 1, n
335  del = step(i) * eps
336  xmin(i) = xmin(i) + del
337  z = fn( xmin )
338  icount = icount + 1
339  if ( z < ynewlo ) then
340  ifault = 2
341  exit
342  end if
343  xmin(i) = xmin(i) - del - del
344  z = fn( xmin )
345  icount = icount + 1
346  if ( z < ynewlo ) then
347  ifault = 2
348  exit
349  end if
350  xmin(i) = xmin(i) + del
351  end do
352 
353  if ( ifault == 0 ) then
354  exit
355  end if
356 !
357 ! Restart the procedure.
358 !
359  start(1:n) = xmin(1:n)
360  del = eps
361  numres = numres + 1
362 
363  end do
364 
365  return
366  end