C PROGRAM AMOEBA IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (MAXVAR=47,MAXVEC=48,NVAR=2) CHARACTER*64 NA,STATUS DIMENSION D1(MAXVAR,MAXVEC),D0(MAXVAR),D2(MAXVAR),VALUE(MAXVEC) DATA ITYPE/2/,TOL/1.E-7/,ALPHA/1.0/,GAMMA/2.0/ C C PERFORM A AMOEBA MINIMIZATION OF THE FUNCTION FTBM(D2) C MAXVAR IS THE NUMBER OF CONSTANTS IN FTBM C NVAR IS THE NUMBER OF CONSTANTS TO BE MINIMIZED C NVEC IS THE NUMBER OF INITIAL CONSTANT SETS TO BE READ IN C ITYPE DETERMINES HOW THE INITIAL VECTORS ARE TO BE READ IN C ITMAX IS THE MAXIMUM NUMBER OF REFLECTIONS C WE REQUIRE THAT NVEC=NVAR+1 C PRINT*,' ENTER THE FILE NAME FOR THE STARTING CONSTANTS' READ(*,'(A)')NA STATUS='OLD' CALL MAOPEN(11,NA,STATUS) PRINT*,' ENTER A 1 IF THERE IS ONLY ONE VECTOR OF CONSTANTS' PRINT*,' ENTER A 2 IF THIS IS A CONTINUING RUN AND THERE ARE' PRINT*,' NVAR+1, VECTORS OF CONSTANTS IN THE FILE ',NA(1:12) PRINT*,' AND ALSO ACCURATE VALUES FOR FTBM FOR EACH VECTOR' PRINT*,' ENTER A 0 IF THERE ARE NVAR+1 VECTORS AND VALUES, BUT' PRINT*,' THE VALUES MUST BE RECALCULATED' READ(*,*)ITYPE PRINT*,' ENTER THE OUTPUT FILE NAME' READ(*,'(A)')NA STATUS='UNKNOWN' CALL MAOPEN(12,NA,STATUS) PRINT*,' ENTER THE MAXIMUM NUMBER OF REFLECTIONS' READ(*,*)ITMAX NVEC=NVAR+1 SF=1.*NVEC ITER=0 NREFL=0 C DO 20 I=1,NVEC IF(ITYPE.EQ.0) THEN READ(11,*) (D2(J),J=1,NVAR),VALUE(I) VALUE(I)=FTBM(D2) ENDIF IF(ITYPE.EQ.2) READ(11,*) (D2(J),J=1,NVAR),VALUE(I) IF(ITYPE.EQ.1) THEN IF(I.EQ.1) THEN READ(11,*) (D2(J),J=1,NVAR) ELSE SIGN=-1 DO 25 J=1,NVAR D2(J)=D1(J,1) 25 CONTINUE D2(I-1)=D2(I-1)+0.02*DMAX1(1D0,ABS(D2(I-1)))*SIGN SIGN=-SIGN ENDIF VALUE(I)=FTBM(D2) ENDIF DO 30 J=1,NVAR D1(J,I)=D2(J) 30 CONTINUE WRITE(*,1001) VALUE(I) 1001 FORMAT(' INITIAL',E20.7) WRITE(*,1002) (D2(J),J=1,NVAR) 1002 FORMAT(4E20.7) 20 CONTINUE C C SORT THE VECTORS C IMAX POINTS TO THE VECTOR WITH THE LARGEST VALUE C I2MAX POINTS TO THE VECTOR WITH THE SECOND LARGEST VALUE C IMIN POINTS TO THE VECTOR WITH THE SMALLEST VALUE C 40 CONTINUE IMIN=1 IF(VALUE(1).GT.VALUE(2)) THEN IMAX=1 I2MAX=2 ELSE IMAX=2 I2MAX=1 ENDIF DO 50 I=1,NVEC IF(VALUE(I).LT.VALUE(IMIN)) IMIN=I IF(VALUE(I).GT.VALUE(IMAX)) THEN I2MAX=IMAX IMAX=I ELSE IF(VALUE(I).GT.VALUE(I2MAX)) THEN IF(I.NE.IMAX) I2MAX=I ENDIF 50 CONTINUE C C CHECK IF DONE C TEMP=(VALUE(IMAX)-VALUE(IMIN))/VALUE(IMIN) IF(TEMP.LT.TOL) GOTO 500 IF(ITER.EQ.ITMAX) THEN WRITE(*,*) ' MAXIMUM NUMBER OF ITERATIONS REACHED' GO TO 500 ENDIF C CONVERGENCE IS REALLY SLOW SO SHRINK IF(TEMP.LT.10*TOL.AND.NREFL.GT.5*NVEC) GOTO 400 ITER=ITER+1 C C REFLECTION C NREFL=NREFL+1 DO 100 J=1,NVAR SUM=0.0 DO 110 I=1,NVEC IF(I.EQ.IMAX) GOTO 110 SUM=SUM+D1(J,I) 110 CONTINUE D0(J)=SUM/(NVEC-1) 100 CONTINUE DO 120 J=1,NVAR D2(J)=(1.+ALPHA)*D0(J)-ALPHA*D1(J,IMAX) 120 CONTINUE CVAL=FTBM(D2) WRITE(*,1003) CVAL 1003 FORMAT(' REFLECTION',E20.7) WRITE(*,1002) (D2(J),J=1,NVAR) C VALUE IS WORSE SO DO CONTRACTION IF(CVAL.GE.VALUE(I2MAX)) GOTO 200 C VALUE IS BETTER VALUE(IMAX)=CVAL DO 140 J=1,NVAR D1(J,IMAX)=D2(J) 140 CONTINUE C VALUE IS MUCH BETTER SO DO EXPANSION IF(CVAL.LT.VALUE(IMIN)) GOTO 300 GOTO 40 C C CONTRACTION C 200 CONTINUE BETA=0.75 DO 210 ITEMP=1,3 IF(CVAL.LE.VALUE(IMAX)) THEN VALUE(IMAX)=CVAL DO 220 J=1,NVAR D1(J,IMAX)=D2(J) 220 CONTINUE ENDIF DO 230 J=1,NVAR D2(J)=BETA*D1(J,IMAX)+(1.-BETA)*D0(J) 230 CONTINUE CVAL=FTBM(D2) WRITE(*,1004) BETA,CVAL 1004 FORMAT(' CONTRACTION',F7.3,E20.7) WRITE(*,1002) (D2(J),J=1,NVAR) C VALUE IS BETTER IF(CVAL.LT.VALUE(I2MAX)) THEN VALUE(IMAX)=CVAL DO 250 J=1,NVAR D1(J,IMAX)=D2(J) 250 CONTINUE IF(CVAL.LT.VALUE(IMIN)) WRITE(*,1005) CVAL GOTO 40 ENDIF BETA=BETA-0.25 210 CONTINUE C VALUE IS WORSE SO SHRINK IT GOTO 400 C C EXPANSION C 300 CONTINUE WRITE(*,1005) CVAL 1005 FORMAT(' NEW MIN',E15.6) DO 310 J=1,NVAR D2(J)=GAMMA*D2(J)+(1.-GAMMA)*D0(J) 310 CONTINUE CCVAL=FTBM(D2) WRITE(*,1006) CCVAL 1006 FORMAT(' EXPANSION',E20.7) WRITE(*,1002) (D2(J),J=1,NVAR) C VALUE IS WORSE IF(CCVAL.GT.CVAL) GOTO 40 C VALUE IS BETTER SO USE IT RATHER THAN THE REFLECTED POINT VALUE(IMAX)=CCVAL DO 320 J=1,NVAR D1(J,IMAX)=D2(J) 320 CONTINUE WRITE(*,1005) CCVAL GOTO 40 C C SHRINKAGE C 400 CONTINUE NREFL=0 C THE CENTER OF THE SHRINK IS TAKEN TO BE AT THE AVERAGE OF C THE SQUARES OF THE DISTANCES FROM THE MIN POINT DO 410 J=1,NVAR SUM=0.0 DO 420 I=1,NVEC IF(I.EQ.IMIN) GOTO 420 SUM=SUM+(D1(J,I)-D1(J,IMIN))**2 420 CONTINUE D0(J)=SQRT(SUM/(NVEC-1)) 410 CONTINUE SF=-.75*SF WRITE(*,1002) SF C FORM THE NEW SET OF VECTORS DO 430 I=1,NVEC IF(I.EQ.IMIN) GOTO 430 DO 440 J=1,NVAR D2(J)=D1(J,IMIN) D1(J,I)=D2(J) 440 CONTINUE IF(I.LT.IMIN) THEN ITEMP=I ELSE ITEMP=I-1 ENDIF D2(ITEMP)=SF*D0(ITEMP)+D1(ITEMP,IMIN) D1(ITEMP,I)=D2(ITEMP) VALUE(I)=FTBM(D2) WRITE(*,1007) VALUE(I) 1007 FORMAT(' SHRINKAGE',E20.7) WRITE(*,1002) (D2(J),J=1,NVAR) 430 CONTINUE GOTO 40 C C ALL DONE C 500 CONTINUE DO 510 I=1,NVEC WRITE(12,1002) (D1(J,I),J=1,NVAR),VALUE(I) 510 CONTINUE STOP END SUBROUTINE MAOPEN(L,NA,STATUS) CHARACTER*64 NA,STATUS 5 CONTINUE IF(STATUS.EQ.'BIG')THEN OPEN(L,FILE=NA,STATUS='UNKNOWN',ERR=1197) RETURN ENDIF OPEN(L,FILE=NA,STATUS=STATUS,BLANK='ZERO',ERR=1197) RETURN 1197 PRINT*,' COULD NOT OPEN FILE= ',NA PRINT*,' ENTER A NEW NAME OR STOP' READ(*,'(A)')NA IF(NA.EQ.'STOP'.OR.NA.EQ.'stop')STOP GOTO 5 END SUBROUTINE MAOPUF(L,NA,STATUS) CHARACTER*64 NA,STATUS,NAT NAT='(V:512)'//NA 5 OPEN(L,FILE=NAT,STATUS=STATUS,ERR=1197,FORM='UNFORMATTED') RETURN 1197 PRINT*,' COULD NOT OPEN UNFORMATTED FILE= ',NAT PRINT*,' ENTER A NEW NAME OR STOP' READ(*,'(A)')NAT IF(NAT.EQ.'STOP'.OR.NAT.EQ.'stop')STOP GOTO 5 END SUBROUTINE ERRSET(I1,I2,I3,I4) RETURN END C$INCLUDE FTBM