source: Sophya/trunk/SophyaExt/CodeMinuit/code/mnmigr.F@ 3848

Last change on this file since 3848 was 2403, checked in by cmv, 22 years ago

Creation du module de code source de MINUIT (CERNLIB) extrait par CMV

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