| 1 | ! **************************************************************
|
|---|
| 2 | !
|
|---|
| 3 | ! This file contains the subroutines: mincjg
|
|---|
| 4 | !
|
|---|
| 5 | ! Copyright 2003-2005 Frank Eisenmenger, U.H.E. Hansmann,
|
|---|
| 6 | ! Shura Hayryan, Chin-Ku
|
|---|
| 7 | ! Copyright 2007 Frank Eisenmenger, U.H.E. Hansmann,
|
|---|
| 8 | ! Jan H. Meinke, Sandipan Mohanty
|
|---|
| 9 | !
|
|---|
| 10 | ! ********************************************************************
|
|---|
| 11 | subroutine mincjg(n,mxn,x,f,g,acur,d,xa,ga,dt,yt,gt,maxfun,nfun)
|
|---|
| 12 |
|
|---|
| 13 | ! ....................................................................
|
|---|
| 14 | !
|
|---|
| 15 | ! Conjugate Gradient Minimizer
|
|---|
| 16 | !
|
|---|
| 17 | ! INPUT: X,F,G - variables, value of FUNC, gradient at START/
|
|---|
| 18 | ! ACUR - convergence is assumed if ACUR > SUM ( G(I)**2 )
|
|---|
| 19 | ! MAXFUN - maximum overall number of function calls
|
|---|
| 20 | !
|
|---|
| 21 | ! OUTPUT: X,F,G - variables, value of FUNC, gradient at MINIMUM
|
|---|
| 22 | ! NFUN - overall number of function calls used
|
|---|
| 23 | !
|
|---|
| 24 | ! ARRAYS: D,XA,GA,YT,DT,GT - dimension N
|
|---|
| 25 | !
|
|---|
| 26 | ! CALLS: MOVE - calculate function & its gradients for current X
|
|---|
| 27 | !
|
|---|
| 28 | ! PARAMETERS: AMF - rough estimate of first reduction in F, used
|
|---|
| 29 | ! to guess initial step of 1st line search
|
|---|
| 30 | ! MXFCON - see 'ier=4'
|
|---|
| 31 | ! MAXLIN -
|
|---|
| 32 | !
|
|---|
| 33 | ! DIAGNOSTICS (ier)
|
|---|
| 34 | !
|
|---|
| 35 | ! = 0: minimization completed successfully
|
|---|
| 36 | ! = 1: number of steps reached MAXFUN
|
|---|
| 37 | ! = 2: line search was abandoned
|
|---|
| 38 | ! = 3: search direction is uphill
|
|---|
| 39 | ! = 4: two consecutive line searches failed to reduce F
|
|---|
| 40 | ! ....................................................................
|
|---|
| 41 |
|
|---|
| 42 | double precision amf, tol, eps, gsqrd, fmin, f, t, g, ga, d, xa, x
|
|---|
| 43 | double precision gsq2, gnew, dfpr, stmin, fuit, gdit, gt, gdmi
|
|---|
| 44 | double precision sbound, stepch, step, wo, sum, fch, acur, ddspln
|
|---|
| 45 | double precision gspln, beta, gama, yt, dt, gamden, di
|
|---|
| 46 |
|
|---|
| 47 | integer mxfcon, maxlin, mxn, ier, iter, iterfm, iterrs, nfun
|
|---|
| 48 | integer nfopt, i, n, nfbeg, iretry, maxfun
|
|---|
| 49 |
|
|---|
| 50 | parameter (AMF = 10.d0,
|
|---|
| 51 | & MXFCON = 2,
|
|---|
| 52 | & MAXLIN = 5,
|
|---|
| 53 | & TOL = 1.d-7, ! controls 'stepch'
|
|---|
| 54 | & EPS = .7d0)
|
|---|
| 55 |
|
|---|
| 56 | dimension x(mxn),g(mxn),
|
|---|
| 57 | & d(mxn),xa(mxn),ga(mxn),dt(mxn),yt(mxn),gt(mxn)
|
|---|
| 58 |
|
|---|
| 59 | character(255) logString
|
|---|
| 60 |
|
|---|
| 61 | ier = 0
|
|---|
| 62 | iter = 0
|
|---|
| 63 | iterfm = 0
|
|---|
| 64 | iterrs = 0
|
|---|
| 65 | nfun = 0
|
|---|
| 66 | nfopt = nfun
|
|---|
| 67 |
|
|---|
| 68 | gsqrd = 0.d0
|
|---|
| 69 | fmin = f
|
|---|
| 70 |
|
|---|
| 71 | do i=1,n
|
|---|
| 72 | t = g(i)
|
|---|
| 73 | ga(i) = t
|
|---|
| 74 | d(i) = -t
|
|---|
| 75 | gsqrd = gsqrd + t**2
|
|---|
| 76 | xa(i) = x(i)
|
|---|
| 77 | enddo
|
|---|
| 78 |
|
|---|
| 79 | gsq2 = .2d0 * gsqrd
|
|---|
| 80 | gnew = -gsqrd
|
|---|
| 81 | dfpr = AMF
|
|---|
| 82 | stmin = AMF / gsqrd
|
|---|
| 83 |
|
|---|
| 84 | 1 iter = iter + 1
|
|---|
| 85 |
|
|---|
| 86 | fuit = f
|
|---|
| 87 | gdit = 0.d0
|
|---|
| 88 |
|
|---|
| 89 | do i=1,n
|
|---|
| 90 | t = g(i)
|
|---|
| 91 | gt(i) = t
|
|---|
| 92 | gdit = gdit + d(i) * t
|
|---|
| 93 | enddo
|
|---|
| 94 |
|
|---|
| 95 | if ( gdit .ge. 0.d0 ) then
|
|---|
| 96 | write (logString, *) ' mincjg> search direction is uphill'
|
|---|
| 97 | ier = 3
|
|---|
| 98 | goto 6
|
|---|
| 99 | endif
|
|---|
| 100 |
|
|---|
| 101 | gdmi = gdit
|
|---|
| 102 | sbound = -1.d0
|
|---|
| 103 | nfbeg = nfun
|
|---|
| 104 | iretry = -1
|
|---|
| 105 |
|
|---|
| 106 | stepch = min( stmin, abs ( (dfpr/gdit) ) )
|
|---|
| 107 | stmin = 0.d0
|
|---|
| 108 |
|
|---|
| 109 | 2 step = stmin + stepch
|
|---|
| 110 | wo = 0.d0
|
|---|
| 111 |
|
|---|
| 112 | do i=1,n
|
|---|
| 113 | t = stepch * d(i)
|
|---|
| 114 | wo = max( wo, abs( t ) )
|
|---|
| 115 | x(i) = xa(i) + t
|
|---|
| 116 | enddo
|
|---|
| 117 |
|
|---|
| 118 | if ( wo .gt. TOL ) then
|
|---|
| 119 |
|
|---|
| 120 | nfun = nfun + 1
|
|---|
| 121 | call move(nfun,n,f,x,g)
|
|---|
| 122 |
|
|---|
| 123 | gnew = 0.d0
|
|---|
| 124 | sum = 0.d0
|
|---|
| 125 |
|
|---|
| 126 | do i=1,n
|
|---|
| 127 | t = g(i)
|
|---|
| 128 | gnew = gnew + d(i) * t
|
|---|
| 129 | sum = sum + t**2
|
|---|
| 130 | enddo
|
|---|
| 131 |
|
|---|
| 132 | fch = f - fmin
|
|---|
| 133 |
|
|---|
| 134 | if ( fch .le. 0.d0 ) then
|
|---|
| 135 |
|
|---|
| 136 | if ( fch .lt. 0.d0 .or. (gnew/gdmi) .ge. -1.d0 ) then
|
|---|
| 137 |
|
|---|
| 138 | fmin = f
|
|---|
| 139 | gsqrd = sum
|
|---|
| 140 | gsq2 = .2d0 * gsqrd
|
|---|
| 141 | nfopt = nfun
|
|---|
| 142 |
|
|---|
| 143 | do i=1,n
|
|---|
| 144 | xa(i) = x(i)
|
|---|
| 145 | ga(i) = g(i)
|
|---|
| 146 | enddo
|
|---|
| 147 |
|
|---|
| 148 | endif
|
|---|
| 149 |
|
|---|
| 150 | if ( sum .le. ACUR ) return ! normal end
|
|---|
| 151 |
|
|---|
| 152 | endif
|
|---|
| 153 |
|
|---|
| 154 | if ( nfun .eq. MAXFUN ) then
|
|---|
| 155 | ier = 1
|
|---|
| 156 | return
|
|---|
| 157 | endif
|
|---|
| 158 |
|
|---|
| 159 | else ! stepch is effectively zero
|
|---|
| 160 |
|
|---|
| 161 | if ( nfun .gt. (nfbeg + 1) .or.
|
|---|
| 162 | & abs(gdmi/gdit) .gt. EPS ) then
|
|---|
| 163 |
|
|---|
| 164 | ier=2
|
|---|
| 165 | write (logString, *)
|
|---|
| 166 | & ' mincjg> too small step in search direction'
|
|---|
| 167 | endif
|
|---|
| 168 |
|
|---|
| 169 | goto 6
|
|---|
| 170 |
|
|---|
| 171 | endif
|
|---|
| 172 |
|
|---|
| 173 | wo = (fch + fch) / stepch - gnew - gdmi
|
|---|
| 174 | ddspln = (gnew - gdmi) / stepch
|
|---|
| 175 |
|
|---|
| 176 | if ( nfun .le. nfopt ) then
|
|---|
| 177 |
|
|---|
| 178 | if ( gdmi * gnew .le. 0.d0 ) sbound = stmin
|
|---|
| 179 |
|
|---|
| 180 | stmin = step
|
|---|
| 181 | gdmi = gnew
|
|---|
| 182 | stepch = -stepch
|
|---|
| 183 |
|
|---|
| 184 | else
|
|---|
| 185 |
|
|---|
| 186 | sbound = step
|
|---|
| 187 |
|
|---|
| 188 | endif
|
|---|
| 189 |
|
|---|
| 190 | if ( fch .ne. 0.d0 ) ddspln = ddspln + (wo + wo) / stepch
|
|---|
| 191 |
|
|---|
| 192 | if ( gdmi .eq. 0.d0 ) goto 6
|
|---|
| 193 |
|
|---|
| 194 | if ( nfun .le. (nfbeg + 1) ) goto 4
|
|---|
| 195 |
|
|---|
| 196 | if ( abs(gdmi/gdit) .le. EPS ) goto 6
|
|---|
| 197 |
|
|---|
| 198 | 3 if ( nfun .ge. (nfopt + MAXLIN) ) then
|
|---|
| 199 |
|
|---|
| 200 | ier = 2
|
|---|
| 201 | goto 6
|
|---|
| 202 |
|
|---|
| 203 | endif
|
|---|
| 204 |
|
|---|
| 205 | 4 if ( sbound .lt. -.5d0 ) then
|
|---|
| 206 | stepch = 9.d0 * stmin
|
|---|
| 207 | else
|
|---|
| 208 | stepch = .5d0 * ( sbound - stmin )
|
|---|
| 209 | endif
|
|---|
| 210 |
|
|---|
| 211 | gspln = gdmi + stepch * ddspln
|
|---|
| 212 |
|
|---|
| 213 | if ( (gdmi * gspln) .lt. 0.d0 ) stepch = stepch * gdmi /
|
|---|
| 214 | & (gdmi - gspln)
|
|---|
| 215 |
|
|---|
| 216 | goto 2
|
|---|
| 217 |
|
|---|
| 218 | 5 sum = 0.d0
|
|---|
| 219 | do i=1,n
|
|---|
| 220 | sum = sum + g(i) * gt(i)
|
|---|
| 221 | enddo
|
|---|
| 222 |
|
|---|
| 223 | beta = (gsqrd - sum) / (gdmi - gdit)
|
|---|
| 224 |
|
|---|
| 225 | if ( abs(beta * gdmi) .gt. gsq2) then
|
|---|
| 226 | iretry = iretry + 1
|
|---|
| 227 | if (iretry .le. 0) goto 3
|
|---|
| 228 | endif
|
|---|
| 229 |
|
|---|
| 230 | if ( f .lt. fuit ) iterfm = iter
|
|---|
| 231 |
|
|---|
| 232 | if ( iter .ge. (iterfm + MXFCON) ) then
|
|---|
| 233 | ier = 4
|
|---|
| 234 | write (logString, *)
|
|---|
| 235 | & ' mincjg> line search failed to reduce function'
|
|---|
| 236 | return
|
|---|
| 237 | endif
|
|---|
| 238 |
|
|---|
| 239 | dfpr = stmin * gdit
|
|---|
| 240 |
|
|---|
| 241 | if ( iretry .gt. 0 ) then
|
|---|
| 242 |
|
|---|
| 243 | do i=1,n
|
|---|
| 244 | d(i) = -g(i)
|
|---|
| 245 | enddo
|
|---|
| 246 |
|
|---|
| 247 | iterrs = 0
|
|---|
| 248 |
|
|---|
| 249 | else
|
|---|
| 250 |
|
|---|
| 251 | if (iterrs .ne. 0 .and. (iter - iterrs) .lt. (n-1) .and.
|
|---|
| 252 | & abs(sum) .lt. gsq2 ) then
|
|---|
| 253 |
|
|---|
| 254 | gama = 0.d0
|
|---|
| 255 | sum = 0.d0
|
|---|
| 256 |
|
|---|
| 257 | do i=1,n
|
|---|
| 258 | t = g(i)
|
|---|
| 259 | gama = gama + t * yt(i)
|
|---|
| 260 | sum = sum + t * dt(i)
|
|---|
| 261 | enddo
|
|---|
| 262 |
|
|---|
| 263 | gama = gama / gamden
|
|---|
| 264 |
|
|---|
| 265 | if ( abs( (beta * gdmi + gama * sum) ) .lt. gsq2) then
|
|---|
| 266 |
|
|---|
| 267 | do i=1,n
|
|---|
| 268 | d(i) = -g(i) + beta * d(i) + gama * dt(i)
|
|---|
| 269 | enddo
|
|---|
| 270 |
|
|---|
| 271 | goto 1
|
|---|
| 272 |
|
|---|
| 273 | endif
|
|---|
| 274 |
|
|---|
| 275 | endif
|
|---|
| 276 |
|
|---|
| 277 | gamden = gdmi - gdit
|
|---|
| 278 |
|
|---|
| 279 | do i=1,n
|
|---|
| 280 | t = g(i)
|
|---|
| 281 | di = d(i)
|
|---|
| 282 | dt(i) = di
|
|---|
| 283 | yt(i) = t - gt(i)
|
|---|
| 284 | d(i) = -t + beta * di
|
|---|
| 285 | enddo
|
|---|
| 286 |
|
|---|
| 287 | iterrs = iter
|
|---|
| 288 |
|
|---|
| 289 | endif
|
|---|
| 290 |
|
|---|
| 291 | goto 1
|
|---|
| 292 |
|
|---|
| 293 | 6 if ( nfun .ne. nfopt ) then
|
|---|
| 294 |
|
|---|
| 295 | f = fmin
|
|---|
| 296 |
|
|---|
| 297 | do i=1,n
|
|---|
| 298 | x(i) = xa(i)
|
|---|
| 299 | g(i) = ga(i)
|
|---|
| 300 | enddo
|
|---|
| 301 |
|
|---|
| 302 | endif
|
|---|
| 303 |
|
|---|
| 304 | if ( ier .eq. 0 ) goto 5
|
|---|
| 305 |
|
|---|
| 306 | nfun = nfun + 1
|
|---|
| 307 | call move(nfun,n,f,x,g)
|
|---|
| 308 |
|
|---|
| 309 | write (logString, *) ' mincjg> ier = ',ier
|
|---|
| 310 |
|
|---|
| 311 | return
|
|---|
| 312 | end
|
|---|
| 313 |
|
|---|
