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 |
|
---|