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