*
* $Id: mnmigr.f,v 1.1.1.1 2003/11/25 22:45:45 dorman Exp $
*
* $Log: mnmigr.f,v $
* Revision 1.1.1.1  2003/11/25 22:45:45  dorman
* Xspec12.0 Export Version 11-25-2003
*
* Revision 7.0  2002/07/12 16:03:24  dorman
*
* version number change
*
* Revision 1.1  2002/06/26 14:52:12  dorman
*
*
* minuit library
*
* Revision 1.1  2001/12/28 02:53:55  kaa
* Added MINUIT source code.
*
* 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)
      INCLUDE "d506dp.inc"
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
      INCLUDE "d506cm.inc"
      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,
     +'.  CONVERGENCE WHEN EDM .LT.',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
