mnhess.f File Reference

Go to the source code of this file.

Functions/Subroutines
subroutine	mnhess (FCN, FUTIL)

Function/Subroutine Documentation

mnhess()

subroutine mnhess	(	external	FCN,
		external	FUTIL
	)

Definition at line 10 of file mnhess.f.

References a, b, c, d, e, f, g, h, softenpik::half, i, j, m, maxp, mnamin(), mndxdi(), mnhes1(), mninex(), mnpsdf(), mnvert(), mnwarn(), n, o, p, parameter(), up(), x, xs, and z.

Referenced by mncntr(), mncuve(), mnexcm(), and mnmigr().

*
* $Id: d506dp.inc,v 1.1.1.1 1996/03/07 14:31:32 mclareni Exp $
*
* $Log: d506dp.inc,v $
* Revision 1.1.1.1  1996/03/07 14:31:32  mclareni
* Minuit
*
*
*
*
* d506dp.inc
*
C ************ DOUBLE PRECISION VERSION *************
      IMPLICIT DOUBLE PRECISION (a-h,o-z)
CC        Calculates the full second-derivative matrix of FCN
CC        by taking finite differences. When calculating diagonal
CC        elements, it may iterate so that step size is nearly that
CC        which gives function change= UP/10. The first derivatives
CC        of course come as a free side effect, but with a smaller
CC        step size in order to obtain a known accuracy.
CC
*
* $Id: d506cm.inc,v 1.1.1.1 1996/03/07 14:31:32 mclareni Exp $
*
* $Log: d506cm.inc,v $
* Revision 1.1.1.1  1996/03/07 14:31:32  mclareni
* Minuit
*
*
*
*
* d506cm.inc
*
      parameter(mne=100 , mni=50)
      parameter(mnihl=mni*(mni+1)/2)
      CHARACTER*10 cpnam
      COMMON
     1/mn7nam/ cpnam(mne)
     2/mn7ext/ u(mne)     ,alim(mne)  ,blim(mne)
     3/mn7err/ erp(mni)   ,ern(mni)   ,werr(mni)  ,globcc(mni)
     4/mn7inx/ nvarl(mne) ,niofex(mne),nexofi(mni)
     5/mn7int/ x(mni)     ,xt(mni)    ,dirin(mni)
     6/mn7fx2/ xs(mni)    ,xts(mni)   ,dirins(mni)
     7/mn7der/ grd(mni)   ,g2(mni)    ,gstep(mni) ,gin(mne) ,dgrd(mni)
     8/mn7fx3/ grds(mni)  ,g2s(mni)   ,gsteps(mni)
     9/mn7fx1/ ipfix(mni) ,npfix
     a/mn7var/ vhmat(mnihl)
     b/mn7vat/ vthmat(mnihl)
     c/mn7sim/ p(mni,mni+1),pstar(mni),pstst(mni) ,pbar(mni),prho(mni)
C
      parameter(maxdbg=10, maxstk=10, maxcwd=20, maxp=30, maxcpt=101)
      parameter(zero=0.0,  one=1.0,   half=0.5)
      COMMON
     d/mn7npr/ maxint ,npar   ,maxext ,nu
     e/mn7iou/ isysrd ,isyswr ,isyssa ,npagwd ,npagln ,newpag
     e/mn7io2/ istkrd(maxstk) ,nstkrd ,istkwr(maxstk) ,nstkwr
     f/mn7tit/ cfrom  ,cstatu ,ctitl  ,cword  ,cundef ,cvrsn ,covmes
     g/mn7flg/ isw(7) ,idbg(0:maxdbg) ,nblock ,icomnd
     h/mn7min/ amin   ,up     ,edm    ,fval3  ,epsi   ,apsi  ,dcovar
     i/mn7cnv/ nfcn   ,nfcnmx ,nfcnlc ,nfcnfr ,itaur,istrat,nwrmes(2)
     j/mn7arg/ word7(maxp)
     k/mn7log/ lwarn  ,lrepor ,limset ,lnolim ,lnewmn ,lphead
     l/mn7cns/ epsmac ,epsma2 ,vlimlo ,vlimhi ,undefi ,bigedm,updflt
     m/mn7rpt/ xpt(maxcpt)    ,ypt(maxcpt)
     n/mn7cpt/ chpt(maxcpt)
     o/mn7xcr/ xmidcr ,ymidcr ,xdircr ,ydircr ,ke1cr  ,ke2cr
      CHARACTER ctitl*50, cword*(maxcwd), cundef*10, cfrom*8,
     +          cvrsn*6,  covmes(0:3)*22, cstatu*10, chpt*1
      LOGICAL   lwarn, lrepor, limset, lnolim, lnewmn, lphead
      EXTERNAL fcn,futil
      dimension yy(mni)
      LOGICAL ldebug
      CHARACTER cbf1*22
C
      ldebug = (idbg(3) .GE. 1)
      IF (amin .EQ. undefi)  CALL mnamin(fcn,futil)
      IF (istrat .LE. 0) THEN
         ncyc = 3
         tlrstp = 0.5
         tlrg2  = 0.1
      ELSE IF (istrat .EQ. 1) THEN
         ncyc = 5
         tlrstp = 0.3
         tlrg2  = 0.05
      ELSE
         ncyc = 7
         tlrstp = 0.1
         tlrg2  = 0.02
      ENDIF
      IF (isw(5).GE.2 .OR. ldebug)  WRITE (isyswr,'(A)')
     +   '   START COVARIANCE MATRIX CALCULATION.'
      cfrom = 'HESSE   '
      nfcnfr = nfcn
      cstatu= 'OK        '
      npard = npar
C                 make sure starting at the right place
      CALL mninex(x)
      nparx = npar
      CALL fcn(nparx,gin,fs1,u,4,futil)
      nfcn = nfcn + 1
      IF (fs1 .NE. amin) THEN
         df = amin - fs1
         WRITE (cbf1(1:12),'(G12.3)') df
         CALL mnwarn('D','MNHESS',
     +       'function value differs from AMIN by '//cbf1(1:12) )
      ENDIF
      amin = fs1
      IF (ldebug) WRITE (isyswr,'(A,A)') ' PAR D   GSTEP          ',
     +' D          G2         GRD         SAG    '
C                                        . . . . . . diagonal elements .
C         ISW(2) = 1 if approx, 2 if not posdef, 3 if ok
C         AIMSAG is the sagitta we are aiming for in second deriv calc.
      aimsag = sqrt(epsma2)*(abs(amin)+up)
C         Zero the second derivative matrix
      npar2 = npar*(npar+1)/2
      DO 10 i= 1,npar2
   10 vhmat(i) = 0.
C
C         Loop over variable parameters for second derivatives
      idrv = 2
      DO 100 id= 1, npard
      i = id + npar - npard
      iext = nexofi(i)
      IF (g2(i) .EQ. zero) THEN
           WRITE (cbf1(1:4),'(I4)') iext
           CALL mnwarn('W','HESSE',
     +      'Second derivative enters zero, param '//cbf1(1:4) )
        wint = werr(i)
        IF (nvarl(iext) .GT. 1) THEN
           CALL mndxdi(x(i),i,dxdi)
           IF (abs(dxdi) .LT. .001) THEN
              wint = .01
           ELSE
              wint = wint/abs(dxdi)
           ENDIF
        ENDIF
        g2(i) = up/wint**2
      ENDIF
      xtf = x(i)
      dmin = 8.*epsma2*abs(xtf)
C
C                               find step which gives sagitta = AIMSAG
      d = abs(gstep(i))
      DO 40 icyc= 1, ncyc
C                               loop here only if SAG=0.
      DO 25 multpy= 1, 5
C           take two steps
         x(i) = xtf + d
         CALL mninex(x)
         nparx = npar
         CALL fcn(nparx,gin,fs1,u,4,futil)
         nfcn = nfcn + 1
         x(i) = xtf - d
         CALL mninex(x)
         CALL fcn(nparx,gin,fs2,u,4,futil)
         nfcn = nfcn + 1
         x(i) = xtf
         sag = 0.5*(fs1+fs2-2.0*amin)
         IF (sag .NE. zero) GO TO 30
         IF (gstep(i) .LT. zero) THEN
           IF (d .GE. .5)  GO TO 26
           d = 10.*d
           IF (d .GT. 0.5)  d = 0.51
           GO TO 25
         ENDIF
         d = 10.*d
   25 CONTINUE
   26      WRITE (cbf1(1:4),'(I4)') iext
           CALL mnwarn('W','HESSE',
     +      'Second derivative zero for parameter'//cbf1(1:4) )
           GO TO 390
C                             SAG is not zero
   30 g2bfor = g2(i)
      g2(i) = 2.*sag/d**2
      grd(i) = (fs1-fs2)/(2.*d)
      IF (ldebug) WRITE (isyswr,31) i,idrv,gstep(i),d,g2(i),grd(i),sag
   31 FORMAT (i4,i2,6g12.5)
      gstep(i) = sign(d,gstep(i))
      dirin(i) = d
      yy(i) = fs1
      dlast = d
      d = sqrt(2.0*aimsag/abs(g2(i)))
C         if parameter has limits, max int step size = 0.5
      stpinm = 0.5
      IF (gstep(i) .LT. zero)  d = min(d,stpinm)
      IF (d .LT. dmin)  d = dmin
C           see if converged
      IF (abs((d-dlast)/d)          .LT. tlrstp)  GO TO 50
      IF (abs((g2(i)-g2bfor)/g2(i)) .LT. tlrg2 )  GO TO 50
      d = min(d, 10.*dlast)
      d = max(d, 0.1*dlast)
   40 CONTINUE
C                       end of step size loop
      WRITE (cbf1,'(I2,2E10.2)') iext,sag,aimsag
      CALL mnwarn('D','MNHESS','Second Deriv. SAG,AIM= '//cbf1)
C
   50 CONTINUE
      ndex = i*(i+1)/2
      vhmat(ndex) = g2(i)
  100 CONTINUE
C                              end of diagonal second derivative loop
      CALL mninex(x)
C                                     refine the first derivatives
      IF (istrat .GT. 0) CALL mnhes1(fcn,futil)
      isw(2) = 3
      dcovar = 0.
C                                        . . . .  off-diagonal elements
      IF (npar .EQ. 1)  GO TO 214
      DO 200 i= 1, npar
      DO 180 j= 1, i-1
      xti = x(i)
      xtj = x(j)
      x(i) = xti + dirin(i)
      x(j) = xtj + dirin(j)
      CALL mninex(x)
      CALL fcn(nparx,gin,fs1,u,4,futil)
      nfcn = nfcn + 1
      x(i) = xti
      x(j) = xtj
      elem = (fs1+amin-yy(i)-yy(j)) / (dirin(i)*dirin(j))
      ndex = i*(i-1)/2 + j
      vhmat(ndex) = elem
  180 CONTINUE
  200 CONTINUE
  214 CALL mninex(x)
C                  verify matrix positive-definite
      CALL mnpsdf
      DO 220 i= 1, npar
      DO 219 j= 1, i
      ndex = i*(i-1)/2 + j
      p(i,j) = vhmat(ndex)
  219 p(j,i) = p(i,j)
  220 CONTINUE
      CALL mnvert(p,maxint,maxint,npar,ifail)
      IF (ifail .GT. 0)  THEN
        CALL mnwarn('W','HESSE', 'Matrix inversion fails.')
        GO TO 390
      ENDIF
C                                        . . . . . . .  calculate  e d m
      edm = 0.
        DO 230 i= 1, npar
C                              off-diagonal elements
        ndex = i*(i-1)/2
          DO 225 j= 1, i-1
          ndex = ndex + 1
          ztemp = 2.0 * p(i,j)
          edm = edm + grd(i)*ztemp*grd(j)
  225     vhmat(ndex) = ztemp
C                              diagonal elements
        ndex = ndex + 1
        vhmat(ndex) = 2.0 * p(i,i)
        edm = edm  + p(i,i) * grd(i)**2
  230   CONTINUE
      IF (isw(5).GE.1 .AND. isw(2).EQ.3 .AND. itaur.EQ.0)
     + WRITE(isyswr,'(A)')' COVARIANCE MATRIX CALCULATED SUCCESSFULLY'
      GO TO 900
C                              failure to invert 2nd deriv matrix
  390 isw(2) = 1
      dcovar = 1.
      cstatu = 'FAILED    '
      IF (isw(5) .GE. 0) WRITE (isyswr,'(A)')
     +        '  MNHESS FAILS AND WILL RETURN DIAGONAL MATRIX. '
      DO 395 i= 1, npar
      ndex = i*(i-1)/2
      DO 394 j= 1, i-1
      ndex = ndex + 1
  394 vhmat(ndex) = 0.0
      ndex = ndex +1
      g2i = g2(i)
      IF (g2i .LE. zero)  g2i = 1.0
  395 vhmat(ndex) = 2.0/g2i
  900 RETURN

Here is the call graph for this function:

Here is the caller graph for this function: