/[ascend]/trunk/minos54/mi60srch.f
ViewVC logotype

Contents of /trunk/minos54/mi60srch.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1 - (show annotations) (download)
Fri Oct 29 20:54:12 2004 UTC (19 years, 7 months ago) by aw0a
File size: 87614 byte(s)
Setting up web subdirectory in repository
1 ****************************************************************
2 *
3 * file mi60srch fortran.
4 *
5 * m6dmmy m6fcon m6fobj m6fun m6fun1 m6grd m6grd1
6 * m6dobj m6dcon m6srch srchc srchq
7 *
8 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
9
10 subroutine m6dmmy( n, g )
11
12 implicit double precision (a-h,o-z)
13 double precision g(n)
14
15 * ------------------------------------------------------------------
16 * m6dmmy sets the elements of g to a dummy value.
17 * When funcon or funobj are subsequently called, we can
18 * then determine which elements of g are known.
19 * ------------------------------------------------------------------
20
21 common /m8diff/ difint(2),gdummy,lderiv,lvldif,knowng(2)
22
23 call dload ( n, gdummy, g, 1 )
24
25 * end of m6dmmy
26 end
27
28 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
29
30 subroutine m6fcon( modefg, nncon, nnjac, njac, f, g,
31 $ ne, nka, ha, ka,
32 $ x, z, nwcore )
33
34 implicit double precision (a-h,o-z)
35 integer*4 ha(ne)
36 integer ka(nka)
37 double precision f(nncon), g(njac), x(nnjac), z(nwcore)
38
39 * ------------------------------------------------------------------
40 * m6fcon calls the user-written routine funcon to evaluate
41 * the nonlinear constraints and possibly their gradients.
42 *
43 * 08-Jun-1987: Rewritten for GAMS during Tony Brooke's visit.
44 * 20-Mar-1988: Tony found it hadn't been installed (gasp!).
45 * 12 Sep 1991: More if-then-elses.
46 * 19 Nov 1991: ha, ka passed in as parameters to help with minoss.
47 * 11 Jul 1992: Clock 4 added.
48 * ------------------------------------------------------------------
49
50 common /m1file/ iread,iprint,isumm
51 parameter ( ntime = 5 )
52 common /m1tim / tlast(ntime), tsum(ntime), numt(ntime), ltime
53 common /m3len / m ,n ,nb ,nscl
54 common /m3loc / lascal,lbl ,lbu ,lbbl ,lbbu ,
55 $ lhrtyp,lhs ,lkb
56 common /m3scal/ sclobj,scltol,lscale
57 common /m5loc / lpi ,lpi2 ,lw ,lw2 ,
58 $ lx ,lx2 ,ly ,ly2 ,
59 $ lgsub ,lgsub2,lgrd ,lgrd2 ,
60 $ lr ,lrg ,lrg2 ,lxn
61 common /m5log1/ idebug,ierr,lprint
62 logical prnt0 ,prnt1 ,summ0 ,summ1 ,newhed
63 common /m5log4/ prnt0 ,prnt1 ,summ0 ,summ1 ,newhed
64 common /m8loc / lfcon ,lfcon2,lfdif ,lfdif2,lfold ,
65 $ lblslk,lbuslk,lxlam ,lrhs ,
66 $ lgcon ,lgcon2,lxdif ,lxold
67 common /m8diff/ difint(2),gdummy,lderiv,lvldif,knowng(2)
68 common /m8func/ nfcon(4),nfobj(4),nprob,nstat1,nstat2
69
70 logical scaled
71
72 if (ltime .ge. 2) call m1time( 4,0 )
73 scaled = lscale .eq. 2
74 mode = modefg
75 nfcon(1) = nfcon(1) + 1
76 if (mode .eq. 2) nfcon(2) = nfcon(2) + 1
77
78 * ------------------------------------------------------------------
79 * On first entry, count how many Jacobian elements are known.
80 * ------------------------------------------------------------------
81 if (nstat2 .eq. 1) then
82 call m6dmmy( njac, g )
83 call funcon( mode, nncon, nnjac, njac, x, f, g,
84 $ nstat2, nprob, z, nwcore )
85 nstat2 = 0
86 if (mode .lt. 0) go to 900
87
88 ngrad = 0
89 do 100 l = 1, njac
90 if (g(l) .ne. gdummy) ngrad = ngrad + 1
91 100 continue
92
93 knowng(2) = ngrad
94 if (iprint .gt. 0) then
95 write(iprint, 1100)
96 write(iprint, 2000) ngrad, njac
97 end if
98 if (isumm .gt. 0) then
99 write(isumm , 2000) ngrad, njac
100 end if
101 if (ngrad .eq. njac) go to 990
102
103 * Some Jacobian elements are missing.
104 * If the Jacobian is supposedly known (lderiv .ge. 2), the
105 * missing elements are assumed to be constant, and must be
106 * restored. They have been saved from the initial a(*) by
107 * m8augl( 1, ... ) and saved in gcon2.
108
109 if (lderiv .ge. 2) call dcopy ( njac, z(lgcon2), 1, g, 1 )
110 nfcon(2) = nfcon(2) + 1
111 go to 400
112 end if
113
114 * Test for last entry.
115 * Scaling will have been disabled, so it's ok to
116 * skip the subsequent part.
117
118 if (nstat2 .ge. 2) then
119 if (iprint .gt. 0) write(iprint, 3000) nstat2
120 if (isumm .gt. 0) write(isumm , 3000) nstat2
121 go to 400
122 end if
123
124 * ------------------------------------------------------------------
125 * Normal entry.
126 * If the Jacobian is known and scaled, and if there are some
127 * constant elements saved in gcon2, we have to copy them into g.
128 * ------------------------------------------------------------------
129
130 if (scaled) then
131 call dcopy ( nnjac, x , 1, z(lx2), 1 )
132 call ddscl ( nnjac, z(lascal), 1, x, 1 )
133 if (lderiv .ge. 2 .and. knowng(2) .lt. njac
134 $ .and. modefg .eq. 2 ) then
135 call dcopy ( njac, z(lgcon2), 1, g, 1 )
136 end if
137 end if
138
139 400 continue
140 call funcon( mode, nncon, nnjac, njac, x, f, g,
141 $ nstat2, nprob, z, nwcore )
142
143 if (scaled) then
144 call dcopy ( nnjac, z(lx2) , 1, x, 1 )
145 call dddiv ( nncon, z(lascal+n), 1, f, 1 )
146 if (modefg .eq. 2) then
147 call m8sclj( 2, nncon, nnjac, njac, n, nb, ne, nka,
148 $ z(lascal), ha, ka, g )
149 end if
150 end if
151
152 if (mode .lt. 0) go to 900
153 go to 990
154
155 * -----------------------
156 * The user wants to stop.
157 * -----------------------
158 900 ierr = 6
159 call m1page( 2 )
160 if (iprint .gt. 0) write(iprint, 1000)
161 if (isumm .gt. 0) write(isumm , 1000)
162
163 990 if (ltime .ge. 2) call m1time(-4,0 )
164 return
165
166 1000 format(' EXIT -- Termination requested by User',
167 $ ' in subroutine funcon')
168 1100 format(//)
169 2000 format(' funcon sets', i8, ' out of', i8,
170 $ ' constraint gradients.')
171 3000 format(/ ' funcon called with nstate =', i4)
172
173 * end of m6fcon
174 end
175
176 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
177
178 subroutine m6fobj( modefg, n, f, g, x, z, nwcore )
179
180 implicit double precision (a-h,o-z)
181 double precision g(n), x(n), z(nwcore)
182
183 * ------------------------------------------------------------------
184 * m6fobj calls the user-written routine funobj to evaluate
185 * the nonlinear objective function and possibly its gradient.
186 * 12 Sep 1991: More if-then-elses.
187 * 11 Jul 1992: Clock 5 added.
188 * ------------------------------------------------------------------
189
190 common /m1file/ iread,iprint,isumm
191 parameter ( ntime = 5 )
192 common /m1tim / tlast(ntime), tsum(ntime), numt(ntime), ltime
193 common /m3loc / lascal,lbl ,lbu ,lbbl ,lbbu ,
194 $ lhrtyp,lhs ,lkb
195 common /m3scal/ sclobj,scltol,lscale
196 common /m5loc / lpi ,lpi2 ,lw ,lw2 ,
197 $ lx ,lx2 ,ly ,ly2 ,
198 $ lgsub ,lgsub2,lgrd ,lgrd2 ,
199 $ lr ,lrg ,lrg2 ,lxn
200 common /m5log1/ idebug,ierr,lprint
201 logical prnt0 ,prnt1 ,summ0 ,summ1 ,newhed
202 common /m5log4/ prnt0 ,prnt1 ,summ0 ,summ1 ,newhed
203 common /m8diff/ difint(2),gdummy,lderiv,lvldif,knowng(2)
204 common /m8func/ nfcon(4),nfobj(4),nprob,nstat1,nstat2
205
206 logical scaled
207
208 if (ltime .ge. 2) call m1time( 5,0 )
209 scaled = lscale .eq. 2
210 mode = modefg
211 nfobj(1) = nfobj(1) + 1
212 if (mode .eq. 2) nfobj(2) = nfobj(2) + 1
213
214 * Test for first or last entry.
215
216 if (nstat1 .eq. 1) then
217 call m6dmmy( n, g )
218 else if (nstat1 .ge. 2) then
219 if (iprint .gt. 0) write(iprint, 3000) nstat1
220 if (isumm .gt. 0) write(isumm , 3000) nstat1
221 end if
222
223 * Normal entry.
224
225 if (scaled) then
226 call dcopy ( n, x , 1, z(lx2), 1 )
227 call ddscl ( n, z(lascal), 1, x, 1 )
228 end if
229
230 call funobj( mode, n, x, f, g, nstat1, nprob, z, nwcore )
231
232 if (scaled) then
233 call dcopy ( n, z(lx2), 1, x, 1 )
234 if (modefg .eq. 2) then
235 call m8sclj( 1, 1, n, n, 1, n, 1, 1, z(lascal), z, z, g )
236 end if
237 end if
238
239 if (mode .lt. 0) go to 900
240
241 * ------------------------------------------------------------------
242 * On first entry, count components to be estimated by differences.
243 * ------------------------------------------------------------------
244 if (nstat1 .eq. 1) then
245 nstat1 = 0
246 ngrad = 0
247 do 100 l = 1, n
248 if (g(l) .ne. gdummy) ngrad = ngrad + 1
249 100 continue
250
251 knowng(1) = ngrad
252 if (iprint .gt. 0) then
253 write(iprint, 1100)
254 write(iprint, 2000) ngrad, n
255 end if
256 if (isumm .gt. 0) then
257 write(isumm , 2000) ngrad, n
258 end if
259
260 if (ngrad .lt. n) then
261 if (lderiv .eq. 1 .or. lderiv .eq. 3) then
262 lderiv = lderiv - 1
263 if (iprint .gt. 0) write(iprint, 2200) lderiv
264 if (isumm .gt. 0) write(isumm , 2200) lderiv
265 end if
266 end if
267 end if
268
269 go to 990
270
271 * -----------------------
272 * The user wants to stop.
273 * -----------------------
274 900 ierr = 6
275 call m1page( 2 )
276 if (iprint .gt. 0) write(iprint, 1000) nfobj(1)
277 if (isumm .gt. 0) write(isumm , 1000) nfobj(1)
278
279 990 if (ltime .ge. 2) call m1time(-5,0 )
280 return
281
282 1000 format(' EXIT -- Termination requested by User',
283 $ ' in subroutine funobj after', i8, ' calls')
284 1100 format(//)
285 2000 format(' funobj sets', i8, ' out of', i8,
286 $ ' objective gradients.')
287 2200 format(/ ' Derivative Level now reduced to', i3)
288 3000 format(/ ' funobj called with nstate =', i4)
289
290 * end of m6fobj
291 end
292
293 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
294
295 subroutine m6fun ( mode, modefg, n, nb, ms, fsub,
296 $ ne, nka, a, ha, ka,
297 $ x, xn, z, nwcore )
298
299 implicit double precision (a-h,o-z)
300 integer*4 ha(ne)
301 integer ka(nka)
302 double precision a(ne), x(ms), xn(nb), z(nwcore)
303
304 * ------------------------------------------------------------------
305 * m6fun is called mainly by the linesearch, m6srch,
306 * to compute the subproblem function value fsub.
307 * It is also called by m5frmc when a feasible point is first found,
308 * and by m6rgit when the linesearch says the max step is too small.
309 *
310 * If mode = 0, m6fun sets objective and constraint gradients to
311 * dummy values (if necessary) in preparation for calls with mode=1.
312 *
313 * If mode = 1, m6fun returns the subproblem function value in fsub
314 * for the current value of the basic and superbasic vars in x.
315 * The user functions funcon and/or funobj are called
316 * via m6fcon and m6fobj, using modefg to control the gradients
317 * as follows:
318 *
319 * If modefg = 2, gradients will be requested.
320 * If modefg = 0, gradients will NOT be requested. This is used
321 * by m6srch during a function-only linesearch. An extra call
322 * with modefg = 2 is needed once the linesearch finds a good
323 * point, but this may be cheaper than getting gradients every
324 * time, e.g. if the user estimates gradients by differencing.
325 *
326 * 12 Sep 1991: More if-then-elses.
327 * 02 Oct 1991: modefg added following request from Hilary Hoynes,
328 * Hoover Institution (Hoynes@Hoover.bitnet, 723-9511).
329 * 14 Oct 1991: a, ha, ka added as parameters.
330 * ------------------------------------------------------------------
331
332 common /m3loc / lascal,lbl ,lbu ,lbbl ,lbbu ,
333 $ lhrtyp,lhs ,lkb
334 common /m5len / maxr ,maxs ,mbs ,nn ,nn0 ,nr ,nx
335 common /m5loc / lpi ,lpi2 ,lw ,lw2 ,
336 $ lx ,lx2 ,ly ,ly2 ,
337 $ lgsub ,lgsub2,lgrd ,lgrd2 ,
338 $ lr ,lrg ,lrg2 ,lxn
339 common /m7len / fobj ,fobj2 ,nnobj ,nnobj0
340 common /m7loc / lgobj ,lgobj2
341 common /m8len / njac ,nncon ,nncon0,nnjac
342 common /m8loc / lfcon ,lfcon2,lfdif ,lfdif2,lfold ,
343 $ lblslk,lbuslk,lxlam ,lrhs ,
344 $ lgcon ,lgcon2,lxdif ,lxold
345 common /m8al1 / penpar,rowtol,ncom,nden,nlag,nmajor,nminor
346 common /m8diff/ difint(2),gdummy,lderiv,lvldif,knowng(2)
347
348 logical grdcon, grdobj
349
350 if (mode .eq. 0) then
351 * ---------------------------------------------------------------
352 * mode = 0. See if any gradients will have to be estimated.
353 * This call is from m5frmc after a refactorize,
354 * or from m6srch at the beginning of a linesearch.
355 * ---------------------------------------------------------------
356 grdcon = nncon .eq. 0 .or. lderiv .ge. 2 .or. nlag .eq. 0
357 grdobj = nnobj .eq. 0 .or. lderiv .eq. 1 .or. lderiv .eq. 3
358 if (.not. grdcon) call m6dmmy( njac , z(lgcon) )
359 if (.not. grdobj) call m6dmmy( nnobj, z(lgobj) )
360 else
361 * ---------------------------------------------------------------
362 * mode = 1. Evaluate the subproblem objective.
363 * Copy the free variables x into xn first.
364 * ---------------------------------------------------------------
365 call m5bsx ( 1, ms, nb, z(lkb), x, xn )
366
367 call m6fun1( modefg, ms, n, nn,
368 $ ne, nka, a, ha, ka,
369 $ nncon, nncon0, nnobj, nnobj0, nnjac, njac,
370 $ z(lfcon), z(lgcon), fobj, z(lgobj), fsub,
371 $ z(lfdif), z(lfold), z(lxlam),
372 $ x, xn, z(lxdif), z(lxold), z, nwcore )
373 end if
374
375 * end of m6fun
376 end
377
378 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
379
380 subroutine m6fun1( modefg, ms, n, nn,
381 $ ne, nka, a, ha, ka,
382 $ nncon, nncon0, nnobj, nnobj0, nnjac, njac,
383 $ fcon, gcon, fobj, gobj, fsub,
384 $ fdif, fold, xlam,
385 $ x, xn, xdif, xold, z, nwcore )
386
387 implicit double precision (a-h,o-z)
388 integer*4 ha(ne)
389 integer ka(nka)
390 double precision a(ne),
391 $ fcon(nncon0), gcon(njac ), gobj(nnobj0),
392 $ fdif(nncon0), fold(nncon0), xlam(nncon0),
393 $ x(ms), xn(nn), xdif(nn), xold(nn), z(nwcore)
394
395 * ------------------------------------------------------------------
396 * m6fun1 calls the constraint and objective functions, and then
397 * computes the subproblem objective function from them.
398 *
399 * 03 Sep 1992: Penalty parameter (penpar) is now used as a
400 * multiple of a default value, pendef = 100 / nncon.
401 * ------------------------------------------------------------------
402
403 common /m1file/ iread,iprint,isumm
404 common /m3scal/ sclobj,scltol,lscale
405 common /m5lobj/ sinf,wtobj,minimz,ninf,iobj,jobj,kobj
406 common /m5log1/ idebug,ierr,lprint
407 common /m8al1 / penpar,rowtol,ncom,nden,nlag,nmajor,nminor
408
409 logical lagrng
410 parameter ( zero = 0.0d+0, dhalf = 0.5d+0 )
411
412 * Evaluate the functions (constraints first, then objective).
413 * If Lagrangian = No, we only evaluate the true objective.
414
415 lagrng = nnjac .gt. 0 .and. nlag .gt. 0
416
417 if (lagrng) then
418 call m6fcon( modefg, nncon, nnjac, njac, fcon, gcon,
419 $ ne, nka, ha, ka,
420 $ xn, z, nwcore )
421 if (ierr .ne. 0) return
422 end if
423
424 fobj = zero
425 if (nnobj .gt. 0) then
426 call m6fobj( modefg, nnobj, fobj, gobj, xn, z, nwcore )
427 if (ierr .ne. 0) return
428 end if
429
430 * Include the linear objective.
431
432 fsub = fobj
433 if (iobj .ne. 0) fsub = fsub - x(iobj) * sclobj
434 if (minimz .lt. 0) fsub = - fsub
435
436 * Compute the augmented Lagrangian terms:
437 * fdif = (fcon - fold) - (Jacobian)*(xn - xold).
438
439 if (lagrng) then
440 do 110 i = 1, nncon
441 fdif(i) = fcon(i) - fold(i)
442 110 continue
443
444 do 120 j = 1, nnjac
445 xdif(j) = xn(j) - xold(j)
446 120 continue
447
448 call m2aprd( 7, xdif, nnjac, fdif, nncon,
449 $ ne, nka, a, ha, ka,
450 $ z, nwcore )
451
452 fsub = fsub - ddot( nncon, xlam, 1, fdif, 1 )
453 if (penpar .gt. zero) then
454 pendef = 100.0d+0 / nncon
455 rho = penpar * pendef
456 fdnorm = dnrm2( nncon, fdif, 1 )
457 fsub = fsub + dhalf * rho * fdnorm**2
458 end if
459 end if
460
461 900 return
462
463 * end of m6fun1
464 end
465
466 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
467
468 subroutine m6grd ( ms, nb, nn, gsub, grd,
469 $ ne, nka, a, ha, ka,
470 $ xn, z, nwcore )
471
472 implicit double precision (a-h,o-z)
473 integer*4 ha(ne)
474 integer ka(nka)
475 double precision a(ne), gsub(nn), grd(ms), xn(nb), z(nwcore)
476
477 * ------------------------------------------------------------------
478 * m6grd returns the gradient of the subproblem objective in gsub.
479 *
480 * 14 Oct 1991: a, ha, ka added as parameters.
481 * ------------------------------------------------------------------
482
483 common /m3loc / lascal,lbl ,lbu ,lbbl ,lbbu ,
484 $ lhrtyp,lhs ,lkb
485 common /m5loc / lpi ,lpi2 ,lw ,lw2 ,
486 $ lx ,lx2 ,ly ,ly2 ,
487 $ lgsub ,lgsub2,lgrd ,lgrd2 ,
488 $ lr ,lrg ,lrg2 ,lxn
489 common /m7len / fobj ,fobj2 ,nnobj ,nnobj0
490 common /m7loc / lgobj ,lgobj2
491 common /m8len / njac ,nncon ,nncon0,nnjac
492 common /m8loc / lfcon ,lfcon2,lfdif ,lfdif2,lfold ,
493 $ lblslk,lbuslk,lxlam ,lrhs ,
494 $ lgcon ,lgcon2,lxdif ,lxold
495
496 call m6grd1( n, nn, fobj,
497 $ nncon, nncon0, nnobj, nnobj0, nnjac, njac,
498 $ ne, nka, a, ha, ka,
499 $ z(lfcon), z(lfcon2), z(lfdif ), z(lxlam),
500 $ z(lgcon), z(lgcon2), z(lgobj ), z(lgobj2), gsub,
501 $ xn, z(ly), z, nwcore )
502
503 call m7bsg ( ms, nn, z(lkb), gsub, grd )
504
505 * end of m6grd
506 end
507
508 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
509
510 subroutine m6grd1( n, nn, fobj,
511 $ nncon, nncon0, nnobj, nnobj0, nnjac, njac,
512 $ ne, nka, a, ha, ka,
513 $ fcon, fcon2, fdif, xlam,
514 $ gcon, gcon2, gobj, gobj2, gsub,
515 $ xn, y, z, nwcore )
516
517 implicit double precision (a-h,o-z)
518 integer*4 ha(ne)
519 integer ka(nka)
520 double precision a(ne), fcon(nncon0), fcon2(nncon0),
521 $ fdif(nncon0), xlam (nncon0),
522 $ gcon(njac ), gcon2(njac ),
523 $ gobj(nnobj0), gobj2(nnobj0), gsub(nn),
524 $ xn(nn), y(nncon0), z(nwcore)
525
526 * ------------------------------------------------------------------
527 * m6grd1 computes gsub(j), the gradient of the subproblem objective.
528 *
529 * NOTE -- m6dcon overwrites the first nncon elements of y
530 * if central differences are needed.
531 *
532 * 03 Sep 1992: Penalty parameter (penpar) is now used as a
533 * multiple of a default value, pendef = 100 / nncon.
534 * ------------------------------------------------------------------
535
536 common /m5lobj/ sinf,wtobj,minimz,ninf,iobj,jobj,kobj
537 common /m5log1/ idebug,ierr,lprint
538 common /m8al1 / penpar,rowtol,ncom,nden,nlag,nmajor,nminor
539 common /m8diff/ difint(2),gdummy,lderiv,lvldif,knowng(2)
540
541 logical grdcon, grdobj, lagrng
542 parameter ( zero = 0.0d+0, one = 1.0d+0 )
543
544 grdcon = nncon .eq. 0 .or. lderiv .ge. 2 .or. nlag .eq. 0
545 grdobj = nnobj .eq. 0 .or. lderiv .eq. 1 .or. lderiv .eq. 3
546 lagrng = nnjac .gt. 0 .and. nlag .gt. 0
547
548 * ------------------------------------------------------
549 * Fill in any missing constraint or objective gradients.
550 * ------------------------------------------------------
551
552 if (.not. grdcon) then
553 call m6dcon( nncon, nnjac, njac,
554 $ ne, nka, ha, ka,
555 $ fcon, fcon2, gcon, gcon2, xn, y, z, nwcore )
556 if (ierr .ne. 0) return
557 end if
558
559 if (.not. grdobj) then
560 call m6dobj( nnobj, fobj, gobj, gobj2, xn, z, nwcore )
561 if (ierr .ne. 0) return
562 end if
563
564 * -----------------------------------------
565 * Compute the required elements of gsub.
566 * -----------------------------------------
567
568 * Set gsub = ( gobj, 0 ) and allow for the sign of the objective.
569
570 l = nn - nnobj
571 if (nnobj .gt. 0) call dcopy ( nnobj , gobj, 1, gsub, 1 )
572 if (l .gt. 0) call dload ( l , zero, gsub(nnobj+1), 1 )
573 if (minimz .lt. 0) call dscal ( nnobj0, (- one), gsub, 1 )
574
575 * Compute the augmented Lagrangian terms.
576 * m6fun1 has set up fdif = (fcon - fold) - Jac * (x - xold).
577 * To save work, set fdif = rho * fdif - xlam.
578
579 if (lagrng) then
580 pendef = 100.0d+0 / nncon
581 rho = penpar * pendef
582 call dscal ( nncon, rho, fdif, 1 )
583 call daxpy ( nncon, (- one), xlam, 1, fdif, 1 )
584
585 l = 0
586 do 600 j = 1, nnjac
587 g = gsub(j)
588 k1 = ka(j)
589 k2 = ka(j+1) - 1
590
591 do 500 k = k1, k2
592 ir = ha(k)
593 if (ir .gt. nncon) go to 590
594 l = l + 1
595 g = g + (gcon(l) - a(k)) * fdif(ir)
596 500 continue
597
598 590 gsub(j) = g
599 600 continue
600 end if
601
602 * end of m6grd1
603 end
604
605 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
606
607 subroutine m6dobj( nnobj, fobj, gobj, gobj2, xn, z, nwcore )
608
609 implicit double precision (a-h,o-z)
610 double precision gobj(nnobj), gobj2(nnobj)
611 double precision xn(nnobj), z(nwcore)
612
613 * ------------------------------------------------------------------
614 * m6dobj estimates missing objective gradients
615 * using finite differences of the objective function fobj.
616 * ------------------------------------------------------------------
617
618 common /m5log1/ idebug,ierr,lprint
619 common /m8diff/ difint(2),gdummy,lderiv,lvldif,knowng(2)
620 common /m8func/ nfcon(4),nfobj(4),nprob,nstat1,nstat2
621
622 intrinsic abs
623
624 logical centrl
625 parameter ( one = 1.0d+0 )
626
627 centrl = lvldif .eq. 2
628 delta = difint(lvldif)
629 numf = 0
630 fback = fobj
631
632 do 200 j = 1, nnobj
633 if (gobj(j) .eq. gdummy) then
634 xj = xn(j)
635 dx = delta*(one + abs( xj ))
636 xn(j) = xj + dx
637 numf = numf + 1
638 call m6fobj( 0, nnobj, fforwd, gobj2, xn, z, nwcore )
639 if (ierr .ne. 0) go to 900
640
641 if (centrl) then
642 xn(j) = xj - dx
643 dx = dx + dx
644 numf = numf + 1
645 call m6fobj( 0, nnobj, fback, gobj2, xn, z, nwcore )
646 if (ierr .ne. 0) go to 900
647 end if
648
649 gobj(j) = (fforwd - fback) / dx
650 xn(j) = xj
651 end if
652 200 continue
653
654 900 l = lvldif + 2
655 nfobj(l) = nfobj(l) + numf
656
657 * end of m6dobj
658 end
659
660 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
661
662 subroutine m6dcon( nncon, nnjac, njac,
663 $ ne, nka, ha, ka,
664 $ fcon, fcon2, gcon, gcon2, xn, y, z, nwcore )
665
666 implicit double precision (a-h,o-z)
667 integer*4 ha(ne)
668 integer ka(nka)
669 double precision fcon(nncon), fcon2(nncon),
670 $ gcon(njac ), gcon2(njac ),
671 $ xn(nnjac), y(nncon), z(nwcore)
672
673 * ------------------------------------------------------------------
674 * m6dcon estimates missing elements in the columns of the
675 * Jacobian (the first nncon rows and nnjac cols of A)
676 * using finite differences of the constraint functions fcon.
677 * ------------------------------------------------------------------
678
679 common /m5log1/ idebug,ierr,lprint
680 common /m8diff/ difint(2),gdummy,lderiv,lvldif,knowng(2)
681 common /m8func/ nfcon(4),nfobj(4),nprob,nstat1,nstat2
682
683 intrinsic abs
684
685 logical centrl
686 parameter ( one = 1.0d+0 )
687
688 centrl = lvldif .eq. 2
689 delta = difint(lvldif)
690 numf = 0
691 l = 0
692
693 * If central differences are needed, we first save fcon in y.
694 * When fcon is later loaded with f(x - h) for each pertbn h,
695 * it can be used the same as when forward differences are in effect.
696
697 if (centrl) call dcopy ( nncon, fcon, 1, y, 1 )
698
699 do 200 j = 1, nnjac
700
701 * Look for the first missing element in this column.
702
703 lsave = l
704 k1 = ka(j)
705 k2 = ka(j+1) - 1
706 do 100 k = k1, k2
707 ir = ha(k)
708 if (ir .gt. nncon) go to 200
709 l = l + 1
710 if (gcon(l) .eq. gdummy) go to 120
711 100 continue
712 go to 200
713
714 * Found one. Approximate it and any others by finite differences.
715
716 120 l = lsave
717 xj = xn(j)
718 dx = delta * (one + abs( xj ))
719 xn(j) = xj + dx
720 numf = numf + 1
721 call m6fcon( 0, nncon, nnjac, njac, fcon2, gcon2,
722 $ ne, nka, ha, ka,
723 $ xn, z, nwcore )
724 if (ierr .ne. 0) go to 900
725
726 if (centrl) then
727 xn(j) = xj - dx
728 dx = dx + dx
729 numf = numf + 1
730 call m6fcon( 0, nncon, nnjac, njac, fcon, gcon2,
731 $ ne, nka, ha, ka,
732 $ xn, z, nwcore )
733 if (ierr .ne. 0) go to 900
734 end if
735
736 do 150 k = k1, k2
737 ir = ha(k)
738 if (ir .gt. nncon) go to 190
739 l = l + 1
740 if (gcon(l) .eq. gdummy) then
741 gcon(l) = (fcon2(ir) - fcon(ir)) / dx
742 end if
743 150 continue
744 190 xn(j) = xj
745 200 continue
746
747 900 if (centrl) call dcopy ( nncon, y, 1, fcon, 1 )
748 l = lvldif + 2
749 nfcon(l) = nfcon(l) + numf
750
751 * end of m6dcon
752 end
753
754 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
755
756 subroutine m6srch( ms, ns, n, nb, nn, itn, inform,
757 $ debug, fonly, switch,
758 $ ne, nka, a, ha, ka,
759 $ alfa, alfmax, difint, epsmch, epsrf, eta,
760 $ fsub, gnorm , pnorm , xnorm,
761 $ gsub, grd, p, x, x2, xn, z, nwcore )
762
763 implicit double precision (a-h,o-z)
764 logical debug, fonly, switch
765 integer*4 ha(ne)
766 integer ka(nka)
767 double precision a(ne), gsub(nn), grd(ms),
768 $ p(ms), x(ms), x2(ms), xn(nb), z(nwcore)
769
770 * ------------------------------------------------------------------
771 * m6srch finds a step-length alfa along the search direction p,
772 * such that the (subproblem) objective function F is sufficiently
773 * reduced: F(x + alpha*p) < F(x).
774 *
775 * Input parameters...
776 *
777 * alfa Initial estimate of the step.
778 * alfmax Maximum allowable step.
779 * difint A difference interval that has been used for forward
780 * differences (not central differences) to get p.
781 * (difint is used if switch is true).
782 * epsmch Relative machine precision.
783 * epsrf Relative function precision.
784 * eta Linesearch accuracy parameter in the range (0,1).
785 * 0.001 means very accurate. 0.99 means quite sloppy.
786 * fonly true if function-only search should be used.
787 * switch true if there is a possibility of switching to
788 * central differences to get a better p.
789 * fsub Current value of subproblem objective, F(x).
790 * gsub(*) Current gradient of subproblem objective.
791 * grd(*) gsub reordered like x and p (free variables only).
792 * p(*) Search direction for the free variables.
793 * x(*) Current free variables (basic and superbasic).
794 *
795 * Output parameters...
796 *
797 * alfa
798 * fsub New value of the subproblem objective, F(x + alfa*p).
799 * gsub(*) New gradient of the subproblem objective.
800 * grd(*) New gsub reordered parallel to x and p.
801 * x New free variables (basic and superbasic).
802 * inform = -1 if the user wants to stop.
803 * = 1 to 8 -- see below.
804 *
805 * 02 Oct 1991: fonly search uses modefg = 0 during the search,
806 * and modefg = 2 after the search.
807 * 04 Oct 1991: switch added to suppress request for central diffs
808 * when fonly = true but all gradients are really known.
809 * 16 Oct 1991: getptc, getptq replaced by srchc, srchq.
810 * inform values altered.
811 * ------------------------------------------------------------------
812
813 common /m1file/ iread,iprint,isumm
814 common /m5log1/ idebug,ierr,lprint
815
816 intrinsic abs , max
817 logical done , first , imprvd
818 parameter ( zero = 0.0d+0, one = 1.0d+0, two = 2.0d+0 )
819 parameter ( rmu = 1.0d-4 )
820
821 * ------------------------------------------------------------------
822 * Set the input parameters for srchc or srchq.
823 *
824 * alfsml is used by srchq. If alfa would be less than alfsml,
825 * the search will be terminated early.
826 * If p was computed using forward or backward differences,
827 * alfsml should be positive (and related to the difference
828 * interval used).
829 * If p was computed using central differences (lvldif = 2)
830 * alfsml should be zero.
831 *
832 * epsrf is the relative function precision (specified by user).
833 *
834 * epsaf is the absolute function precision. If F(x1) and F(x2) are
835 * as close as epsaf, we cannot safely conclude from the
836 * function values alone which of x1 or x2 is a better point.
837 *
838 * tolabs is an estimate of the absolute spacing between points
839 * along p. This step should produce a perturbation
840 * of epsaf in the subproblem objective function.
841 *
842 * tolrel is an estimate of the relative spacing between points
843 * along p.
844 *
845 * toltny is the minimum allowable absolute spacing between points
846 * along p.
847 * ------------------------------------------------------------------
848 first = .true.
849 maxf = 10
850 if (fonly) maxf = 15
851 nout = iprint
852
853 alfsml = zero
854 if (fonly .and. switch) then
855 alfsml = difint * (one + xnorm) / pnorm
856 end if
857
858 epsaf = max( epsrf, epsmch ) * (one + abs( fsub ))
859 oldf = fsub
860 oldg = ddot ( ms, grd, 1, p, 1 )
861 eps0 = epsmch**0.8
862 tolax = eps0
863 tolrx = eps0
864 tolrel = max( tolrx, epsmch )
865 t = xnorm * tolrx + tolax
866 tolabs = alfmax
867 if (t .lt. pnorm * alfmax) tolabs = t / pnorm
868
869 t = zero
870 do 10 j = 1, ms
871 s = abs( p(j) )
872 q = abs( x(j) )*tolrx + tolax
873 if (s .gt. q*t) t = s / q
874 10 continue
875
876 if (t*tolabs .gt. one) then
877 toltny = one / t
878 else
879 toltny = tolabs
880 end if
881
882 alfbst = zero
883 fbest = zero
884 gbest = (one - rmu)*oldg
885 targtg = (rmu - eta)*oldg
886 g0 = gbest
887
888 if (debug) write(nout, 1000) itn, pnorm, gnorm
889
890 * ------------------------------------------------------------------
891 * Commence main loop, entering srchc or srchq two or more times.
892 * first = true for the first entry, false for subsequent entries.
893 * done = true indicates termination, in which case the value of
894 * inform gives the result of the search.
895 * inform = 1 if the search is successful and alfa < alfmax.
896 * = 2 if the search is successful and alfa = alfmax.
897 * = 3 if a better point was found but too many functions
898 * were needed (not sufficient decrease).
899 * = 4 if alfmax < tolabs (too small to do a search).
900 * = 5 if alfa < alfsml (srchq only -- maybe want to switch
901 * to central differences to get a better direction).
902 * = 6 if the search found that there is no useful step.
903 * The interval of uncertainty is less than 2*tolabs.
904 * The minimizer is very close to alfa = zero
905 * or the gradients are not sufficiently accurate.
906 * = 7 if there were too many function calls.
907 * = 8 if the input parameters were bad
908 * (alfmax le toltny or oldg ge 0).
909 * ------------------------------------------------------------------
910 * Start loop
911 100 if (fonly) then
912 call srchq ( first , debug , done , imprvd, inform,
913 $ maxf , numf , nout ,
914 $ alfmax, alfsml, epsaf ,
915 $ g0 , targtg, ftry ,
916 $ tolabs, tolrel, toltny,
917 $ alfa , alfbst, fbest )
918 else
919 call srchc ( first , debug , done , imprvd, inform,
920 $ maxf , numf , nout ,
921 $ alfmax, epsaf ,
922 $ g0 , targtg, ftry , gtry ,
923 $ tolabs, tolrel, toltny,
924 $ alfa , alfbst, fbest , gbest )
925 end if
926
927 if (.not. done ) then
928 * ------------------------------------------------------------
929 * Compute the function value fsub at the new trial point
930 * x + alfa*p.
931 * ------------------------------------------------------------
932 do 180 j = 1, ms
933 x2(j) = x(j) + alfa*p(j)
934 180 continue
935
936 * Check for first function request.
937 * Maybe set dummy gradients.
938
939 if (numf .eq. 0) then
940 modefg = 2
941 if ( fonly ) then
942 modefg = 0
943 call m6fun ( 0, modefg, n, nb, ms, fsub,
944 $ ne, nka, a, ha, ka,
945 $ x2, xn, z, nwcore )
946 end if
947 end if
948
949 call m6fun ( 1, modefg, n, nb, ms, fsub,
950 $ ne, nka, a, ha, ka,
951 $ x2, xn, z, nwcore )
952 if (ierr .ne. 0) go to 999
953
954 ftry = fsub - oldf - rmu*oldg*alfa
955
956 * If gradients are being used, also get the
957 * directional derivative gtry.
958
959 if (.not. fonly) then
960 call m6grd ( ms, nb, nn, gsub, grd,
961 $ ne, nka, a, ha, ka,
962 $ xn, z, nwcore )
963 gtp = ddot( ms, grd, 1, p, 1 )
964 gtry = gtp - rmu*oldg
965 end if
966
967 go to 100
968 end if
969 * end loop
970
971 * ------------------------------------------------------------------
972 * The search has finished.
973 * If some function evaluations were made,
974 * set alfa = alfbst and x = x + alfa*p.
975 * ------------------------------------------------------------------
976 if (inform .le. 7) then
977 alfa = alfbst
978 call daxpy ( ms, alfa, p, 1, x, 1 )
979
980 * If only function values were used in the search,
981 * or if the last call to m6fun was not at the best point,
982 * we have to recompute the function and gradient there
983 * (and reset grd).
984 * Beware: if central differences are used, m6grd overwrites p.
985
986 if (fonly .or. .not. imprvd) then
987 modefg = 2
988 call m6fun ( 1, modefg, n, nb, ms, fsub,
989 $ ne, nka, a, ha, ka,
990 $ x, xn, z, nwcore )
991 if (ierr .ne. 0) go to 999
992
993 call m6grd ( ms, nb, nn, gsub, grd,
994 $ ne, nka, a, ha, ka,
995 $ xn, z, nwcore )
996 if (ierr .ne. 0) go to 999
997 end if
998 end if
999
1000 * ------------------------------------------------------------------
1001 * Print a few error msgs if necessary.
1002 * ------------------------------------------------------------------
1003 if (inform .ge. 7) then
1004 if (inform .eq. 7) then
1005 *- if (iprint .gt. 0) write(iprint, 1700) numf
1006 *- if (isumm .gt. 0) write(isumm , 1700) numf
1007 else if (oldg .ge. zero) then
1008 *- if (iprint .gt. 0) write(iprint, 1600) oldg
1009 end if
1010
1011 if (.not. debug) then
1012 if (iprint .gt. 0) then
1013 write(iprint, 1800) alfmax, pnorm, gnorm, oldg, numf
1014 end if
1015 end if
1016 end if
1017
1018 return
1019
1020 * The user wants to stop.
1021
1022 999 inform = -1
1023 return
1024
1025 1000 format(// ' --------------------------------------------'
1026 $ / ' Output from m6srch following iteration', i7,
1027 $ 5x, ' Norm p =', 1p, e11.2, 5x, ' Norm g =', e11.2)
1028 *1600 format(/ ' XXX The search direction is uphill. gtp =', 1p,e9.1)
1029 *1700 format( ' XXX The linesearch has evaluated the functions', i5,
1030 * $ ' times')
1031 1800 format(' alfmax =', 1p, e11.2, ' pnorm =', e11.2,
1032 $ ' gnorm =', e11.2, ' g(t)p =', e11.2,
1033 $ ' numf =', i3)
1034
1035 * end of m6srch
1036 end
1037
1038 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1039
1040 subroutine srchc ( first , debug , done , imprvd, inform,
1041 $ maxf , numf , nout ,
1042 $ alfmax, epsaf ,
1043 $ g0 , targtg, ftry , gtry ,
1044 $ tolabs, tolrel, toltny,
1045 $ alfa , alfbst, fbest , gbest )
1046
1047 implicit double precision (a-h,o-z)
1048 logical first , debug , done , imprvd
1049
1050 ************************************************************************
1051 * srchc finds a sequence of improving estimates of a minimizer of
1052 * the univariate function f(alpha) in the interval (0,alfmax].
1053 * f(alpha) is a smooth function such that f(0) = 0 and f'(0) < 0.
1054 * srchc requires both f(alpha) and f'(alpha) to be evaluated at
1055 * points in the interval. Estimates of the minimizer are computed
1056 * using safeguarded cubic interpolation.
1057 *
1058 * Reverse communication is used to allow the calling program to
1059 * evaluate f and f'. Some of the parameters must be set or tested
1060 * by the calling program. The remainder would ordinarily be local
1061 * variables.
1062 *
1063 * Input parameters (relevant to the calling program)
1064 * --------------------------------------------------
1065 *
1066 * first must be true on the first entry. It is subsequently
1067 * altered by srchc.
1068 *
1069 * debug specifies whether detailed output is wanted.
1070 *
1071 * maxf is an upper limit on the number of times srchc is to
1072 * be entered consecutively with done = false
1073 * (following an initial entry with first = true).
1074 *
1075 * alfa is the first estimate of a minimizer. alfa is
1076 * subsequently altered by srchc (see below).
1077 *
1078 * alfmax is the upper limit of the interval to be searched.
1079 *
1080 * epsaf is an estimate of the absolute precision in the
1081 * computed value of f(0).
1082 *
1083 * ftry, gtry are the values of f, f' at the new point
1084 * alfa = alfbst + xtry.
1085 *
1086 * g0 is the value of f'(0). g0 must be negative.
1087 *
1088 * tolabs,tolrel define a function tol(alfa) = tolrel*alfa + tolabs
1089 * such that if f has already been evaluated at alfa,
1090 * it will not be evaluated closer than tol(alfa).
1091 * These values may be reduced by srchc.
1092 *
1093 * targtg is the target value of abs(f'(alfa)). The search
1094 * is terminated when
1095 * abs(f'(alfa)) le targtg and f(alfa) lt 0.
1096 *
1097 * toltny is the smallest value that tolabs is allowed to be
1098 * reduced to.
1099 *
1100 * Output parameters (relevant to the calling program)
1101 * ---------------------------------------------------
1102 *
1103 * imprvd is true if the previous alfa was the best point so
1104 * far. Any related quantities should be saved by the
1105 * calling program (e.g., gradient arrays) before
1106 * paying attention to the variable done.
1107 *
1108 * done = false means the calling program should evaluate
1109 * ftry = f(alfa), gtry = f'(alfa)
1110 * for the new trial alfa, and re-enter srchc.
1111 *
1112 * done = true means that no new alfa was calculated. The value
1113 * of inform gives the result of the search as follows
1114 *
1115 * inform = 1 means the search has terminated
1116 * successfully with alfbst < alfmax.
1117 *
1118 * inform = 2 means the search has terminated
1119 * successfully with alfbst = alfmax.
1120 *
1121 * inform = 3 means that the search failed to find a
1122 * point of sufficient decrease in maxf
1123 * functions, but a lower point was found.
1124 *
1125 * inform = 4 means alfmax is so small that a search
1126 * should not have been attempted.
1127 *
1128 * inform = 5 is never set by srchc.
1129 *
1130 * inform = 6 means the search has failed to find a
1131 * useful step. The interval of uncertainty
1132 * is [0,b] with b < 2*tolabs. A minimizer
1133 * lies very close to alfa = 0, or f'(0) is
1134 * not sufficiently accurate.
1135 *
1136 * inform = 7 if no better point could be found after
1137 * maxf function calls.
1138 *
1139 * inform = 8 means the input parameters were bad.
1140 * alfmax le toltny or g0 ge zero.
1141 * No function evaluations were made.
1142 *
1143 * numf counts the number of times srchc has been entered
1144 * consecutively with done = false (i.e., with a new
1145 * function value ftry).
1146 *
1147 * alfa is the point at which the next function ftry and
1148 * derivative gtry must be computed.
1149 *
1150 * alfbst should be accepted by the calling program as the
1151 * approximate minimizer, whenever srchc returns
1152 * inform = 1 or 2 (and possibly 3).
1153 *
1154 * fbest, gbest will be the corresponding values of f, f'.
1155 *
1156 *
1157 * The following parameters retain information between entries
1158 * -----------------------------------------------------------
1159 *
1160 * braktd is false if f and f' have not been evaluated at
1161 * the far end of the interval of uncertainty. In this
1162 * case, the point b will be at alfmax + tol(alfmax).
1163 *
1164 * crampd is true if alfmax is very small (le tolabs). If the
1165 * search fails, this indicates that a zero step should
1166 * be taken.
1167 *
1168 * extrap is true if xw lies outside the interval of
1169 * uncertainty. In this case, extra safeguards are
1170 * applied to allow for instability in the polynomial
1171 * fit.
1172 *
1173 * moved is true if a better point has been found, i.e.,
1174 * alfbst gt 0.
1175 *
1176 * wset records whether a second-best point has been
1177 * determined it will always be true when convergence
1178 * is tested.
1179 *
1180 * nsamea is the number of consecutive times that the
1181 * left-hand end point of the interval of uncertainty
1182 * has remained the same.
1183 *
1184 * nsameb similarly for the right-hand end.
1185 *
1186 * a, b, alfbst define the current interval of uncertainty.
1187 * A minimizer lies somewhere in the interval
1188 * [alfbst + a, alfbst + b].
1189 *
1190 * alfbst is the best point so far. It is always at one end
1191 * of the interval of uncertainty. hence we have
1192 * either a lt 0, b = 0 or a = 0, b gt 0.
1193 *
1194 * fbest, gbest are the values of f, f' at the point alfbst.
1195 *
1196 * factor controls the rate at which extrapolated estimates
1197 * of alfa may expand into the interval of uncertainty.
1198 * factor is not used if a minimizer has been bracketed
1199 * (i.e., when the variable braktd is true).
1200 *
1201 * fw, gw are the values of f, f' at the point alfbst + xw.
1202 * they are not defined until wset is true.
1203 *
1204 * xtry is the trial point in the shifted interval (a, b).
1205 *
1206 * xw is such that alfbst + xw is the second-best point.
1207 * it is not defined until wset is true.
1208 * in some cases, xw will replace a previous xw
1209 * that has a lower function but has just been excluded
1210 * from the interval of uncertainty.
1211 *
1212 *
1213 * Systems Optimization Laboratory, Stanford University, California.
1214 * Original version February 1982. Rev. May 1983.
1215 * Original f77 version 22-August-1985.
1216 * This version of srchc dated 24-Oct-91.
1217 ************************************************************************
1218
1219 logical braktd, crampd, extrap, moved , wset
1220 save braktd, crampd, extrap, moved , wset
1221
1222 save nsamea, nsameb
1223 save a , b , factor
1224 save xtry , xw , fw , gw , tolmax
1225
1226 logical badfun, closef, found
1227 logical quitF , quitI
1228 logical fitok , setxw
1229 intrinsic abs , sqrt
1230
1231 parameter ( zero =0.0d+0, point1 =0.1d+0, half =0.5d+0 )
1232 parameter ( one =1.0d+0, two =2.0d+0, three =3.0d+0 )
1233 parameter ( five =5.0d+0, ten =1.0d+1, eleven =1.1d+1 )
1234
1235 * ------------------------------------------------------------------
1236 * Local variables
1237 * ===============
1238 *
1239 * closef is true if the new function ftry is within epsaf of
1240 * fbest (up or down).
1241 *
1242 * found is true if the sufficient decrease conditions hold at
1243 * alfbst.
1244 *
1245 * quitF is true when maxf function calls have been made.
1246 *
1247 * quitI is true when the interval of uncertainty is less than
1248 * 2*tol.
1249 * ---------------------------------------------------------------------
1250
1251 badfun = .false.
1252 quitF = .false.
1253 quitI = .false.
1254 imprvd = .false.
1255
1256 if (first) then
1257 * ---------------------------------------------------------------
1258 * First entry. Initialize various quantities, check input data
1259 * and prepare to evaluate the function at the initial alfa.
1260 * ---------------------------------------------------------------
1261 first = .false.
1262 numf = 0
1263 alfbst = zero
1264 badfun = alfmax .le. toltny .or. g0 .ge. zero
1265 done = badfun
1266 moved = .false.
1267
1268 if (.not. done) then
1269 braktd = .false.
1270 crampd = alfmax .le. tolabs
1271 extrap = .false.
1272 wset = .false.
1273 nsamea = 0
1274 nsameb = 0
1275
1276 tolmax = tolabs + tolrel*alfmax
1277 a = zero
1278 b = alfmax + tolmax
1279 factor = five
1280 tol = tolabs
1281 xtry = alfa
1282 if (debug) then
1283 write (nout, 1000) g0 , tolabs, alfmax,
1284 $ targtg, tolrel, epsaf , crampd
1285 end if
1286 end if
1287 else
1288 * ---------------------------------------------------------------
1289 * Subsequent entries. The function has just been evaluated at
1290 * alfa = alfbst + xtry, giving ftry and gtry.
1291 * ---------------------------------------------------------------
1292 if (debug) write (nout, 1100) alfa, ftry, gtry
1293
1294 numf = numf + 1
1295 nsamea = nsamea + 1
1296 nsameb = nsameb + 1
1297
1298 if (.not. braktd) then
1299 tolmax = tolabs + tolrel*alfmax
1300 b = alfmax - alfbst + tolmax
1301 end if
1302
1303 * See if the new step is better. If alfa is large enough that
1304 * ftry can be distinguished numerically from zero, the function
1305 * is required to be sufficiently negative.
1306
1307 closef = abs( ftry - fbest ) .le. epsaf
1308 if (closef) then
1309 imprvd = abs( gtry ) .le. abs( gbest )
1310 else
1311 imprvd = ftry .lt. fbest
1312 end if
1313
1314 if (imprvd) then
1315
1316 * We seem to have an improvement. The new point becomes the
1317 * origin and other points are shifted accordingly.
1318
1319 fw = fbest
1320 fbest = ftry
1321 gw = gbest
1322 gbest = gtry
1323 alfbst = alfa
1324 moved = .true.
1325
1326 a = a - xtry
1327 b = b - xtry
1328 xw = zero - xtry
1329 wset = .true.
1330 extrap = xw .lt. zero .and. gbest .lt. zero
1331 $ .or. xw .gt. zero .and. gbest .gt. zero
1332
1333 * Decrease the length of the interval of uncertainty.
1334
1335 if (gtry .le. zero) then
1336 a = zero
1337 nsamea = 0
1338 else
1339 b = zero
1340 nsameb = 0
1341 braktd = .true.
1342 end if
1343 else
1344
1345 * The new function value is not better than the best point so
1346 * far. The origin remains unchanged but the new point may
1347 * qualify as xw. xtry must be a new bound on the best point.
1348
1349 if (xtry .le. zero) then
1350 a = xtry
1351 nsamea = 0
1352 else
1353 b = xtry
1354 nsameb = 0
1355 braktd = .true.
1356 end if
1357
1358 * If xw has not been set or ftry is better than fw, update the
1359 * points accordingly.
1360
1361 if (wset) then
1362 setxw = ftry .lt. fw .or. .not. extrap
1363 else
1364 setxw = .true.
1365 end if
1366
1367 if (setxw) then
1368 xw = xtry
1369 fw = ftry
1370 gw = gtry
1371 wset = .true.
1372 extrap = .false.
1373 end if
1374 end if
1375
1376 * ---------------------------------------------------------------
1377 * Check the termination criteria. wset will always be true.
1378 * ---------------------------------------------------------------
1379 tol = tolabs + tolrel*alfbst
1380 truea = alfbst + a
1381 trueb = alfbst + b
1382
1383 found = abs(gbest) .le. targtg
1384 quitF = numf .ge. maxf
1385 quitI = b - a .le. tol + tol
1386
1387 if (quitI .and. .not. moved) then
1388
1389 * The interval of uncertainty appears to be small enough,
1390 * but no better point has been found. Check that changing
1391 * alfa by b-a changes f by less than epsaf.
1392
1393 tol = tol/ten
1394 tolabs = tol
1395 quitI = abs(fw) .le. epsaf .or. tol .le. toltny
1396 end if
1397
1398 done = quitF .or. quitI .or. found
1399
1400 if (debug) then
1401 write (nout, 1200) truea , trueb , b - a , tol ,
1402 $ nsamea , nsameb, numf ,
1403 $ braktd , extrap, closef, imprvd,
1404 $ found , quitI ,
1405 $ alfbst , fbest , gbest ,
1406 $ alfbst+xw, fw , gw
1407 end if
1408
1409 * ---------------------------------------------------------------
1410 * Proceed with the computation of a trial steplength.
1411 * The choices are...
1412 * 1. Parabolic fit using derivatives only, if the f values are
1413 * close.
1414 * 2. Cubic fit for a minimizer, using both f and f'.
1415 * 3. Damped cubic or parabolic fit if the regular fit appears to
1416 * be consistently overestimating the distance to a minimizer.
1417 * 4. Bisection, geometric bisection, or a step of tol if
1418 * choices 2 or 3 are unsatisfactory.
1419 * ---------------------------------------------------------------
1420 if (.not. done) then
1421 xmidpt = half*(a + b)
1422 s = zero
1423 q = zero
1424
1425 if (closef) then
1426 * ---------------------------------------------------------
1427 * Fit a parabola to the two best gradient values.
1428 * ---------------------------------------------------------
1429 s = gbest
1430 q = gbest - gw
1431 if (debug) write (nout, 2200)
1432 else
1433 * ---------------------------------------------------------
1434 * Fit cubic through fbest and fw.
1435 * ---------------------------------------------------------
1436 if (debug) write (nout, 2100)
1437 fitok = .true.
1438 r = three*(fbest - fw)/xw + gbest + gw
1439 absr = abs( r )
1440 s = sqrt( abs( gbest ) ) * sqrt( abs( gw ) )
1441
1442 * Compute q = the square root of r*r - gbest*gw.
1443 * The method avoids unnecessary underflow and overflow.
1444
1445 if ((gw .lt. zero .and. gbest .gt. zero) .or.
1446 $ (gw .gt. zero .and. gbest .lt. zero)) then
1447 scale = absr + s
1448 if (scale .eq. zero) then
1449 q = zero
1450 else
1451 q = scale*sqrt( (absr/scale)**2 + (s/scale)**2 )
1452 end if
1453 else if (absr .ge. s) then
1454 q = sqrt(absr + s)*sqrt(absr - s)
1455 else
1456 fitok = .false.
1457 end if
1458
1459 if (fitok) then
1460
1461 * Compute a minimizer of the fitted cubic.
1462
1463 if (xw .lt. zero) q = - q
1464 s = gbest - r - q
1465 q = gbest - gw - q - q
1466 end if
1467 end if
1468 * ------------------------------------------------------------
1469 * Construct an artificial interval (artifa, artifb) in which
1470 * the new estimate of a minimizer must lie. Set a default
1471 * value of xtry that will be used if the polynomial fit fails.
1472 * ------------------------------------------------------------
1473 artifa = a
1474 artifb = b
1475 if (.not. braktd) then
1476
1477 * A minimizer has not been bracketed. Set an artificial
1478 * upper bound by expanding the interval xw by a suitable
1479 * factor.
1480
1481 xtry = - factor*xw
1482 artifb = xtry
1483 if (alfbst + xtry .lt. alfmax) factor = five*factor
1484
1485 else if (extrap) then
1486
1487 * The points are configured for an extrapolation.
1488 * Set a default value of xtry in the interval (a, b)
1489 * that will be used if the polynomial fit is rejected. In
1490 * the following, dtry and daux denote the lengths of
1491 * the intervals (a, b) and (0, xw) (or (xw, 0), if
1492 * appropriate). The value of xtry is the point at which
1493 * the exponents of dtry and daux are approximately
1494 * bisected.
1495
1496 daux = abs( xw )
1497 dtry = b - a
1498 if (daux .ge. dtry) then
1499 xtry = five*dtry*(point1 + dtry/daux)/eleven
1500 else
1501 xtry = half * sqrt( daux ) * sqrt( dtry )
1502 end if
1503 if (xw .gt. zero) xtry = - xtry
1504 if (debug) write (nout, 2400) xtry, daux, dtry
1505
1506 * Reset the artificial bounds. If the point computed by
1507 * extrapolation is rejected, xtry will remain at the
1508 * relevant artificial bound.
1509
1510 if (xtry .le. zero) artifa = xtry
1511 if (xtry .gt. zero) artifb = xtry
1512 else
1513
1514 * The points are configured for an interpolation. The
1515 * default value xtry bisects the interval of uncertainty.
1516 * the artificial interval is just (a, b).
1517
1518 xtry = xmidpt
1519 if (debug) write (nout, 2300) xtry
1520 if (nsamea .ge. 3 .or. nsameb .ge. 3) then
1521
1522 * If the interpolation appears to be overestimating the
1523 * distance to a minimizer, damp the interpolation.
1524
1525 factor = factor / five
1526 s = factor * s
1527 else
1528 factor = one
1529 end if
1530 end if
1531 * ------------------------------------------------------------
1532 * The polynomial fits give (s/q)*xw as the new step.
1533 * Reject this step if it lies outside (artifa, artifb).
1534 * ------------------------------------------------------------
1535 if (q .ne. zero) then
1536 if (q .lt. zero) s = - s
1537 if (q .lt. zero) q = - q
1538 if (s*xw .ge. q*artifa .and. s*xw .le. q*artifb) then
1539
1540 * Accept the polynomial fit.
1541
1542 if (abs( s*xw ) .ge. q*tol) then
1543 xtry = (s/q)*xw
1544 else
1545 xtry = zero
1546 end if
1547 if (debug) write (nout, 2500) xtry
1548 end if
1549 end if
1550 end if
1551 end if
1552
1553 * ==================================================================
1554
1555 if (.not. done) then
1556 alfa = alfbst + xtry
1557 if (braktd .or. alfa .lt. alfmax - tolmax) then
1558
1559 * The function must not be evaluated too close to a or b.
1560 * (It has already been evaluated at both those points.)
1561
1562 if (xtry .le. a + tol .or. xtry .ge. b - tol) then
1563 if (half*(a + b) .le. zero) then
1564 xtry = - tol
1565 else
1566 xtry = tol
1567 end if
1568 alfa = alfbst + xtry
1569 end if
1570 else
1571
1572 * The step is close to, or larger than alfmax, replace it by
1573 * alfmax to force evaluation of f at the boundary.
1574
1575 braktd = .true.
1576 xtry = alfmax - alfbst
1577 alfa = alfmax
1578 end if
1579 end if
1580
1581 * ------------------------------------------------------------------
1582 * Exit.
1583 * ------------------------------------------------------------------
1584 if (done) then
1585 if (badfun) then
1586 inform = 8
1587 else if (found) then
1588 if (alfbst .lt. alfmax) then
1589 inform = 1
1590 else
1591 inform = 2
1592 end if
1593 else if (moved ) then
1594 inform = 3
1595 else if (quitF) then
1596 inform = 7
1597 else if (crampd) then
1598 inform = 4
1599 else
1600 inform = 6
1601 end if
1602 end if
1603
1604 if (debug) write (nout, 3000)
1605 return
1606
1607 1000 format(/' g0 tolabs alfmax ', 1p, 2e22.14, e16.8
1608 $ /' targtg tolrel epsaf ', 1p, 2e22.14, e16.8
1609 $ /' crampd ', l3)
1610 1100 format(/' alfa ftry gtry ', 1p, 2e22.14, e16.8)
1611 1200 format(/' a b b - a tol ', 1p, 2e22.14, 2e16.8
1612 $ /' nsamea nsameb numf ', 3i3
1613 $ /' braktd extrap closef imprvd', 4l3
1614 $ /' found quitI ', 2l3
1615 $ /' alfbst fbest gbest ', 1p, 3e22.14
1616 $ /' alfaw fw gw ', 1p, 3e22.14)
1617 2100 format( ' Cubic. ')
1618 2200 format( ' Parabola.')
1619 2300 format( ' Bisection. xmidpt', 1p, e22.14)
1620 2400 format( ' Geo. bisection. xtry,daux,dtry', 1p, 3e22.14)
1621 2500 format( ' Polynomial fit accepted. xtry', 1p, e22.14)
1622 3000 format( ' ----------------------------------------------------'/)
1623
1624 * End of srchc .
1625 end
1626 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1627
1628 subroutine srchq ( first , debug , done , imprvd, inform,
1629 $ maxf , numf , nout ,
1630 $ alfmax, alfsml, epsaf ,
1631 $ g0 , targtg, ftry ,
1632 $ tolabs, tolrel, toltny,
1633 $ alfa , alfbst, fbest )
1634
1635 implicit double precision (a-h,o-z)
1636 logical first , debug , done , imprvd
1637
1638 ************************************************************************
1639 * srchq finds a sequence of improving estimates of a minimizer of
1640 * the univariate function f(alpha) in the interval (0,alfmax].
1641 * f(alpha) is a smooth function such that f(0) = 0 and f'(0) < 0.
1642 * srchq requires f(alpha) (but not f'(alpha)) to be evaluated
1643 * in the interval. New estimates of a minimizer are computed using
1644 * safeguarded quadratic interpolation.
1645 *
1646 * Reverse communication is used to allow the calling program to
1647 * evaluate f. Some of the parameters must be set or tested by the
1648 * calling program. The remainder would ordinarily be local
1649 * variables.
1650 *
1651 * Input parameters (relevant to the calling program)
1652 * --------------------------------------------------
1653 *
1654 * first must be true on the first entry. It is subsequently
1655 * altered by srchq.
1656 *
1657 * debug specifies whether detailed output is wanted.
1658 *
1659 * maxf is an upper limit on the number of times srchq is to
1660 * be entered consecutively with done = false
1661 * (following an initial entry with first = true).
1662 *
1663 * alfa is the first estimate of a minimizer. alfa is
1664 * subsequently altered by srchq (see below).
1665 *
1666 * alfmax is the upper limit of the interval to be searched.
1667 *
1668 * alfsml is intended to prevent inefficiency when a minimizer
1669 * is very small, for cases where the calling program
1670 * would prefer to redefine f'(alfa). alfsml is
1671 * allowed to be zero. Early termination will occur if
1672 * srchq determines that a minimizer lies somewhere in
1673 * the interval [0, alfsml) (but not if alfmax is
1674 * smaller that alfsml).
1675 *
1676 * epsaf is an estimate of the absolute precision in the
1677 * computed value of f(0).
1678 *
1679 * ftry the value of f at the new point
1680 * alfa = alfbst + xtry.
1681 *
1682 * g0 is the value of f'(0). g0 must be negative.
1683 *
1684 * tolabs,tolrel define a function tol(alfa) = tolrel*alfa + tolabs
1685 * such that if f has already been evaluated at alfa,
1686 * it will not be evaluated closer than tol(alfa).
1687 * These values may be reduced by srchc.
1688 *
1689 * targtg is the target value of abs(f'(alfa)). The search
1690 * is terminated when
1691 * abs(f'(alfa)) le targtg and f(alfa) lt 0.
1692 *
1693 * toltny is the smallest value that tolabs is allowed to be
1694 * reduced to.
1695 *
1696 * Output parameters (relevant to the calling program)
1697 * ---------------------------------------------------
1698 *
1699 * imprvd is true if the previous alfa was the best point so
1700 * far. Any related quantities should be saved by the
1701 * calling program (e.g., arrays) before paying
1702 * attention to the variable done.
1703 *
1704 * done = false means the calling program should evaluate ftry
1705 * for the new trial step alfa, and reenter srchq.
1706 *
1707 * done = true means that no new alfa was calculated. The value
1708 * of inform gives the result of the search as follows
1709 *
1710 * inform = 1 means the search has terminated
1711 * successfully with alfbst < alfmax.
1712 *
1713 * inform = 2 means the search has terminated
1714 * successfully with alfbst = alfmax.
1715 *
1716 * inform = 3 means that the search failed to find a
1717 * point of sufficient decrease in maxf
1718 * functions, but a lower point was found.
1719 *
1720 * inform = 4 means alfmax is so small that a search
1721 * should not have been attempted.
1722 *
1723 * inform = 5 means that the search was terminated
1724 * because of alfsml (see above).
1725 *
1726 * inform = 6 means the search has failed to find a
1727 * useful step. The interval of uncertainty
1728 * is [0,b] with b < 2*tolabs. A minimizer
1729 * lies very close to alfa = 0, or f'(0) is
1730 * not sufficiently accurate.
1731 *
1732 * inform = 7 if no better point could be found after
1733 * maxf function calls.
1734 *
1735 * inform = 8 means the input parameters were bad.
1736 * alfmax le toltny or g0 ge zero.
1737 * No function evaluations were made.
1738 *
1739 * numf counts the number of times srchq has been entered
1740 * consecutively with done = false (i.e., with a new
1741 * function value ftry).
1742 *
1743 * alfa is the point at which the next function ftry must
1744 * be computed.
1745 *
1746 * alfbst should be accepted by the calling program as the
1747 * approximate minimizer, whenever srchq returns
1748 * inform = 1, 2 or 3.
1749 *
1750 * fbest will be the corresponding value of f.
1751 *
1752 * The following parameters retain information between entries
1753 * -----------------------------------------------------------
1754 *
1755 * braktd is false if f has not been evaluated at the far end
1756 * of the interval of uncertainty. In this case, the
1757 * point b will be at alfmax + tol(alfmax).
1758 *
1759 * crampd is true if alfmax is very small (le tolabs). If the
1760 * search fails, this indicates that a zero step should
1761 * be taken.
1762 *
1763 * extrap is true if alfbst has moved at least once and xv
1764 * lies outside the interval of uncertainty. In this
1765 * case, extra safeguards are applied to allow for
1766 * instability in the polynomial fit.
1767 *
1768 * moved is true if a better point has been found, i.e.,
1769 * alfbst gt 0.
1770 *
1771 * vset records whether a third-best point has been defined.
1772 *
1773 * wset records whether a second-best point has been
1774 * defined. It will always be true by the time the
1775 * convergence test is applied.
1776 *
1777 * nsamea is the number of consecutive times that the
1778 * left-hand end point of the interval of uncertainty
1779 * has remained the same.
1780 *
1781 * nsameb similarly for the right-hand end.
1782 *
1783 * a, b, alfbst define the current interval of uncertainty.
1784 * A minimizer lies somewhere in the interval
1785 * [alfbst + a, alfbst + b].
1786 *
1787 * alfbst is the best point so far. It lies strictly within
1788 * [atrue,btrue] (except when alfbst has not been
1789 * moved, in which case it lies at the left-hand end
1790 * point). Hence we have a .le. 0 and b .gt. 0.
1791 *
1792 * fbest is the value of f at the point alfbst.
1793 *
1794 * fa is the value of f at the point alfbst + a.
1795 *
1796 * factor controls the rate at which extrapolated estimates of
1797 * alfa may expand into the interval of uncertainty.
1798 * Factor is not used if a minimizer has been bracketed
1799 * (i.e., when the variable braktd is true).
1800 *
1801 * fv, fw are the values of f at the points alfbst + xv and
1802 * alfbst + xw. They are not defined until vset or
1803 * wset are true.
1804 *
1805 * xtry is the trial point within the shifted interval
1806 * (a, b). The new trial function value must be
1807 * computed at the point alfa = alfbst + xtry.
1808 *
1809 * xv is such that alfbst + xv is the third-best point.
1810 * It is not defined until vset is true.
1811 *
1812 * xw is such that alfbst + xw is the second-best point.
1813 * It is not defined until wset is true. In some
1814 * cases, xw will replace a previous xw that has a
1815 * lower function but has just been excluded from
1816 * (a,b).
1817 *
1818 * Systems Optimization Laboratory, Stanford University, California.
1819 * Original version February 1982. Rev. May 1983.
1820 * Original F77 version 22-August-1985.
1821 * This version of srchq dated 24-Oct-91.
1822 ************************************************************************
1823
1824 logical braktd, crampd, extrap, moved , vset , wset
1825 save braktd, crampd, extrap, moved , vset , wset
1826
1827 save nsamea, nsameb
1828 save a , b , fa , factor
1829 save xtry , xw , fw , xv , fv , tolmax
1830
1831 logical badfun, closef, found
1832 logical quitF , quitFZ, quitI , quitS
1833 logical setxv , xinxw
1834 intrinsic abs , sqrt
1835
1836 parameter ( zero =0.0d+0, point1 =0.1d+0, half =0.5d+0 )
1837 parameter ( one =1.0d+0, two =2.0d+0, five =5.0d+0 )
1838 parameter ( ten =1.0d+1, eleven =1.1d+1 )
1839
1840 * ------------------------------------------------------------------
1841 * Local variables
1842 * ===============
1843 *
1844 * closef is true if the worst function fv is within epsaf of
1845 * fbest (up or down).
1846 *
1847 * found is true if the sufficient decrease conditions holds at
1848 * alfbst.
1849 *
1850 * quitF is true when maxf function calls have been made.
1851 *
1852 * quitFZ is true when the three best function values are within
1853 * epsaf of each other, and the new point satisfies
1854 * fbest le ftry le fbest+epsaf.
1855 *
1856 * quitI is true when the interval of uncertainty is less than
1857 * 2*tol.
1858 *
1859 * quitS is true as soon as alfa is too small to be useful;
1860 * i.e., btrue le alfsml.
1861 *
1862 * xinxw is true if xtry is in (xw,0) or (0,xw).
1863 * ------------------------------------------------------------------
1864
1865 imprvd = .false.
1866 badfun = .false.
1867 quitF = .false.
1868 quitFZ = .false.
1869 quitS = .false.
1870 quitI = .false.
1871
1872 if (first) then
1873 * ---------------------------------------------------------------
1874 * First entry. Initialize various quantities, check input data
1875 * and prepare to evaluate the function at the initial step alfa.
1876 * ---------------------------------------------------------------
1877 first = .false.
1878 numf = 0
1879 alfbst = zero
1880 badfun = alfmax .le. toltny .or. g0 .ge. zero
1881 done = badfun
1882 moved = .false.
1883
1884 if (.not. done) then
1885 braktd = .false.
1886 crampd = alfmax .le. tolabs
1887 extrap = .false.
1888 vset = .false.
1889 wset = .false.
1890 nsamea = 0
1891 nsameb = 0
1892
1893 tolmax = tolrel*alfmax + tolabs
1894 a = zero
1895 b = alfmax + tolmax
1896 fa = zero
1897 factor = five
1898 tol = tolabs
1899 xtry = alfa
1900 if (debug) then
1901 write (nout, 1000) g0 , tolabs, alfmax,
1902 $ targtg, tolrel, epsaf , crampd
1903 end if
1904 end if
1905 else
1906 * ---------------------------------------------------------------
1907 * Subsequent entries. The function has just been evaluated at
1908 * alfa = alfbst + xtry, giving ftry.
1909 * ---------------------------------------------------------------
1910 if (debug) write (nout, 1100) alfa, ftry
1911
1912 numf = numf + 1
1913 nsamea = nsamea + 1
1914 nsameb = nsameb + 1
1915
1916 if (.not. braktd) then
1917 tolmax = tolabs + tolrel*alfmax
1918 b = alfmax - alfbst + tolmax
1919 end if
1920
1921 * Check if xtry is in the interval (xw,0) or (0,xw).
1922
1923 if (wset) then
1924 xinxw = zero .lt. xtry .and. xtry .le. xw
1925 $ .or. xw .le. xtry .and. xtry .lt. zero
1926 else
1927 xinxw = .false.
1928 end if
1929
1930 imprvd = ftry .lt. fbest
1931 if (vset) then
1932 closef = abs( fbest - fv ) .le. epsaf
1933 else
1934 closef = .false.
1935 end if
1936
1937 if (imprvd) then
1938
1939 * We seem to have an improvement. The new point becomes the
1940 * origin and other points are shifted accordingly.
1941
1942 if (wset) then
1943 xv = xw - xtry
1944 fv = fw
1945 vset = .true.
1946 end if
1947
1948 xw = zero - xtry
1949 fw = fbest
1950 wset = .true.
1951 fbest = ftry
1952 alfbst = alfa
1953 moved = .true.
1954
1955 a = a - xtry
1956 b = b - xtry
1957 extrap = .not. xinxw
1958
1959 * Decrease the length of (a,b).
1960
1961 if (xtry .ge. zero) then
1962 a = xw
1963 fa = fw
1964 nsamea = 0
1965 else
1966 b = xw
1967 nsameb = 0
1968 braktd = .true.
1969 end if
1970 else if (closef .and. ftry - fbest .lt. epsaf) then
1971
1972 * Quit if there has been no progress and ftry, fbest, fw
1973 * and fv are all within epsaf of each other.
1974
1975 quitFZ = .true.
1976 else
1977
1978 * The new function value is no better than the current best
1979 * point. xtry must an end point of the new (a,b).
1980
1981 if (xtry .lt. zero) then
1982 a = xtry
1983 fa = ftry
1984 nsamea = 0
1985 else
1986 b = xtry
1987 nsameb = 0
1988 braktd = .true.
1989 end if
1990
1991 * The origin remains unchanged but xtry may qualify as xw.
1992
1993 if (wset) then
1994 if (ftry .lt. fw) then
1995 xv = xw
1996 fv = fw
1997 vset = .true.
1998
1999 xw = xtry
2000 fw = ftry
2001 if (moved) extrap = xinxw
2002 else if (moved) then
2003 if (vset) then
2004 setxv = ftry .lt. fv .or. .not. extrap
2005 else
2006 setxv = .true.
2007 end if
2008
2009 if (setxv) then
2010 if (vset .and. xinxw) then
2011 xw = xv
2012 fw = fv
2013 end if
2014 xv = xtry
2015 fv = ftry
2016 vset = .true.
2017 end if
2018 else
2019 xw = xtry
2020 fw = ftry
2021 end if
2022 else
2023 xw = xtry
2024 fw = ftry
2025 wset = .true.
2026 end if
2027 end if
2028
2029 * ---------------------------------------------------------------
2030 * Check the termination criteria.
2031 * ---------------------------------------------------------------
2032 tol = tolabs + tolrel*alfbst
2033 truea = alfbst + a
2034 trueb = alfbst + b
2035
2036 found = moved .and. abs(fa - fbest) .le. -a*targtg
2037 quitF = numf .ge. maxf
2038 quitI = b - a .le. tol + tol
2039 quitS = trueb .le. alfsml
2040
2041 if (quitI .and. .not. moved) then
2042
2043 * The interval of uncertainty appears to be small enough,
2044 * but no better point has been found. Check that changing
2045 * alfa by b-a changes f by less than epsaf.
2046
2047 tol = tol/ten
2048 tolabs = tol
2049 quitI = abs(fw) .le. epsaf .or. tol .le. toltny
2050 end if
2051
2052 done = quitF .or. quitFZ .or. quitS .or. quitI
2053 $ .or. found
2054
2055 if (debug) then
2056 write (nout, 1200) truea , trueb , b-a , tol ,
2057 $ nsamea , nsameb, numf ,
2058 $ braktd , extrap, closef, imprvd,
2059 $ found , quitI , quitFZ, quitS ,
2060 $ alfbst , fbest ,
2061 $ alfbst+xw, fw
2062 if (vset) then
2063 write (nout, 1300) alfbst + xv, fv
2064 end if
2065 end if
2066
2067 * ---------------------------------------------------------------
2068 * Proceed with the computation of an estimate of a minimizer.
2069 * The choices are...
2070 * 1. Parabolic fit using function values only.
2071 * 2. Damped parabolic fit if the regular fit appears to be
2072 * consistently overestimating the distance to a minimizer.
2073 * 3. Bisection, geometric bisection, or a step of tol if the
2074 * parabolic fit is unsatisfactory.
2075 * ---------------------------------------------------------------
2076 if (.not. done) then
2077 xmidpt = half*(a + b)
2078 s = zero
2079 q = zero
2080
2081 * ============================================================
2082 * Fit a parabola.
2083 * ============================================================
2084 * See if there are two or three points for the parabolic fit.
2085
2086 gw = (fw - fbest)/xw
2087 if (vset .and. moved) then
2088
2089 * Three points available. Use fbest, fw and fv.
2090
2091 gv = (fv - fbest)/xv
2092 s = gv - (xv/xw)*gw
2093 q = two*(gv - gw)
2094 if (debug) write (nout, 2200)
2095 else
2096
2097 * Only two points available. Use fbest, fw and g0.
2098
2099 if (moved) then
2100 s = g0 - two*gw
2101 else
2102 s = g0
2103 end if
2104 q = two*(g0 - gw)
2105 if (debug) write (nout, 2100)
2106 end if
2107
2108 * ------------------------------------------------------------
2109 * Construct an artificial interval (artifa, artifb) in which
2110 * the new estimate of the steplength must lie. Set a default
2111 * value of xtry that will be used if the polynomial fit is
2112 * rejected. In the following, the interval (a,b) is considered
2113 * the sum of two intervals of lengths dtry and daux, with
2114 * common end point the best point (zero). dtry is the length
2115 * of the interval into which the default xtry will be placed
2116 * and endpnt denotes its non-zero end point. The magnitude of
2117 * xtry is computed so that the exponents of dtry and daux are
2118 * approximately bisected.
2119 * ------------------------------------------------------------
2120 artifa = a
2121 artifb = b
2122 if (.not. braktd) then
2123
2124 * A minimizer has not yet been bracketed.
2125 * Set an artificial upper bound by expanding the interval
2126 * xw by a suitable factor.
2127
2128 xtry = - factor*xw
2129 artifb = xtry
2130 if (alfbst + xtry .lt. alfmax) factor = five*factor
2131 else if (vset .and. moved) then
2132
2133 * Three points exist in the interval of uncertainty.
2134 * Check if the points are configured for an extrapolation
2135 * or an interpolation.
2136
2137 if (extrap) then
2138
2139 * The points are configured for an extrapolation.
2140
2141 if (xw .lt. zero) endpnt = b
2142 if (xw .gt. zero) endpnt = a
2143 else
2144
2145 * If the interpolation appears to be overestimating the
2146 * distance to a minimizer, damp the interpolation step.
2147
2148 if (nsamea .ge. 3 .or. nsameb .ge. 3) then
2149 factor = factor / five
2150 s = factor * s
2151 else
2152 factor = one
2153 end if
2154
2155 * The points are configured for an interpolation. The
2156 * artificial interval will be just (a,b). Set endpnt so
2157 * that xtry lies in the larger of the intervals (a,b)
2158 * and (0,b).
2159
2160 if (xmidpt .gt. zero) then
2161 endpnt = b
2162 else
2163 endpnt = a
2164 end if
2165
2166 * If a bound has remained the same for three iterations,
2167 * set endpnt so that xtry is likely to replace the
2168 * offending bound.
2169
2170 if (nsamea .ge. 3) endpnt = a
2171 if (nsameb .ge. 3) endpnt = b
2172 end if
2173
2174 * Compute the default value of xtry.
2175
2176 dtry = abs( endpnt )
2177 daux = b - a - dtry
2178 if (daux .ge. dtry) then
2179 xtry = five*dtry*(point1 + dtry/daux)/eleven
2180 else
2181 xtry = half*sqrt( daux )*sqrt( dtry )
2182 end if
2183 if (endpnt .lt. zero) xtry = - xtry
2184 if (debug) write (nout, 2500) xtry, daux, dtry
2185
2186 * If the points are configured for an extrapolation set the
2187 * artificial bounds so that the artificial interval lies
2188 * within (a,b). If the polynomial fit is rejected, xtry
2189 * will remain at the relevant artificial bound.
2190
2191 if (extrap) then
2192 if (xtry .le. zero) then
2193 artifa = xtry
2194 else
2195 artifb = xtry
2196 end if
2197 end if
2198 else
2199
2200 * The gradient at the origin is being used for the
2201 * polynomial fit. Set the default xtry to one tenth xw.
2202
2203 if (extrap) then
2204 xtry = - xw
2205 else
2206 xtry = xw/ten
2207 end if
2208 if (debug) write (nout, 2400) xtry
2209 end if
2210
2211 * ------------------------------------------------------------
2212 * The polynomial fits give (s/q)*xw as the new step. Reject
2213 * this step if it lies outside (artifa, artifb).
2214 * ------------------------------------------------------------
2215 if (q .ne. zero) then
2216 if (q .lt. zero) s = - s
2217 if (q .lt. zero) q = - q
2218 if (s*xw .ge. q*artifa .and. s*xw .le. q*artifb) then
2219
2220 * Accept the polynomial fit.
2221
2222 if (abs( s*xw ) .ge. q*tol) then
2223 xtry = (s/q)*xw
2224 else
2225 xtry = zero
2226 end if
2227 if (debug) write (nout, 2600) xtry
2228 end if
2229 end if
2230 end if
2231 end if
2232 * ==================================================================
2233
2234 if (.not. done) then
2235 alfa = alfbst + xtry
2236 if (braktd .or. alfa .lt. alfmax - tolmax) then
2237
2238 * The function must not be evaluated too close to a or b.
2239 * (It has already been evaluated at both those points.)
2240
2241 xmidpt = half*(a + b)
2242 if (xtry .le. a + tol .or. xtry .ge. b - tol) then
2243 if (xmidpt .le. zero) then
2244 xtry = - tol
2245 else
2246 xtry = tol
2247 end if
2248 end if
2249
2250 if (abs( xtry ) .lt. tol) then
2251 if (xmidpt .le. zero) then
2252 xtry = - tol
2253 else
2254 xtry = tol
2255 end if
2256 end if
2257 alfa = alfbst + xtry
2258 else
2259
2260 * The step is close to or larger than alfmax, replace it by
2261 * alfmax to force evaluation of the function at the boundary.
2262
2263 braktd = .true.
2264 xtry = alfmax - alfbst
2265 alfa = alfmax
2266 end if
2267 end if
2268 * ------------------------------------------------------------------
2269 * Exit.
2270 * ------------------------------------------------------------------
2271 if (done) then
2272 if (badfun) then
2273 inform = 8
2274 else if (quitS ) then
2275 inform = 5
2276 else if (found) then
2277 if (alfbst .lt. alfmax) then
2278 inform = 1
2279 else
2280 inform = 2
2281 end if
2282 else if (moved ) then
2283 inform = 3
2284 else if (quitF) then
2285 inform = 7
2286 else if (crampd) then
2287 inform = 4
2288 else
2289 inform = 6
2290 end if
2291 end if
2292
2293 if (debug) write (nout, 3000)
2294 return
2295
2296 1000 format(/' g0 tolabs alfmax ', 1p, 2e22.14, e16.8
2297 $ /' targtg tolrel epsaf ', 1p, 2e22.14, e16.8
2298 $ /' crampd ', l3)
2299 1100 format(/' alfa ftry ', 1p,2e22.14 )
2300 1200 format(/' a b b - a tol ', 1p,2e22.14, 2e16.8
2301 $ /' nsamea nsameb numf ', 3i3
2302 $ /' braktd extrap closef imprvd', 4l3
2303 $ /' found quitI quitFZ quitS ', 4l3
2304 $ /' alfbst fbest ', 1p,2e22.14
2305 $ /' alfaw fw ', 1p,2e22.14)
2306 1300 format( ' alfav fv ', 1p,2e22.14 /)
2307 2100 format( ' Parabolic fit, two points. ')
2308 2200 format( ' Parabolic fit, three points. ')
2309 2400 format( ' Exponent reduced. Trial point', 1p, e22.14)
2310 2500 format( ' Geo. bisection. xtry,daux,dtry', 1p, 3e22.14)
2311 2600 format( ' Polynomial fit accepted. xtry', 1p, e22.14)
2312 3000 format( ' ----------------------------------------------------'/)
2313
2314 * End of srchq .
2315 end

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22