339 lines
8.0 KiB
Fortran
339 lines
8.0 KiB
Fortran
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
|
|
|