SUBROUTINE MIGRAD C -------------------- include 'fit.inc' c Minimization subroutine based on a variable metric method. It is fast c near a minimum or in any 'nearly quadratic' region but slower if the c chi**2 function is badly behaved. It uses the first derivatives of the c chi**2 function, which may either be supplied by the user (subroutine c FNCTN, GG(*)) or estimated by MIGRAD (subroutine DERIVE). REAL*8 VG(MAXPAR), VII(MAXPAR), D, VGI real G(MAXPAR), G2(MAXPAR) real GS(MAXPAR), R(MAXPAR),XXS(MAXPAR), FLNU(MAXPAR) integer iswtr, npfn, iflag, npard, i, ntry, negg2, id, kg, nf, ns integer matgd, j, iter, npargd real parn, rho2, rostop, trace, fs, xtf, fs1, fs2, xbegin, f, ri real gdel, denom, slam, slamin, slamax, tlamin, tlamax, f2 real aa, bb, cc, tlam, f3, gvg, delgam, gami DATA SLAMIN,SLAMAX,TLAMIN,TLAMAX/0.2, 3.0, 0.05, 6.0/ IF (NPAR .LE. 0) RETURN if (isw(2) .eq. 2) isw(2)=1 ISWTR = ISW(5) - ITAUR NPFN = NFCN PARN=NPAR RHO2 = 10.*APSI ! ROSTOP = 1.0E-5 * APSI ROSTOP = 1.0E-3 * APSI TRACE=1. IFLAG=4 IF (ISW(3) .EQ. 1) IFLAG = 2 FS = AMIN IF (ITAUR .LT. 1 .and. isyswr.ne.0) 1 WRITE (ISYSWR,470) ROSTOP, APSI, VTEST GO TO 2 470 FORMAT (/' Start MIGRAD minimization. ' 1/5X,'Convergence criteria: Estimated distance to minimum edm <' 1/,E9.2/5X,'or edm <',E9.2 1,' and fractional change in variance matrix <',E9.2) 1 if (isyswr .ne. 0) WRITE (ISYSWR,520) C. . . . STEP SIZES DIRIN . . . 2 NPARD = NPAR DO 3 I= 1, NPAR D = 0.02* ABS(DIRIN(I)) IF (ISW(2) .GE. 1) D = 0.02* SQRT(ABS(V(I,I))*UP) IF (D .LT. 1.0E-5 *ABS(X(I))) D = 1.0E-5 * X(I) ! was 1E-6 (M.Z. 1.9.95) 3 DIRIN(I) = D C. . . . . . STARTING GRADIENT NTRY = 0 4 NEGG2 = 0 DO 10 ID= 1, NPARD I = ID + NPAR - NPARD D = DIRIN(I) XTF = X(I) X(I) = XTF + D CALL FNCTN(X,FS1) NFCN = NFCN + 1 X(I) = XTF - D CALL FNCTN(X,FS2) NFCN = NFCN + 1 X(I) = XTF GS(I) = (FS1-FS2)/(2.0 * D) G2(I) = (FS1 + FS2 - 2.0*AMIN) / D**2 IF (G2(I) .GT. 1.0E-30) GO TO 10 C . . . SEARCH IF G2 .LE. 0. . . if (isyswr .ne. 0) WRITE (ISYSWR,520) NEGG2 = NEGG2 + 1 NTRY = NTRY + 1 IF (NTRY .GT. 4) GO TO 230 D = 50.*ABS(DIRIN(I)) XBEGIN = XTF IF (GS(I) .LT. 0.) DIRIN(I) = -DIRIN(I) KG = 0 NF = 0 NS = 0 5 X(I) = XTF + D CALL FNCTN(X,F) NFCN = NFCN + 1 IF (F .LE. AMIN) GO TO 6 C FAILURE IF (KG .EQ. 1) GO TO 8 KG = -1 NF = NF + 1 D = -0.4*D IF (NF .LT. 10) GO TO 5 D = 1000.*D GO TO 7 C SUCCESS 6 XTF = X(I) D = 3.0*D AMIN = F KG = 1 NS = NS + 1 IF (NS .LT. 10) GO TO 5 IF (AMIN .LT. FS) GO TO 8 D = 0.001*D 7 XTF = XBEGIN G2(I) = 1.0 NEGG2 = NEGG2 - 1 8 X(I) = XTF DIRIN(I) = 0.1*D FS = AMIN 10 CONTINUE IF (NEGG2 .GE. 1) GO TO 4 NTRY = 0 MATGD = 1 C . . . . . . DIAGONAL MATRIX IF (ISW(2) .GT. 1) GO TO 15 11 NTRY = 1 MATGD = 0 DO 13 I= 1, NPAR DO 12 J= 1, NPAR 12 V(I,J) = 0. 13 V(I,I) = 2.0/G2(I) C. . . GET SIGMA AND SET UP LOOP 15 SIGMA = 0. DO 18 I= 1, NPAR IF (V(I,I) .LE. 0.) GO TO 11 RI = 0. DO 17 J= 1, NPAR XXS(I) = X(I) 17 RI= RI+ V(I,J) * GS(J) 18 SIGMA = SIGMA + GS(I) *RI *0.5 IF (SIGMA .GE. 0.) GO TO 20 if (isyswr .ne. 0) WRITE (ISYSWR,520) IF (NTRY.EQ.0) GO TO 11 ISW(2) = 0 GO TO 230 20 ISW(2) = 1 ITER = 0 CALL INTOEX(X) IF (ISWTR .GE. 1) CALL FIT_PRINT(0) C . . . . . START MAIN LOOP 24 CONTINUE GDEL = 0. DO 30 I=1,NPAR RI = 0. DO 25 J=1,NPAR 25 RI = RI + V(I,J) *GS(J) DIRIN(I) = -0.5*RI GDEL = GDEL + DIRIN(I)*GS(I) C.LINEAR SEARCH ALONG -VG . . . 30 X(I) =XXS(I) + DIRIN(I) CALL FNCTN(X, F) NFCN=NFCN+1 C. QUADR INTERP USING SLOPE GDEL DENOM = 2.0*(F-AMIN-GDEL) IF (DENOM .LE. 0.) GO TO 35 SLAM = -GDEL/DENOM IF (SLAM .GT. SLAMAX) GO TO 35 IF (SLAM .LT. SLAMIN) SLAM=SLAMIN GO TO 40 35 SLAM = SLAMAX 40 IF (ABS(SLAM-1.0) .LT. 0.1) GO TO 70 DO 45 I= 1, NPAR 45 X(I) =XXS(I) + SLAM*DIRIN(I) CALL FNCTN(X,F2) NFCN = NFCN + 1 C. QUADR INTERP USING 3 POINTS AA = FS/SLAM BB = F/(1.0-SLAM) CC = F2/ (SLAM*(SLAM-1.0)) DENOM = 2.0*(AA+BB+CC) IF (DENOM .LE. 0.) GO TO 48 TLAM = (AA*(SLAM+1.0) + BB*SLAM + CC)/DENOM IF (TLAM .GT. TLAMAX) GO TO 48 IF (TLAM .LT. TLAMIN) TLAM=TLAMIN GO TO 50 48 TLAM = TLAMAX 50 CONTINUE DO 51 I= 1, NPAR 51 X(I) = XXS(I)+TLAM*DIRIN(I) CALL FNCTN(X,F3) NFCN = NFCN + 1 IF (F.GE.AMIN .AND. F2.GE.AMIN .AND. F3.GE.AMIN) GO TO 200 IF (F .LT. F2 .AND. F .LT. F3) GO TO 61 IF (F2 .LT. F3) GO TO 58 55 F = F3 SLAM = TLAM GO TO 65 58 F = F2 GO TO 65 61 SLAM = 1.0 65 DO 67 I= 1, NPAR DIRIN(I) = DIRIN(I)*SLAM 67 X(I) = XXS(I) + DIRIN(I) 70 AMIN = F ISW(2) = 2 IF (SIGMA+FS-AMIN .LT. ROSTOP) GO TO 170 IF (SIGMA+RHO2+FS-AMIN .GT. APSI) GO TO 75 IF (TRACE .LT. VTEST) GO TO 170 75 CONTINUE IF (NFCN-NPFN .GE. NFCNMX) GO TO 190 ITER = ITER + 1 IF (ISWTR.GE. 3 .OR.(ISWTR.EQ. 2 .AND. MOD(ITER,10) .EQ.1)) 1 CALL FIT_PRINT(0) C. . . GET GRADIENT AND SIGMA . IF (ISW(3) .NE. 1) GO TO 80 CALL FNCTN(X,AMIN) NFCN = NFCN + 1 80 CALL DERIVE(G,G2) RHO2 = SIGMA SIGMA = 0. GVG = 0. DELGAM = 0. DO 100 I= 1, NPAR RI = 0. VGI = 0. DO 90 J= 1, NPAR VGI = VGI + V(I,J)*(G(J)-GS(J)) 90 RI = RI + V(I,J) *G (J) R(I) = RI * 0.5 VG(I) = VGI*0.5 GAMI = G(I) - GS(I) GVG = GVG + GAMI*VG(I) DELGAM = DELGAM + DIRIN(I)*GAMI 100 SIGMA = SIGMA + G(I)*R(I) IF (SIGMA .LT. 0.) GO TO 1 IF (GVG .LE. 0.) GO TO 105 IF (DELGAM .LE. 0.) GO TO 105 GO TO 107 105 IF (SIGMA .LT. 0.1*ROSTOP) GO TO 170 GO TO 1 107 CONTINUE C. UPDATE COVARIANCE MATRIX TRACE=0. DO 120 I= 1, NPAR VII(I) = V(I,I) DO 120 J=1,NPAR D = DIRIN(I)*DIRIN(J)/DELGAM - VG(I)*VG(J)/GVG 120 V(I,J) = V(I,J) + 2.0*D 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, NPAR 130 V(I,J) = V(I,J) + 2.0*GVG*FLNU(I)*FLNU(J) 135 CONTINUE DO 140 I= 1, NPAR 140 TRACE = TRACE + ((V(I,I)-VII(I))/(V(I,I)+VII(I)))**2 TRACE = SQRT(TRACE/PARN) CALL UCOPY(X,XXS,NPAR) CALL UCOPY(G,GS,NPAR) FS = F GO TO 24 C . . . . . END MAIN LOOP 170 if (isyswr .ne. 0) WRITE(ISYSWR,500) 500 FORMAT (/' MIGRAD minimization has converged') c do i=1,npar c type '(X,8F10.5)',(v(i,j)/sqrt(v(i,i)*v(j,j)),j=1,npar) c enddo ISW(2) = 3 IF(ISWTR .GE. 0) CALL FIT_PRINT(1-ITAUR) IF (ITAUR .GT. 0) GO TO 435 IF (MATGD .GT. 0) GO TO 435 NPARGD = NPAR*(NPAR+5)/2 IF (NFCN-NPFN .GE. NPARGD) GO TO 435 if (isyswr .ne. 0) 1 WRITE (ISYSWR,'(X,A)') 'Covariance matrix inaccurate' IF (ISW(2) .GE. 2) ISW(2) = 3 GO TO 435 190 ISW(1) = 1 if (nfcnmx .gt. 0) GO TO 230 if (isyswr .ne. 0) write(isyswr, '(/a)') ' MIGRAD aborted' goto 231 200 if (isyswr .ne.0) WRITE (ISYSWR,650) 650 FORMAT (/' MIGRAD fails to find improvement') CALL UCOPY(XXS,X,NPAR) ISW(2) = 1 IF (SIGMA .LT. ROSTOP) GO TO 170 IF (MATGD .GT. 0) GO TO 2 if (isw(2) .eq. 1) isw(2)=2 230 if (isyswr .ne. 0) WRITE (ISYSWR,510) 510 FORMAT (/' MIGRAD terminated without convergence') 231 CALL FNCTN(X,AMIN) CALL FIT_PRINT(1-ITAUR) 435 RETURN 520 FORMAT (/' Covariance matrix is not positive-definite') END SUBROUTINE DERIVE(GG,GG2) C ---------------------------- include 'fit.inc' real GG(*),GG2(*) integer iflag, i, lc real eps, xtf, fs1, fs2, dd IF (ISW(3) .EQ. 1) GO TO 100 IFLAG = 4 DO 46 I=1,NPAR EPS = 0.1 * ABS(DIRIN(I)) IF (ISW(2) .GE. 1) EPS = EPS + 0.005*SQRT(V(I,I)*UP) IF (EPS .LT. 1.0E-8*ABS(X(I))) EPS = 1.0E-8*X(I) XTF = X(I) X(I) = XTF + EPS CALL FNCTN(X,FS1) NFCN=NFCN+1 X(I) = XTF - EPS CALL FNCTN(X,FS2) NFCN=NFCN+1 C. . . . . . . FIRST DERIVATIVE GG(I)= (FS1-FS2)/(2.0*EPS) C. . . ERROR ON FIRST DERIVATIVE GG2(I)= (FS1+FS2-2.0*AMIN)/(2.0*EPS) X(I) = XTF 46 CONTINUE CALL INTOEX(X) GO TO 200 C. DERIVATIVES CALC BY FCN 100 DO 150 I= 1, NU LC=LCORSP(I) IF (LC .LT. 1) GO TO 150 IF (LCODE(I) .GT. 1) GO TO 120 GG(LC)=GG(I) GO TO 150 120 DD = (BLIM(I)-ALIM(I))*0.5 *COS(X(LC)) GG(LC)=GG(I)*DD 150 CONTINUE 200 RETURN END SUBROUTINE UCOPY(FROM,TO,N) C ----------------------------- integer n, i real FROM(N),TO(N) IF(N.LT.1)RETURN DO 100 I=1,N 100 TO(I)=FROM(I) RETURN END