source: Sophya/trunk/SophyaExt/CodeMinuit/code/mnhess.F@ 2403

Last change on this file since 2403 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: 6.9 KB
RevLine 
[2403]1*
2* $Id: mnhess.F,v 1.1.1.1 2003-06-11 14:18:28 cmv Exp $
3*
4* $Log: not supported by cvs2svn $
5* Revision 1.1.1.1 1996/03/07 14:31:30 mclareni
6* Minuit
7*
8*
9#include "minuit/pilot.h"
10 SUBROUTINE MNHESS(FCN,FUTIL)
11#include "minuit/d506dp.inc"
12CC Calculates the full second-derivative matrix of FCN
13CC by taking finite differences. When calculating diagonal
14CC elements, it may iterate so that step size is nearly that
15CC which gives function change= UP/10. The first derivatives
16CC of course come as a free side effect, but with a smaller
17CC step size in order to obtain a known accuracy.
18CC
19#include "minuit/d506cm.inc"
20 EXTERNAL FCN,FUTIL
21 DIMENSION YY(MNI)
22 LOGICAL LDEBUG
23 CHARACTER CBF1*22
24C
25 LDEBUG = (IDBG(3) .GE. 1)
26 IF (AMIN .EQ. UNDEFI) CALL MNAMIN(FCN,FUTIL)
27 IF (ISTRAT .LE. 0) THEN
28 NCYC = 3
29 TLRSTP = 0.5
30 TLRG2 = 0.1
31 ELSE IF (ISTRAT .EQ. 1) THEN
32 NCYC = 5
33 TLRSTP = 0.3
34 TLRG2 = 0.05
35 ELSE
36 NCYC = 7
37 TLRSTP = 0.1
38 TLRG2 = 0.02
39 ENDIF
40 IF (ISW(5).GE.2 .OR. LDEBUG) WRITE (ISYSWR,'(A)')
41 + ' START COVARIANCE MATRIX CALCULATION.'
42 CFROM = 'HESSE '
43 NFCNFR = NFCN
44 CSTATU= 'OK '
45 NPARD = NPAR
46C make sure starting at the right place
47 CALL MNINEX(X)
48 NPARX = NPAR
49 CALL FCN(NPARX,GIN,FS1,U,4,FUTIL)
50 NFCN = NFCN + 1
51 IF (FS1 .NE. AMIN) THEN
52 DF = AMIN - FS1
53 WRITE (CBF1(1:12),'(G12.3)') DF
54 CALL MNWARN('D','MNHESS',
55 + 'function value differs from AMIN by '//CBF1(1:12) )
56 ENDIF
57 AMIN = FS1
58 IF (LDEBUG) WRITE (ISYSWR,'(A,A)') ' PAR D GSTEP ',
59 +' D G2 GRD SAG '
60C . . . . . . diagonal elements .
61C ISW(2) = 1 if approx, 2 if not posdef, 3 if ok
62C AIMSAG is the sagitta we are aiming for in second deriv calc.
63 AIMSAG = SQRT(EPSMA2)*(ABS(AMIN)+UP)
64C Zero the second derivative matrix
65 NPAR2 = NPAR*(NPAR+1)/2
66 DO 10 I= 1,NPAR2
67 10 VHMAT(I) = 0.
68C
69C Loop over variable parameters for second derivatives
70 IDRV = 2
71 DO 100 ID= 1, NPARD
72 I = ID + NPAR - NPARD
73 IEXT = NEXOFI(I)
74 IF (G2(I) .EQ. ZERO) THEN
75 WRITE (CBF1(1:4),'(I4)') IEXT
76 CALL MNWARN('W','HESSE',
77 + 'Second derivative enters zero, param '//CBF1(1:4) )
78 WINT = WERR(I)
79 IF (NVARL(IEXT) .GT. 1) THEN
80 CALL MNDXDI(X(I),I,DXDI)
81 IF (ABS(DXDI) .LT. .001) THEN
82 WINT = .01
83 ELSE
84 WINT = WINT/ABS(DXDI)
85 ENDIF
86 ENDIF
87 G2(I) = UP/WINT**2
88 ENDIF
89 XTF = X(I)
90 DMIN = 8.*EPSMA2*ABS(XTF)
91C
92C find step which gives sagitta = AIMSAG
93 D = ABS(GSTEP(I))
94 DO 40 ICYC= 1, NCYC
95C loop here only if SAG=0.
96 DO 25 MULTPY= 1, 5
97C take two steps
98 X(I) = XTF + D
99 CALL MNINEX(X)
100 NPARX = NPAR
101 CALL FCN(NPARX,GIN,FS1,U,4,FUTIL)
102 NFCN = NFCN + 1
103 X(I) = XTF - D
104 CALL MNINEX(X)
105 CALL FCN(NPARX,GIN,FS2,U,4,FUTIL)
106 NFCN = NFCN + 1
107 X(I) = XTF
108 SAG = 0.5*(FS1+FS2-2.0*AMIN)
109 IF (SAG .NE. ZERO) GO TO 30
110 IF (GSTEP(I) .LT. ZERO) THEN
111 IF (D .GE. .5) GO TO 26
112 D = 10.*D
113 IF (D .GT. 0.5) D = 0.51
114 GO TO 25
115 ENDIF
116 D = 10.*D
117 25 CONTINUE
118 26 WRITE (CBF1(1:4),'(I4)') IEXT
119 CALL MNWARN('W','HESSE',
120 + 'Second derivative zero for parameter'//CBF1(1:4) )
121 GO TO 390
122C SAG is not zero
123 30 G2BFOR = G2(I)
124 G2(I) = 2.*SAG/D**2
125 GRD(I) = (FS1-FS2)/(2.*D)
126 IF (LDEBUG) WRITE (ISYSWR,31) I,IDRV,GSTEP(I),D,G2(I),GRD(I),SAG
127 31 FORMAT (I4,I2,6G12.5)
128 GSTEP(I) = SIGN(D,GSTEP(I))
129 DIRIN(I) = D
130 YY(I) = FS1
131 DLAST = D
132 D = SQRT(2.0*AIMSAG/ABS(G2(I)))
133C if parameter has limits, max int step size = 0.5
134 STPINM = 0.5
135 IF (GSTEP(I) .LT. ZERO) D = MIN(D,STPINM)
136 IF (D .LT. DMIN) D = DMIN
137C see if converged
138 IF (ABS((D-DLAST)/D) .LT. TLRSTP) GO TO 50
139 IF (ABS((G2(I)-G2BFOR)/G2(I)) .LT. TLRG2 ) GO TO 50
140 D = MIN(D, 10.*DLAST)
141 D = MAX(D, 0.1*DLAST)
142 40 CONTINUE
143C end of step size loop
144 WRITE (CBF1,'(I2,2E10.2)') IEXT,SAG,AIMSAG
145 CALL MNWARN('D','MNHESS','Second Deriv. SAG,AIM= '//CBF1)
146C
147 50 CONTINUE
148 NDEX = I*(I+1)/2
149 VHMAT(NDEX) = G2(I)
150 100 CONTINUE
151C end of diagonal second derivative loop
152 CALL MNINEX(X)
153C refine the first derivatives
154 IF (ISTRAT .GT. 0) CALL MNHES1(FCN,FUTIL)
155 ISW(2) = 3
156 DCOVAR = 0.
157C . . . . off-diagonal elements
158 IF (NPAR .EQ. 1) GO TO 214
159 DO 200 I= 1, NPAR
160 DO 180 J= 1, I-1
161 XTI = X(I)
162 XTJ = X(J)
163 X(I) = XTI + DIRIN(I)
164 X(J) = XTJ + DIRIN(J)
165 CALL MNINEX(X)
166 CALL FCN(NPARX,GIN,FS1,U,4,FUTIL)
167 NFCN = NFCN + 1
168 X(I) = XTI
169 X(J) = XTJ
170 ELEM = (FS1+AMIN-YY(I)-YY(J)) / (DIRIN(I)*DIRIN(J))
171 NDEX = I*(I-1)/2 + J
172 VHMAT(NDEX) = ELEM
173 180 CONTINUE
174 200 CONTINUE
175 214 CALL MNINEX(X)
176C verify matrix positive-definite
177 CALL MNPSDF
178 DO 220 I= 1, NPAR
179 DO 219 J= 1, I
180 NDEX = I*(I-1)/2 + J
181 P(I,J) = VHMAT(NDEX)
182 219 P(J,I) = P(I,J)
183 220 CONTINUE
184 CALL MNVERT(P,MAXINT,MAXINT,NPAR,IFAIL)
185 IF (IFAIL .GT. 0) THEN
186 CALL MNWARN('W','HESSE', 'Matrix inversion fails.')
187 GO TO 390
188 ENDIF
189C . . . . . . . calculate e d m
190 EDM = 0.
191 DO 230 I= 1, NPAR
192C off-diagonal elements
193 NDEX = I*(I-1)/2
194 DO 225 J= 1, I-1
195 NDEX = NDEX + 1
196 ZTEMP = 2.0 * P(I,J)
197 EDM = EDM + GRD(I)*ZTEMP*GRD(J)
198 225 VHMAT(NDEX) = ZTEMP
199C diagonal elements
200 NDEX = NDEX + 1
201 VHMAT(NDEX) = 2.0 * P(I,I)
202 EDM = EDM + P(I,I) * GRD(I)**2
203 230 CONTINUE
204 IF (ISW(5).GE.1 .AND. ISW(2).EQ.3 .AND. ITAUR.EQ.0)
205 + WRITE(ISYSWR,'(A)')' COVARIANCE MATRIX CALCULATED SUCCESSFULLY'
206 GO TO 900
207C failure to invert 2nd deriv matrix
208 390 ISW(2) = 1
209 DCOVAR = 1.
210 CSTATU = 'FAILED '
211 IF (ISW(5) .GE. 0) WRITE (ISYSWR,'(A)')
212 + ' MNHESS FAILS AND WILL RETURN DIAGONAL MATRIX. '
213 DO 395 I= 1, NPAR
214 NDEX = I*(I-1)/2
215 DO 394 J= 1, I-1
216 NDEX = NDEX + 1
217 394 VHMAT(NDEX) = 0.0
218 NDEX = NDEX +1
219 G2I = G2(I)
220 IF (G2I .LE. ZERO) G2I = 1.0
221 395 VHMAT(NDEX) = 2.0/G2I
222 900 RETURN
223 END
Note: See TracBrowser for help on using the repository browser.