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

Last change on this file 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
Line 
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.