| 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 |
Parent Directory
| 
Revision Log