234 lines
5.9 KiB
Fortran
234 lines
5.9 KiB
Fortran
SUBROUTINE SIMPLEX
|
|
C ------------------
|
|
|
|
c minimization subroutine using a simple method by nelder and mead *).
|
|
c it is reasonably fast when far from minimum, and may also be used to
|
|
c converge to the exact minimum, but is mostly slower than "migrad".
|
|
c
|
|
c *) - nelder, j.a., mead, r. (1965), comput. j. 7, 308
|
|
c - james, f., roos, m. (1985), cern documentation d506
|
|
|
|
include 'fit.inc'
|
|
|
|
integer jh,jl
|
|
real p(maxpar,maxpar+1),y(maxpar+1)
|
|
common/comsim/jh,jl,y,p
|
|
|
|
integer iflag, i, kg, ns, nf, k, ignal, ncycl, j, npfn, nparp1
|
|
integer jhold
|
|
real rho1, rho2, alpha, beta, gamma, rhomin, rhomax, wg
|
|
real ynpp1, absmin, aming, bestx, f, sig2, pb, ystar, ystst
|
|
real y1, y2, rho, yrho, ypbar
|
|
|
|
REAL PSTAR(MAXPAR) ,PSTST(MAXPAR) ,PBAR(MAXPAR) ,PRHO(MAXPAR)
|
|
DATA ALPHA,BETA,GAMMA,RHOMIN,RHOMAX / 1.0, 0.5, 2.0, 4.0, 8.0/
|
|
|
|
IF (NPAR .LE. 0) RETURN
|
|
NPFN=NFCN
|
|
NPARP1=NPAR+1
|
|
RHO1 = 1.0 + ALPHA
|
|
RHO2 = RHO1 + ALPHA*GAMMA
|
|
WG = 1.0/FLOAT(NPAR)
|
|
IFLAG=4
|
|
if (isyswr .ne. 0) WRITE(ISYSWR,100) EPSI
|
|
100 FORMAT(/' Start SIMPLEX minimization'
|
|
1/4X,'Convergence criterion: estimated distance to minimum edm <'
|
|
1,E10.2)
|
|
|
|
DO 2 I= 1, NPAR
|
|
IF (ISW(2) .GE. 1) DIRIN(I) = SQRT(V(I,I)*UP)
|
|
IF (ABS(DIRIN(I)) .LT. 1.0E-5*ABS(X(I))) DIRIN(I)=1.0E-5*X(I) ! was 1e-8 M.Z. 6.02
|
|
IF(ITAUR.LT. 1) V(I,I) = DIRIN(I)**2/UP
|
|
2 CONTINUE
|
|
IF (ITAUR .LT. 1) ISW(2) = 1
|
|
C** CHOOSE THE INITIAL SIMPLEX USING SINGLE-PARAMETER SEARCHES
|
|
1 CONTINUE
|
|
YNPP1 = AMIN
|
|
JL = NPARP1
|
|
Y(NPARP1) = AMIN
|
|
ABSMIN = AMIN
|
|
DO 10 I= 1, NPAR
|
|
AMING = AMIN
|
|
PBAR(I) = X(I)
|
|
BESTX = X(I)
|
|
KG = 0
|
|
NS = 0
|
|
NF = 0
|
|
4 X(I) = BESTX + DIRIN(I)
|
|
CALL FNCTN(X,F)
|
|
NFCN = NFCN + 1
|
|
IF (F .LE. AMING) GO TO 6
|
|
C FAILURE
|
|
IF (KG .EQ. 1) GO TO 8
|
|
KG = -1
|
|
NF = NF + 1
|
|
DIRIN(I) = DIRIN(I) * (-0.4)
|
|
IF (NF .LT. 3) GO TO 4
|
|
NS = 6
|
|
C SUCCESS
|
|
6 BESTX = X(I)
|
|
DIRIN(I) = DIRIN(I) * 3.0
|
|
AMING = F
|
|
KG = 1
|
|
NS = NS + 1
|
|
IF (NS .LT. 6) GO TO 4
|
|
C LOCAL MINIMUM FOUND IN ITH DIRECTION
|
|
8 Y(I) = AMING
|
|
IF (AMING .LT. ABSMIN) JL = I
|
|
IF (AMING .LT. ABSMIN) ABSMIN = AMING
|
|
X(I) = BESTX
|
|
DO 9 K= 1, NPAR
|
|
9 P(K,I) = X(K)
|
|
10 CONTINUE
|
|
JH = NPARP1
|
|
AMIN=Y(JL)
|
|
CALL RAZZIA(YNPP1,PBAR)
|
|
DO 20 I= 1, NPAR
|
|
20 X(I) = P(I,JL)
|
|
CALL INTOEX(X)
|
|
IF (ISW(5) .GE. 1) CALL FIT_PRINT(0)
|
|
SIGMA = SIGMA * 10.
|
|
SIG2 = SIGMA
|
|
IGNAL = 0
|
|
NCYCL=0
|
|
C. . . . . START MAIN LOOP
|
|
50 CONTINUE
|
|
IF (IGNAL .GE. 10) GO TO 1
|
|
IF (SIG2 .LT. EPSI .AND. SIGMA.LT.EPSI) GO TO 76
|
|
SIG2 = SIGMA
|
|
IF ((NFCN-NPFN) .GT. NFCNMX) GO TO 78
|
|
C CALCULATE NEW POINT * BY REFLECTION
|
|
DO 60 I= 1, NPAR
|
|
PB = 0.
|
|
DO 59 J= 1, NPARP1
|
|
59 PB = PB + WG * P(I,J)
|
|
PBAR(I) = PB - WG * P(I,JH)
|
|
60 PSTAR(I)=(1.+ALPHA)*PBAR(I)-ALPHA*P(I,JH)
|
|
CALL FNCTN(PSTAR,YSTAR)
|
|
NFCN=NFCN+1
|
|
IF(YSTAR.GE.AMIN) GO TO 70
|
|
C POINT * BETTER THAN JL, CALCULATE NEW POINT **
|
|
DO 61 I=1,NPAR
|
|
61 PSTST(I)=GAMMA*PSTAR(I)+(1.-GAMMA)*PBAR(I)
|
|
CALL FNCTN(PSTST,YSTST)
|
|
NFCN=NFCN+1
|
|
C TRY A PARABOLA THROUGH PH, PSTAR, PSTST. MIN = PRHO
|
|
Y1 = (YSTAR-Y(JH)) * RHO2
|
|
Y2 = (YSTST-Y(JH)) * RHO1
|
|
RHO = 0.5 * (RHO2*Y1 -RHO1*Y2) / (Y1 -Y2)
|
|
IF (RHO .LT. RHOMIN) GO TO 66
|
|
IF (RHO .GT. RHOMAX) RHO = RHOMAX
|
|
DO 64 I= 1, NPAR
|
|
64 PRHO(I) = RHO*PSTAR(I) + (1. -RHO)*P(I,JH)
|
|
CALL FNCTN(PRHO,YRHO)
|
|
NFCN = NFCN + 1
|
|
IF (YRHO .LT. Y(JL) .AND. YRHO .LT. YSTST) GO TO 65
|
|
IF (YSTST .LT. Y(JL)) GO TO 67
|
|
IF (YRHO .GT. Y(JL)) GO TO 66
|
|
C ACCEPT MINIMUM POINT OF PARABOLA, PRHO
|
|
65 CALL RAZZIA (YRHO,PRHO)
|
|
IGNAL = MAX0(IGNAL-2, 0)
|
|
GO TO 68
|
|
66 IF (YSTST .LT. Y(JL)) GO TO 67
|
|
IGNAL = MAX0(IGNAL-1, 0)
|
|
CALL RAZZIA(YSTAR,PSTAR)
|
|
GO TO 68
|
|
67 IGNAL = MAX0(IGNAL-2, 0)
|
|
675 CALL RAZZIA(YSTST,PSTST)
|
|
68 NCYCL=NCYCL+1
|
|
IF (ISW(5) .LT. 2) GO TO 50
|
|
IF (ISW(5) .GE. 3 .OR. MOD(NCYCL, 10) .EQ. 0) CALL FIT_PRINT(0)
|
|
GO TO 50
|
|
C POINT * IS NOT AS GOOD AS JL
|
|
70 IF (YSTAR .GE. Y(JH)) GO TO 73
|
|
JHOLD = JH
|
|
CALL RAZZIA(YSTAR,PSTAR)
|
|
IF (JHOLD .NE. JH) GO TO 50
|
|
C CALCULATE NEW POINT **
|
|
73 DO 74 I=1,NPAR
|
|
74 PSTST(I)=BETA*P(I,JH)+(1.-BETA)*PBAR(I)
|
|
CALL FNCTN (PSTST,YSTST)
|
|
NFCN=NFCN+1
|
|
IF(YSTST.GT.Y(JH)) GO TO 1
|
|
C POINT ** IS BETTER THAN JH
|
|
IF (YSTST .LT. AMIN) GO TO 675
|
|
IGNAL = IGNAL + 1
|
|
CALL RAZZIA(YSTST,PSTST)
|
|
GO TO 50
|
|
C. . . . . . END MAIN LOOP
|
|
76 if (isyswr .ne. 0) WRITE(ISYSWR,120)
|
|
120 format(/' SIMPLEX minimization has converged')
|
|
GO TO 80
|
|
78 if (nfcnmx .eq. 0) then
|
|
if (isyswr .ne. 0) write(isyswr,'(/a)') ' SIMPLEX aborted'
|
|
else
|
|
if (isyswr .ne. 0) WRITE(ISYSWR,130)
|
|
130 format(/' SIMPLEX terminates without convergence')
|
|
endif
|
|
ISW(1) = 1
|
|
80 DO 82 I=1,NPAR
|
|
PB = 0.
|
|
DO 81 J=1,NPARP1
|
|
81 PB = PB + WG * P(I,J)
|
|
82 PBAR(I) = PB - WG * P(I,JH)
|
|
CALL FNCTN(PBAR,YPBAR)
|
|
NFCN=NFCN+1
|
|
IF (YPBAR .LT. AMIN) CALL RAZZIA(YPBAR,PBAR)
|
|
IF (NFCNMX+NPFN-NFCN .LT. 3*NPAR) GO TO 90
|
|
IF (SIGMA .GT. 2.0*EPSI) GO TO 1
|
|
90 CALL FNCTN(X,AMIN)
|
|
NFCN = NFCN + 1
|
|
CALL FIT_PRINT(1-ITAUR)
|
|
RETURN
|
|
END
|
|
|
|
SUBROUTINE RAZZIA(YNEW,PNEW)
|
|
C ------------------------------
|
|
|
|
include 'fit.inc'
|
|
|
|
integer jh,jl
|
|
real p(maxpar,maxpar+1),y(maxpar+1)
|
|
common/comsim/jh,jl,y,p
|
|
|
|
integer i, nparp1, j
|
|
real ynew, us, pbig, plit
|
|
real PNEW(MAXPAR)
|
|
|
|
DO 10 I=1,NPAR
|
|
10 P(I,JH)=PNEW(I)
|
|
Y(JH)=YNEW
|
|
IF(YNEW.GE.AMIN) GO TO 18
|
|
DO 15 I=1,NPAR
|
|
15 X(I)=PNEW(I)
|
|
CALL INTOEX(X)
|
|
AMIN=YNEW
|
|
JL=JH
|
|
18 CONTINUE
|
|
JH=1
|
|
NPARP1=NPAR+1
|
|
20 DO 25 J=2,NPARP1
|
|
IF (Y(J) .GT. Y(JH)) JH = J
|
|
25 CONTINUE
|
|
SIGMA = Y(JH) - Y(JL)
|
|
IF (SIGMA .LE. 0.) GO TO 45
|
|
US = 1.0/SIGMA
|
|
DO 35 I= 1, NPAR
|
|
PBIG = P(I,1)
|
|
PLIT = PBIG
|
|
DO 30 J= 2, NPARP1
|
|
IF (P(I,J) .GT. PBIG) PBIG = P(I,J)
|
|
IF (P(I,J) .LT. PLIT) PLIT = P(I,J)
|
|
30 CONTINUE
|
|
DIRIN(I) = PBIG - PLIT
|
|
IF (ITAUR .LT. 1 ) V(I,I) = 0.5*(V(I,I) +US*DIRIN(I)**2)
|
|
35 CONTINUE
|
|
40 RETURN
|
|
45 WRITE (ISYSWR, 1000) NPAR
|
|
GO TO 40
|
|
1000 FORMAT (
|
|
1 ' ***** FUNCTION VALUE DOES NOT SEEM TO DEPEND ON ANY OF THE '
|
|
1, I3,' VARIABLE PARAMETERS'/15X,'VERIFY THAT STEP SIZES ARE'
|
|
1,' BIG ENOUGH AND CHECK FUN LOGIC.'/1X, 79('*')/1X,79('H')///)
|
|
END
|