mnmigr.f

Go to the documentation of this file.

*
* $Id: mnmigr.F,v 1.2 1996/03/15 18:02:49 james Exp $
*
* $Log: mnmigr.F,v $
* Revision 1.2  1996/03/15 18:02:49  james
*     Modified Files:
* mnderi.F eliminate possible division by zero
* mnexcm.F suppress print on STOP when print flag=-1
*          set FVAL3 to flag if FCN already called with IFLAG=3
* mninit.F set version 96.03
* mnlims.F remove arguments, not needed
* mnmigr.F VLEN -> LENV in debug print statement
* mnparm.F move call to MNRSET to after NPAR redefined, to zero all
* mnpsdf.F eliminate possible division by zero
* mnscan.F suppress printout when print flag =-1
* mnset.F  remove arguments in call to MNLIMS
* mnsimp.F fix CSTATU so status is PROGRESS only if new minimum
* mnvert.F eliminate possible division by zero
*
* Revision 1.1.1.1  1996/03/07 14:31:30  mclareni
* Minuit
*
*
      SUBROUTINE mnmigr(FCN,FUTIL)
*
* $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        Performs a local function minimization using basically the
CC        method of Davidon-Fletcher-Powell as modified by Fletcher
CC        ref. -- Fletcher, Comp.J. 13,317 (1970)   "switching method"
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 gs(mni), step(mni),  xxs(mni), flnu(mni), vg(mni)
      LOGICAL LDEBUG
      parameter(toler=0.05)
      IF (npar .LE. 0)  RETURN
      IF (amin .EQ. undefi)  CALL mnamin(fcn,futil)
      ldebug = (idbg(4) .GE. 1)
      cfrom = 'MIGRAD  '
      nfcnfr = nfcn
      nfcnmg = nfcn
      cstatu= 'INITIATE  '
      iswtr = isw(5) - 2*itaur
      npfn = nfcn
      nparx = npar
      lenv = npar*(npar+1)/2
      nrstrt = 0
      npsdf = 0
      lined2 = 0
      isw(4) = -1
      rhotol = 1.0e-3*apsi
      IF (iswtr .GE. 1)  WRITE (isyswr,470) istrat,rhotol
  470 FORMAT (' START MIGRAD MINIMIZATION.  STRATEGY',i2,
     +.LT.'.  CONVERGENCE WHEN EDM ',e9.2)
C                                           initialization strategy
      IF (istrat.LT.2 .OR. isw(2).GE.3)  GO TO 2
C                                come (back) here to restart completely
    1 CONTINUE
      IF (nrstrt .GT. istrat)  THEN
         cstatu= 'FAILED    '
         isw(4) = -1
         GO TO 230
         ENDIF
C                                      . get full covariance and gradient
      CALL mnhess(fcn,futil)
      CALL mnwerr
      npsdf = 0
      IF (isw(2) .GE. 1)  GO TO 10
C                                        . get gradient at start point
    2 CONTINUE
      CALL mninex(x)
      IF (isw(3) .EQ. 1) THEN
          CALL fcn(nparx,gin,fzero,u,2,futil)
          nfcn = nfcn + 1
      ENDIF
      CALL mnderi(fcn,futil)
      IF (isw(2) .GE. 1)  GO TO 10
C                                   sometimes start with diagonal matrix
      DO 3 i= 1, npar
         xxs(i) = x(i)
         step(i) = zero
    3 CONTINUE
C                           do line search if second derivative negative
      lined2 = lined2 + 1
      IF (lined2 .LT. (istrat+1)*npar) THEN
      DO 5 i= 1, npar
         IF (g2(i) .GT. zero)  GO TO 5
         step(i) = -sign(gstep(i),grd(i))
         gdel = step(i)*grd(i)
         fs = amin
         CALL mnline(fcn,xxs,fs,step,gdel,toler,futil)
         CALL mnwarn('D','MNMIGR','Negative G2 line search')
         iext = nexofi(i)
         IF (ldebug) WRITE (isyswr,'(A,I3,2G13.3)')
     +    ' Negative G2 line search, param ',iext,fs,amin
         GO TO 2
    5 CONTINUE
      ENDIF
C                           make diagonal error matrix
      DO 8 i=1,npar
         ndex = i*(i-1)/2
           DO 7 j=1,i-1
           ndex = ndex + 1
    7      vhmat(ndex) = 0.
         ndex = ndex + 1
         IF (g2(i) .LE. zero)  g2(i) = 1.
         vhmat(ndex) = 2./g2(i)
    8 CONTINUE
      dcovar = 1.
      IF (ldebug) WRITE (isyswr,'(A,A/(1X,10G10.2))') ' DEBUG MNMIGR,',
     +  ' STARTING MATRIX DIAGONAL,  VHMAT=', (vhmat(kk),kk=1,lenv)
C                                         ready to start first iteration
   10 CONTINUE
      nrstrt = nrstrt + 1
      IF (nrstrt .GT. istrat+1)  THEN
         cstatu= 'FAILED    '
         GO TO 230
         ENDIF
      fs = amin
C                                        . . . get EDM and set up loop
      edm = 0.
         DO 18 i= 1, npar
         gs(i) = grd(i)
         xxs(i) = x(i)
         ndex = i*(i-1)/2
           DO 17 j= 1, i-1
           ndex = ndex + 1
   17      edm = edm + gs(i)*vhmat(ndex)*gs(j)
         ndex = ndex + 1
   18    edm = edm + 0.5 * gs(i)**2 *vhmat(ndex)
      edm = edm * 0.5 * (1.0+3.0*dcovar)
        IF (edm .LT. zero)  THEN
        CALL mnwarn('W','MIGRAD','STARTING MATRIX NOT POS-DEFINITE.')
        isw(2) = 0
        dcovar = 1.
        GO TO 2
        ENDIF
      IF (isw(2) .EQ. 0)  edm=bigedm
      iter = 0
      CALL mninex(x)
      CALL mnwerr
      IF (iswtr .GE. 1)  CALL mnprin(3,amin)
      IF (iswtr .GE. 2)  CALL mnmatu(0)
C                                        . . . . .  start main loop
   24 CONTINUE
      IF (nfcn-npfn .GE. nfcnmx)  GO TO 190
      gdel = 0.
      gssq = 0.
         DO 30  i=1,npar
         ri = 0.
         gssq = gssq + gs(i)**2
           DO 25 j=1,npar
           m = max(i,j)
           n = min(i,j)
           ndex = m*(m-1)/2 + n
   25      ri = ri + vhmat(ndex) *gs(j)
         step(i) = -0.5*ri
   30    gdel = gdel + step(i)*gs(i)
      IF (gssq .EQ. zero)  THEN
          CALL mnwarn('D','MIGRAD',
     +             ' FIRST DERIVATIVES OF FCN ARE ALL ZERO')
          GO TO 300
      ENDIF
C                 if gdel positive, V not posdef
      IF (gdel .GE. zero)  THEN
         CALL mnwarn('D','MIGRAD',' NEWTON STEP NOT DESCENT.')
         IF (npsdf .EQ. 1)  GO TO 1
         CALL mnpsdf
         npsdf = 1
         GO TO 24
         ENDIF
C                                        . . . . do line search
      CALL mnline(fcn,xxs,fs,step,gdel,toler,futil)
      IF (amin .EQ. fs) GO TO 200
      cfrom  = 'MIGRAD  '
      nfcnfr = nfcnmg
      cstatu= 'PROGRESS  '
C                                        . get gradient at new point
      CALL mninex(x)
      IF (isw(3) .EQ. 1) THEN
          CALL fcn(nparx,gin,fzero,u,2,futil)
          nfcn = nfcn + 1
      ENDIF
      CALL mnderi(fcn,futil)
C                                         . calculate new EDM
      npsdf = 0
   81 edm = 0.
      gvg = 0.
      delgam = 0.
      gdgssq = 0.
         DO 100 i= 1, npar
         ri = 0.
         vgi = 0.
           DO 90 j= 1, npar
           m = max(i,j)
           n = min(i,j)
           ndex = m*(m-1)/2 + n
           vgi = vgi + vhmat(ndex)*(grd(j)-gs(j))
   90      ri  =  ri + vhmat(ndex)* grd(j)
      vg(i) = vgi*0.5
      gami = grd(i) - gs(i)
      gdgssq = gdgssq + gami**2
      gvg = gvg + gami*vg(i)
      delgam = delgam + dirin(i)*gami
  100 edm = edm + grd(i)*ri*0.5
      edm = edm * 0.5 * (1.0 + 3.0*dcovar)
C                          . if EDM negative,  not positive-definite
      IF (edm .LT. zero .OR. gvg .LE. zero)  THEN
         CALL mnwarn('D','MIGRAD','NOT POS-DEF. EDM OR GVG NEGATIVE.')
         cstatu = 'NOT POSDEF'
         IF (npsdf .EQ. 1)  GO TO 230
         CALL mnpsdf
         npsdf = 1
         GO TO 81
      ENDIF
C                            print information about this iteration
      iter = iter + 1
      IF (iswtr.GE.3 .OR. (iswtr.EQ.2.AND.mod(iter,10).EQ.1)) THEN
         CALL mnwerr
         CALL mnprin(3,amin)
      ENDIF
      IF (gdgssq .EQ. zero)  CALL mnwarn('D','MIGRAD',
     +           'NO CHANGE IN FIRST DERIVATIVES OVER LAST STEP')
      IF (delgam .LT. zero) CALL mnwarn('D','MIGRAD',
     +          'FIRST DERIVATIVES INCREASING ALONG SEARCH LINE')
C                                        .  update covariance matrix
      cstatu = 'IMPROVEMNT'
        IF (ldebug) WRITE (isyswr,'(A,(1X,10G10.3))') ' VHMAT 1 =',
     +             (vhmat(kk),kk=1,10)
      dsum = 0.
      vsum = 0.
         DO  120  i=1, npar
           DO  120  j=1, i
           d = dirin(i)*dirin(j)/delgam - vg(i)*vg(j)/gvg
           dsum = dsum + abs(d)
           ndex = i*(i-1)/2 + j
           vhmat(ndex) = vhmat(ndex) + 2.0*d
           vsum = vsum + abs(vhmat(ndex))
  120      CONTINUE
C                smooth local fluctuations by averaging DCOVAR
      dcovar = 0.5*(dcovar + dsum/vsum)
      IF (iswtr.GE.3 .OR. ldebug) WRITE (isyswr,'(A,F5.1,A)')
     +      ' RELATIVE CHANGE IN COV. MATRIX=',dcovar*100.,'%'
      IF (ldebug) WRITE (isyswr,'(A,(1X,10G10.3))') ' VHMAT 2 =',
     +             (vhmat(kk),kk=1,10)
      IF (delgam .LE. gvg)  GO TO 135
      DO 125 i= 1, npar
  125 flnu(i) = dirin(i)/delgam - vg(i)/gvg
      DO 130 i= 1, npar
      DO 130 j= 1, i
      ndex = i*(i-1)/2 + j
  130 vhmat(ndex) = vhmat(ndex) + 2.0*gvg*flnu(i)*flnu(j)
  135 CONTINUE
C                                              and see if converged
      IF (edm .LT. 0.1*rhotol)  GO TO 300
C                                    if not, prepare next iteration
      DO 140 i= 1, npar
      xxs(i) = x(i)
      gs(i) = grd(i)
  140 CONTINUE
      fs = amin
      IF (isw(2) .EQ. 0  .AND. dcovar.LT. 0.5 )  isw(2) = 1
      IF (isw(2) .EQ. 3  .AND. dcovar.GT. 0.1 )  isw(2) = 1
      IF (isw(2) .EQ. 1  .AND. dcovar.LT. 0.05)  isw(2) = 3
      GO TO 24
C                                        . . . . .  end main loop
C                                         . . call limit in MNMIGR
  190 isw(1) = 1
      IF (isw(5) .GE. 0)
     +     WRITE (isyswr,'(A)')  ' CALL LIMIT EXCEEDED IN MIGRAD.'
      cstatu = 'CALL LIMIT'
      GO TO 230
C                                         . . fails to improve . .
  200 IF (iswtr .GE. 1)  WRITE (isyswr,'(A)')
     +           ' MIGRAD FAILS TO FIND IMPROVEMENT'
      DO 210 i= 1, npar
  210 x(i) = xxs(i)
      IF (edm .LT. rhotol)  GO TO 300
      IF (edm .LT. abs(epsma2*amin))  THEN
         IF (iswtr .GE. 0)  WRITE (isyswr, '(A)')
     +      ' MACHINE ACCURACY LIMITS FURTHER IMPROVEMENT.'
         GO TO 300
         ENDIF
      IF (istrat .LT. 1)  THEN
         IF (isw(5) .GE. 0) WRITE (isyswr, '(A)')
     +    ' MIGRAD FAILS WITH STRATEGY=0.   WILL TRY WITH STRATEGY=1.'
         istrat = 1
      ENDIF
         GO TO 1
C                                         . . fails to converge
  230 IF (iswtr .GE. 0)  WRITE (isyswr,'(A)')
     +    ' MIGRAD TERMINATED WITHOUT CONVERGENCE.'
      IF (isw(2) .EQ. 3)  isw(2) = 1
      isw(4) = -1
      GO TO 400
C                                         . . apparent convergence
  300 IF (iswtr .GE. 0) WRITE(isyswr,'(/A)')
     +   ' MIGRAD MINIMIZATION HAS CONVERGED.'
      IF (itaur .EQ. 0) THEN
        IF (istrat .GE. 2 .OR. (istrat.EQ.1.AND.isw(2).LT.3)) THEN
           IF (isw(5) .GE. 0)  WRITE (isyswr, '(/A)')
     +      ' MIGRAD WILL VERIFY CONVERGENCE AND ERROR MATRIX.'
           CALL mnhess(fcn,futil)
           CALL mnwerr
           npsdf = 0
           IF (edm .GT. rhotol) GO TO 10
        ENDIF
      ENDIF
      cstatu='CONVERGED '
      isw(4) = 1
C                                           come here in any case
  400 CONTINUE
      cfrom = 'MIGRAD  '
      nfcnfr = nfcnmg
      CALL  mninex(x)
      CALL mnwerr
      IF (iswtr .GE. 0)  CALL mnprin (3,amin)
      IF (iswtr .GE. 1)  CALL mnmatu(1)
      RETURN
      END