Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

mnmigr.F@ 3891

Visit:

Last change on this file since 3891 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"
27	CC Performs a local function minimization using basically the
28	CC method of Davidon-Fletcher-Powell as modified by Fletcher
29	CC ref. -- Fletcher, Comp.J. 13,317 (1970) "switching method"
30	CC
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)
55	C initialization strategy
56	IF (ISTRAT.LT.2 .OR. ISW(2).GE.3) GO TO 2
57	C 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
64	C . 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
69	C . 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
78	C sometimes start with diagonal matrix
79	DO 3 I= 1, NPAR
80	XXS(I) = X(I)
81	STEP(I) = ZERO
82	3 CONTINUE
83	C 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
99	C 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)
112	C 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
120	C . . . 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)
144	C . . . . . 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
164	C 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
172	C . . . . 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 '
178	C . 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)
185	C . 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)RI0.5
206	EDM = EDM * 0.5 * (1.0 + 3.0*DCOVAR)
207	C . 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
216	C 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')
226	C . 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
240	C 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.0GVGFLNU(I)*FLNU(J)
253	135 CONTINUE
254	C and see if converged
255	IF (EDM .LT. 0.1*RHOTOL) GO TO 300
256	C 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
266	C . . . . . end main loop
267	C . . 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
273	C . . 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
290	C . . 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
296	C . . 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
311	C 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.

Download in other formats:

Original Format