| [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"
|
|---|
| 12 | CC Calculates the full second-derivative matrix of FCN
|
|---|
| 13 | CC by taking finite differences. When calculating diagonal
|
|---|
| 14 | CC elements, it may iterate so that step size is nearly that
|
|---|
| 15 | CC which gives function change= UP/10. The first derivatives
|
|---|
| 16 | CC of course come as a free side effect, but with a smaller
|
|---|
| 17 | CC step size in order to obtain a known accuracy.
|
|---|
| 18 | CC
|
|---|
| 19 | #include "minuit/d506cm.inc"
|
|---|
| 20 | EXTERNAL FCN,FUTIL
|
|---|
| 21 | DIMENSION YY(MNI)
|
|---|
| 22 | LOGICAL LDEBUG
|
|---|
| 23 | CHARACTER CBF1*22
|
|---|
| 24 | C
|
|---|
| 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
|
|---|
| 46 | C 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 '
|
|---|
| 60 | C . . . . . . diagonal elements .
|
|---|
| 61 | C ISW(2) = 1 if approx, 2 if not posdef, 3 if ok
|
|---|
| 62 | C AIMSAG is the sagitta we are aiming for in second deriv calc.
|
|---|
| 63 | AIMSAG = SQRT(EPSMA2)*(ABS(AMIN)+UP)
|
|---|
| 64 | C Zero the second derivative matrix
|
|---|
| 65 | NPAR2 = NPAR*(NPAR+1)/2
|
|---|
| 66 | DO 10 I= 1,NPAR2
|
|---|
| 67 | 10 VHMAT(I) = 0.
|
|---|
| 68 | C
|
|---|
| 69 | C 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)
|
|---|
| 91 | C
|
|---|
| 92 | C find step which gives sagitta = AIMSAG
|
|---|
| 93 | D = ABS(GSTEP(I))
|
|---|
| 94 | DO 40 ICYC= 1, NCYC
|
|---|
| 95 | C loop here only if SAG=0.
|
|---|
| 96 | DO 25 MULTPY= 1, 5
|
|---|
| 97 | C 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
|
|---|
| 122 | C 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)))
|
|---|
| 133 | C 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
|
|---|
| 137 | C 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
|
|---|
| 143 | C end of step size loop
|
|---|
| 144 | WRITE (CBF1,'(I2,2E10.2)') IEXT,SAG,AIMSAG
|
|---|
| 145 | CALL MNWARN('D','MNHESS','Second Deriv. SAG,AIM= '//CBF1)
|
|---|
| 146 | C
|
|---|
| 147 | 50 CONTINUE
|
|---|
| 148 | NDEX = I*(I+1)/2
|
|---|
| 149 | VHMAT(NDEX) = G2(I)
|
|---|
| 150 | 100 CONTINUE
|
|---|
| 151 | C end of diagonal second derivative loop
|
|---|
| 152 | CALL MNINEX(X)
|
|---|
| 153 | C refine the first derivatives
|
|---|
| 154 | IF (ISTRAT .GT. 0) CALL MNHES1(FCN,FUTIL)
|
|---|
| 155 | ISW(2) = 3
|
|---|
| 156 | DCOVAR = 0.
|
|---|
| 157 | C . . . . 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)
|
|---|
| 176 | C 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
|
|---|
| 189 | C . . . . . . . calculate e d m
|
|---|
| 190 | EDM = 0.
|
|---|
| 191 | DO 230 I= 1, NPAR
|
|---|
| 192 | C 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
|
|---|
| 199 | C 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
|
|---|
| 207 | C 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
|
|---|
![(please configure the [header_logo] section in trac.ini)](/PI/chrome/site/your_project_logo.png){kind=logo}
