-
Notifications
You must be signed in to change notification settings - Fork 4
/
mnmigr.F
323 lines (323 loc) · 10.9 KB
/
mnmigr.F
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
*
* $Id: mnmigr.F,v 1.1.1.1 2000/06/08 11:19:20 andras Exp $
*
* $Log: mnmigr.F,v $
* Revision 1.1.1.1 2000/06/08 11:19:20 andras
* import of MINUIT from CERNlib 2000
*
* Revision 1.2 1996/03/15 18:02:49 james
* Modified Files:
* mnderi.F eliminate possible division by zero
* mnexcm.F suppress print on STOP when print flag=-1
* set FVAL3 to flag if FCN already called with IFLAG=3
* mninit.F set version 96.03
* mnlims.F remove arguments, not needed
* mnmigr.F VLEN -> LENV in debug print statement
* mnparm.F move call to MNRSET to after NPAR redefined, to zero all
* mnpsdf.F eliminate possible division by zero
* mnscan.F suppress printout when print flag =-1
* mnset.F remove arguments in call to MNLIMS
* mnsimp.F fix CSTATU so status is PROGRESS only if new minimum
* mnvert.F eliminate possible division by zero
*
* Revision 1.1.1.1 1996/03/07 14:31:30 mclareni
* Minuit
*
*
#include "minuit/pilot.h"
SUBROUTINE MNMIGR(FCN,FUTIL)
#include "minuit/d506dp.inc"
CC Performs a local function minimization using basically the
CC method of Davidon-Fletcher-Powell as modified by Fletcher
CC ref. -- Fletcher, Comp.J. 13,317 (1970) "switching method"
CC
#include "minuit/d506cm.inc"
EXTERNAL FCN,FUTIL
DIMENSION GS(MNI), STEP(MNI), XXS(MNI), FLNU(MNI), VG(MNI)
LOGICAL LDEBUG
PARAMETER (TOLER=0.05)
IF (NPAR .LE. 0) RETURN
IF (AMIN .EQ. UNDEFI) CALL MNAMIN(FCN,FUTIL)
LDEBUG = (IDBG(4) .GE. 1)
CFROM = 'MIGRAD '
NFCNFR = NFCN
NFCNMG = NFCN
CSTATU= 'INITIATE '
ISWTR = ISW(5) - 2*ITAUR
NPFN = NFCN
NPARX = NPAR
LENV = NPAR*(NPAR+1)/2
NRSTRT = 0
NPSDF = 0
LINED2 = 0
ISW(4) = -1
RHOTOL = 1.0E-3*APSI
IF (ISWTR .GE. 1) WRITE (ISYSWR,470) ISTRAT,RHOTOL
470 FORMAT (' START MIGRAD MINIMIZATION. STRATEGY',I2,
+'. CONVERGENCE WHEN EDM .LT.',E9.2)
C initialization strategy
IF (ISTRAT.LT.2 .OR. ISW(2).GE.3) GO TO 2
C come (back) here to restart completely
1 CONTINUE
IF (NRSTRT .GT. ISTRAT) THEN
CSTATU= 'FAILED '
ISW(4) = -1
GO TO 230
ENDIF
C . get full covariance and gradient
CALL MNHESS(FCN,FUTIL)
CALL MNWERR
NPSDF = 0
IF (ISW(2) .GE. 1) GO TO 10
C . get gradient at start point
2 CONTINUE
CALL MNINEX(X)
IF (ISW(3) .EQ. 1) THEN
CALL FCN(NPARX,GIN,FZERO,U,2,FUTIL)
NFCN = NFCN + 1
ENDIF
CALL MNDERI(FCN,FUTIL)
IF (ISW(2) .GE. 1) GO TO 10
C sometimes start with diagonal matrix
DO 3 I= 1, NPAR
XXS(I) = X(I)
STEP(I) = ZERO
3 CONTINUE
C do line search if second derivative negative
LINED2 = LINED2 + 1
IF (LINED2 .LT. (ISTRAT+1)*NPAR) THEN
DO 5 I= 1, NPAR
IF (G2(I) .GT. ZERO) GO TO 5
STEP(I) = -SIGN(GSTEP(I),GRD(I))
GDEL = STEP(I)*GRD(I)
FS = AMIN
CALL MNLINE(FCN,XXS,FS,STEP,GDEL,TOLER,FUTIL)
CALL MNWARN('D','MNMIGR','Negative G2 line search')
IEXT = NEXOFI(I)
IF (LDEBUG) WRITE (ISYSWR,'(A,I3,2G13.3)')
+ ' Negative G2 line search, param ',IEXT,FS,AMIN
GO TO 2
5 CONTINUE
ENDIF
C make diagonal error matrix
DO 8 I=1,NPAR
NDEX = I*(I-1)/2
DO 7 J=1,I-1
NDEX = NDEX + 1
7 VHMAT(NDEX) = 0.
NDEX = NDEX + 1
IF (G2(I) .LE. ZERO) G2(I) = 1.
VHMAT(NDEX) = 2./G2(I)
8 CONTINUE
DCOVAR = 1.
IF (LDEBUG) WRITE (ISYSWR,'(A,A/(1X,10G10.2))') ' DEBUG MNMIGR,',
+ ' STARTING MATRIX DIAGONAL, VHMAT=', (VHMAT(KK),KK=1,LENV)
C ready to start first iteration
10 CONTINUE
NRSTRT = NRSTRT + 1
IF (NRSTRT .GT. ISTRAT+1) THEN
CSTATU= 'FAILED '
GO TO 230
ENDIF
FS = AMIN
C . . . get EDM and set up loop
EDM = 0.
DO 18 I= 1, NPAR
GS(I) = GRD(I)
XXS(I) = X(I)
NDEX = I*(I-1)/2
DO 17 J= 1, I-1
NDEX = NDEX + 1
17 EDM = EDM + GS(I)*VHMAT(NDEX)*GS(J)
NDEX = NDEX + 1
18 EDM = EDM + 0.5 * GS(I)**2 *VHMAT(NDEX)
EDM = EDM * 0.5 * (1.0+3.0*DCOVAR)
IF (EDM .LT. ZERO) THEN
CALL MNWARN('W','MIGRAD','STARTING MATRIX NOT POS-DEFINITE.')
ISW(2) = 0
DCOVAR = 1.
GO TO 2
ENDIF
IF (ISW(2) .EQ. 0) EDM=BIGEDM
ITER = 0
CALL MNINEX(X)
CALL MNWERR
IF (ISWTR .GE. 1) CALL MNPRIN(3,AMIN)
IF (ISWTR .GE. 2) CALL MNMATU(0)
C . . . . . start main loop
24 CONTINUE
IF (NFCN-NPFN .GE. NFCNMX) GO TO 190
GDEL = 0.
GSSQ = 0.
DO 30 I=1,NPAR
RI = 0.
GSSQ = GSSQ + GS(I)**2
DO 25 J=1,NPAR
M = MAX(I,J)
N = MIN(I,J)
NDEX = M*(M-1)/2 + N
25 RI = RI + VHMAT(NDEX) *GS(J)
STEP(I) = -0.5*RI
30 GDEL = GDEL + STEP(I)*GS(I)
IF (GSSQ .EQ. ZERO) THEN
CALL MNWARN('D','MIGRAD',
+ ' FIRST DERIVATIVES OF FCN ARE ALL ZERO')
GO TO 300
ENDIF
C if gdel positive, V not posdef
IF (GDEL .GE. ZERO) THEN
CALL MNWARN('D','MIGRAD',' NEWTON STEP NOT DESCENT.')
IF (NPSDF .EQ. 1) GO TO 1
CALL MNPSDF
NPSDF = 1
GO TO 24
ENDIF
C . . . . do line search
CALL MNLINE(FCN,XXS,FS,STEP,GDEL,TOLER,FUTIL)
IF (AMIN .EQ. FS) GO TO 200
CFROM = 'MIGRAD '
NFCNFR = NFCNMG
CSTATU= 'PROGRESS '
C . get gradient at new point
CALL MNINEX(X)
IF (ISW(3) .EQ. 1) THEN
CALL FCN(NPARX,GIN,FZERO,U,2,FUTIL)
NFCN = NFCN + 1
ENDIF
CALL MNDERI(FCN,FUTIL)
C . calculate new EDM
NPSDF = 0
81 EDM = 0.
GVG = 0.
DELGAM = 0.
GDGSSQ = 0.
DO 100 I= 1, NPAR
RI = 0.
VGI = 0.
DO 90 J= 1, NPAR
M = MAX(I,J)
N = MIN(I,J)
NDEX = M*(M-1)/2 + N
VGI = VGI + VHMAT(NDEX)*(GRD(J)-GS(J))
90 RI = RI + VHMAT(NDEX)* GRD(J)
VG(I) = VGI*0.5
GAMI = GRD(I) - GS(I)
GDGSSQ = GDGSSQ + GAMI**2
GVG = GVG + GAMI*VG(I)
DELGAM = DELGAM + DIRIN(I)*GAMI
100 EDM = EDM + GRD(I)*RI*0.5
EDM = EDM * 0.5 * (1.0 + 3.0*DCOVAR)
C . if EDM negative, not positive-definite
IF (EDM .LT. ZERO .OR. GVG .LE. ZERO) THEN
CALL MNWARN('D','MIGRAD','NOT POS-DEF. EDM OR GVG NEGATIVE.')
CSTATU = 'NOT POSDEF'
IF (NPSDF .EQ. 1) GO TO 230
CALL MNPSDF
NPSDF = 1
GO TO 81
ENDIF
C print information about this iteration
ITER = ITER + 1
IF (ISWTR.GE.3 .OR. (ISWTR.EQ.2.AND.MOD(ITER,10).EQ.1)) THEN
CALL MNWERR
CALL MNPRIN(3,AMIN)
ENDIF
IF (GDGSSQ .EQ. ZERO) CALL MNWARN('D','MIGRAD',
+ 'NO CHANGE IN FIRST DERIVATIVES OVER LAST STEP')
IF (DELGAM .LT. ZERO) CALL MNWARN('D','MIGRAD',
+ 'FIRST DERIVATIVES INCREASING ALONG SEARCH LINE')
C . update covariance matrix
CSTATU = 'IMPROVEMNT'
IF (LDEBUG) WRITE (ISYSWR,'(A,(1X,10G10.3))') ' VHMAT 1 =',
+ (VHMAT(KK),KK=1,10)
DSUM = 0.
VSUM = 0.
DO 120 I=1, NPAR
DO 120 J=1, I
D = DIRIN(I)*DIRIN(J)/DELGAM - VG(I)*VG(J)/GVG
DSUM = DSUM + ABS(D)
NDEX = I*(I-1)/2 + J
VHMAT(NDEX) = VHMAT(NDEX) + 2.0*D
VSUM = VSUM + ABS(VHMAT(NDEX))
120 CONTINUE
C smooth local fluctuations by averaging DCOVAR
DCOVAR = 0.5*(DCOVAR + DSUM/VSUM)
IF (ISWTR.GE.3 .OR. LDEBUG) WRITE (ISYSWR,'(A,F5.1,A)')
+ ' RELATIVE CHANGE IN COV. MATRIX=',DCOVAR*100.,'%'
IF (LDEBUG) WRITE (ISYSWR,'(A,(1X,10G10.3))') ' VHMAT 2 =',
+ (VHMAT(KK),KK=1,10)
IF (DELGAM .LE. GVG) GO TO 135
DO 125 I= 1, NPAR
125 FLNU(I) = DIRIN(I)/DELGAM - VG(I)/GVG
DO 130 I= 1, NPAR
DO 130 J= 1, I
NDEX = I*(I-1)/2 + J
130 VHMAT(NDEX) = VHMAT(NDEX) + 2.0*GVG*FLNU(I)*FLNU(J)
135 CONTINUE
C and see if converged
IF (EDM .LT. 0.1*RHOTOL) GO TO 300
C if not, prepare next iteration
DO 140 I= 1, NPAR
XXS(I) = X(I)
GS(I) = GRD(I)
140 CONTINUE
FS = AMIN
IF (ISW(2) .EQ. 0 .AND. DCOVAR.LT. 0.5 ) ISW(2) = 1
IF (ISW(2) .EQ. 3 .AND. DCOVAR.GT. 0.1 ) ISW(2) = 1
IF (ISW(2) .EQ. 1 .AND. DCOVAR.LT. 0.05) ISW(2) = 3
GO TO 24
C . . . . . end main loop
C . . call limit in MNMIGR
190 ISW(1) = 1
IF (ISW(5) .GE. 0)
+ WRITE (ISYSWR,'(A)') ' CALL LIMIT EXCEEDED IN MIGRAD.'
CSTATU = 'CALL LIMIT'
GO TO 230
C . . fails to improve . .
200 IF (ISWTR .GE. 1) WRITE (ISYSWR,'(A)')
+ ' MIGRAD FAILS TO FIND IMPROVEMENT'
DO 210 I= 1, NPAR
210 X(I) = XXS(I)
IF (EDM .LT. RHOTOL) GO TO 300
IF (EDM .LT. ABS(EPSMA2*AMIN)) THEN
IF (ISWTR .GE. 0) WRITE (ISYSWR, '(A)')
+ ' MACHINE ACCURACY LIMITS FURTHER IMPROVEMENT.'
GO TO 300
ENDIF
IF (ISTRAT .LT. 1) THEN
IF (ISW(5) .GE. 0) WRITE (ISYSWR, '(A)')
+ ' MIGRAD FAILS WITH STRATEGY=0. WILL TRY WITH STRATEGY=1.'
ISTRAT = 1
ENDIF
GO TO 1
C . . fails to converge
230 IF (ISWTR .GE. 0) WRITE (ISYSWR,'(A)')
+ ' MIGRAD TERMINATED WITHOUT CONVERGENCE.'
IF (ISW(2) .EQ. 3) ISW(2) = 1
ISW(4) = -1
GO TO 400
C . . apparent convergence
300 IF (ISWTR .GE. 0) WRITE(ISYSWR,'(/A)')
+ ' MIGRAD MINIMIZATION HAS CONVERGED.'
IF (ITAUR .EQ. 0) THEN
IF (ISTRAT .GE. 2 .OR. (ISTRAT.EQ.1.AND.ISW(2).LT.3)) THEN
IF (ISW(5) .GE. 0) WRITE (ISYSWR, '(/A)')
+ ' MIGRAD WILL VERIFY CONVERGENCE AND ERROR MATRIX.'
CALL MNHESS(FCN,FUTIL)
CALL MNWERR
NPSDF = 0
IF (EDM .GT. RHOTOL) GO TO 10
ENDIF
ENDIF
CSTATU='CONVERGED '
ISW(4) = 1
C come here in any case
400 CONTINUE
CFROM = 'MIGRAD '
NFCNFR = NFCNMG
CALL MNINEX(X)
CALL MNWERR
IF (ISWTR .GE. 0) CALL MNPRIN (3,AMIN)
IF (ISWTR .GE. 1) CALL MNMATU(1)
RETURN
END