C ====================================================================== SUBROUTINE CREATE(M,S,V,N,XL,XR,FVAL,WA) C ====================================================================== * CREATES A POPULATION OF M MEMBERS. * EACH MEMBER IS AN ARRAY WITH N COMPONENTS (I.E. A POINT IN R^N) * XL(I), XR(I) ARE THE LOWER AND UPPER BOUNDS FOR THE I-TH COMPONENT. * S(N,M) CONTAINS THE POPULATION, I.E.: M, N-VECTORS. * V(N,M) CONTAINS THE VELOCITY COMPONENTS FOR EACH MEMBER. * FVAL(J) CONTAINS THE VALUE OF THE OBJECTIVE FUNCTION EVALUATED AT THE * POINT CORRESPONDING TO J-TH COLUMN OF S. * WA IS A WORK ARRAY. ** INPUT: M,N,XL,XR ** OUTPUT: S,FVAL IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION S(N,M),V(N,M) DIMENSION XL(N),XR(N),WA(N),FVAL(M) DO 1 I=1,M DO 2 J=1,N S(J,I) = (XR(J)-XL(J))*RANM()+XL(J) WA(J) = S(J,I) V(J,I) = 2*RANM()-1.D0 2 CONTINUE FVAL(I) = FUNMIN(WA,N) 1 CONTINUE END C ====================================================================== SUBROUTINE EQUATE(POP,N,M,S) C ====================================================================== IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION POP(N,M), S(N,M) DO 1 I = 1,N DO 2 J = 1,M POP(I,J) = S(I,J) 2 CONTINUE 1 CONTINUE END C ====================================================================== SUBROUTINE EVALUATE(POP,N,M,FVAL,WA) C ====================================================================== IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION POP(N,M),FVAL(M),WA(N) DO 1 I=1,M DO 2 J=1,N WA(J) = POP(J,I) 2 CONTINUE FVAL(I) = FUNMIN(WA,N) 1 CONTINUE END C ====================================================================== SUBROUTINE PICKG(POP,N,M,ICL,P) C ====================================================================== * SELECTS THE ICL COLUMN OF POP, AND STORES IT TO P IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION POP(N,M),P(N) DO 1 I=1,N P(I) = POP(I,ICL) 1 CONTINUE END C ====================================================================== SUBROUTINE MINEL(FVAL,M,IEL) c ====================================================================== * INPUT : M, FVAL(M) * OUTPUT: IEL, THE POSITION OF THE SMALLEST VALUE IN FVAL(M) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION FVAL(M) FMIN = FVAL(1) IEL = 1 DO 1 I = 1,M IF (FVAL(I) .LE. FMIN) THEN FMIN = FVAL(I) IEL = I ENDIF 1 CONTINUE END C ====================================================================== SUBROUTINE ADDTWO(POP,XL,XR,N,M,VEL) C ====================================================================== * UPDATES THE BIRD POSITIONS AND CHECKS THAT BIRDS REMAIN IN BOUNDS. * INPUT : N, M, POP(N,M), VEL(N,M), XL(N), XR(N) * OUTPUT: POP(N,M) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION POP(N,M),VEL(N,M),XL(N),XR(N) DO 1 J=1,M DO 2 I=1,N POP(I,J) = POP(I,J) + VEL(I,J) IF (POP(I,J).LT.XL(I)) POP(I,J) = XL(I) IF (POP(I,J).GT.XR(I)) POP(I,J) = XR(I) 2 CONTINUE 1 CONTINUE END C ====================================================================== SUBROUTINE UPDBES(POP,N,M,BESTP,FVAL,BESTV) C ====================================================================== * KEEPS THE HISTORICALY BEST POSITION OF EVERY BIRD. * INPUT : N, M, POP(N,M), BESTP(N,M), FVAL(M), BESTV(M) * OUTPUT: BESTP(N,M), BESTV(M) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION POP(N,M),BESTP(N,M) DIMENSION FVAL(M), BESTV(M) DO 1 J=1,M IF (FVAL(J).LT.BESTV(J)) THEN BESTV(J) = FVAL(J) DO 2 I=1,N BESTP(I,J) = POP(I,J) 2 CONTINUE ENDIF 1 CONTINUE END C ====================================================================== SUBROUTINE UPDVEL(POP,N,M,BESTP,XB,VEL,C1,C2) C ====================================================================== * UPDATES THE BIRD VELOCITIES. * INPUT : N, M, POP(N,M), BESTP(N,M), VEL(N,M), XB(N) * OUTPUT: VEL(N,M) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION POP(N,M),BESTP(N,M),VEL(N,M) DIMENSION XB(N) PHI = C1 + C2 CHI = 2.D0/ABS(2-PHI-SQRT(PHI**2-4*PHI)) DO 1 J=1,M * if you move the loop below (DO 2) up here, you pick different * random numbers for each component XI1 = RANM() XI2 = RANM() DO 2 I = 1,N VEL(I,J) = CHI*(VEL(I,J) + C1*XI1*(BESTP(I,J)-POP(I,J)) + > C2*XI2*(XB(I) - POP(I,J))) 2 CONTINUE 1 CONTINUE END C ====================================================================== PROGRAM SMEENOS C ====================================================================== IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (N = 2, M = 150, NITER = 20) PARAMETER (CFL = 2.D0, CFG = 5.D0) DIMENSION POP(N,M),S(N,M),FVAL(M),WA(N),XB(N) DIMENSION VEL(N,M),V(N,M) DIMENSION XL(N),XR(N) DIMENSION INDX(M),SCR(M) DIMENSION BESTP(N,M),BESTV(M) CHARACTER CFO*50 DATA XL/(N)*-10.D0/ DATA XR/(N)*10.D0/ C CFO='(T10,I7,T20,D14.7,T40,D14.7,10000(/T40,D14.7))' WRITE(CFO(29:33),'(I5)')N-1 WRITE(*,'(57(''=''))') WRITE(*,*) 'SMEENOS . . . ''FLIES'' WITH THE FOLLOWING DATA !!! :' WRITE(*,*) 'NUMBER OF VARIABLES: ',N WRITE(*,*) 'SWARM SIZE: ',M WRITE(*,*) 'MAX.# OF MOVES: ',NITER WRITE(*,*) 'LOCAL FACTOR: ',CFL WRITE(*,*) 'GLOBAL FACTOR: ',CFG WRITE(*,'(57(''=''))') C CREATE INITIAL BIRD STATE. (POSITIONS & VELOCITIES). CALL CREATE(M,POP,VEL,N,XL,XR,FVAL,WA) C INITIALIZE BEST POSITIONS AND VALUES CALL EQUATE(BESTP,N,M,POP) CALL EQUATE(BESTV,1,M,FVAL) MOVE = 1 30 CONTINUE C CHECK FOR TERMINATION. IF (MOVE .EQ. NITER) GO TO 31 C PREPARE NEXT MOVE MOVE = MOVE + 1 C UPDATE POSITIONS. CALL ADDTWO(POP,XL,XR,N,M,VEL) C CALCULATE THE (OBJECTIVE) VALUES AT THE NEW POSITIONS. CALL EVALUATE(POP,N,M,FVAL,WA) C UPDATE BEST POSITIONS. CALL UPDBES(POP,N,M,BESTP,FVAL,BESTV) C DETERMINE THE CURRENTLY GLOBALY BEST POSITION AND ASSIGN IT TO XB CALL MINEL(BESTV,M,IGLO) CALL PICKG(BESTP,N,M,IGLO,XB) C1 = CFL C2 = CFG C UPDATE THE VELOCITIES. CALL UPDVEL(POP,N,M,BESTP,XB,VEL,C1,C2) C START A NEW ITERATION (BIRD MOVE) GO TO 30 31 CONTINUE WRITE(*,88)'PSO:','CALLS','MINIMUM','POINT' 88 FORMAT(1X,A8,T13,A5,T19,A13,T32,A17) WRITE(*,CFO) M*NITER,BESTV(IGLO),(BESTP(J,IGLO),J=1,N) WRITE(*,'(57(''-''))') C END C --------------------------------------------------------------------- SUBROUTINE RANLUX(RVEC,LENV) C --------------------------------------------------------------------- C Subtract-and-borrow random number generator proposed by C Marsaglia and Zaman, implemented by F. James with the name C RCARRY in 1991, and later improved by Martin Luescher C in 1993 to produce "Luxury Pseudorandom Numbers". C Fortran 77 coded by F. James, 1993 C C references: C M. Luscher, Computer Physics Communications 79 (1994) 100 C F. James, Computer Physics Communications 79 (1994) 111 C C LUXURY LEVELS. C ------ ------ The available luxury levels are: C C level 0 (p=24): equivalent to the original RCARRY of Marsaglia C and Zaman, very long period, but fails many tests. C level 1 (p=48): considerable improvement in quality over level 0, C now passes the gap test, but still fails spectral test. C level 2 (p=97): passes all known tests, but theoretically still C defective. C level 3 (p=223): DEFAULT VALUE. Any theoretically possible C correlations have very small chance of being observed. C level 4 (p=389): highest possible luxury, all 24 bits chaotic. C C!!! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C!!! Calling sequences for RANLUX: ++ C!!! CALL RANLUX (RVEC, LEN) returns a vector RVEC of LEN ++ C!!! 32-bit random floating point numbers between ++ C!!! zero (not included) and one (also not incl.). ++ C!!! CALL RLUXGO(LUX,INT,K1,K2) initializes the generator from ++ C!!! one 32-bit integer INT and sets Luxury Level LUX ++ C!!! which is integer between zero and MAXLEV, or if ++ C!!! LUX .GT. 24, it sets p=LUX directly. K1 and K2 ++ C!!! should be set to zero unless restarting at a break++ C!!! point given by output of RLUXAT (see RLUXAT). ++ C!!! CALL RLUXAT(LUX,INT,K1,K2) gets the values of four integers++ C!!! which can be used to restart the RANLUX generator ++ C!!! at the current point by calling RLUXGO. K1 and K2++ C!!! specify how many numbers were generated since the ++ C!!! initialization with LUX and INT. The restarting ++ C!!! skips over K1+K2*E9 numbers, so it can be long.++ C!!! A more efficient but less convenient way of restarting is by: ++ C!!! CALL RLUXIN(ISVEC) restarts the generator from vector ++ C!!! ISVEC of 25 32-bit integers (see RLUXUT) ++ C!!! CALL RLUXUT(ISVEC) outputs the current values of the 25 ++ C!!! 32-bit integer seeds, to be used for restarting ++ C!!! ISVEC must be dimensioned 25 in the calling program ++ C!!! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION RVEC(LENV) DIMENSION SEEDS(24), ISEEDS(24), ISDEXT(25) PARAMETER (MAXLEV=4, LXDFLT=3) DIMENSION NDSKIP(0:MAXLEV) DIMENSION NEXT(24) PARAMETER (TWOP12=4096., IGIGA=1000000000,JSDFLT=314159265) PARAMETER (ITWO24=2**24, ICONS=2147483563) SAVE NOTYET, I24, J24, CARRY, SEEDS, TWOM24, TWOM12, LUXLEV SAVE NSKIP, NDSKIP, IN24, NEXT, KOUNT, MKOUNT, INSEED INTEGER LUXLEV LOGICAL NOTYET DATA NOTYET, LUXLEV, IN24, KOUNT, MKOUNT /.TRUE., LXDFLT, 0,0,0/ DATA I24,J24,CARRY/24,10,0./ C default C Luxury Level 0 1 2 *3* 4 DATA NDSKIP/0, 24, 73, 199, 365 / Corresponds to p=24 48 97 223 389 C time factor 1 2 3 6 10 on slow workstation C 1 1.5 2 3 5 on fast mainframe C C NOTYET is .TRUE. if no initialization has been performed yet. C Default Initialization by Multiplicative Congruential IF (NOTYET) THEN NOTYET = .FALSE. JSEED = JSDFLT INSEED = JSEED CC WRITE(6,'(A,I12)') ' RANLUX DEFAULT INITIALIZATION: ',JSEED LUXLEV = LXDFLT NSKIP = NDSKIP(LUXLEV) IN24 = 0 KOUNT = 0 MKOUNT = 0 CC LP = NSKIP + 24 CC WRITE(6,'(A,I2,A,I4)') ' RANLUX DEFAULT LUXURY LEVEL = ', CC + LUXLEV,' p =',LP TWOM24 = 1. DO 25 I= 1, 24 TWOM24 = TWOM24 * 0.5 K = JSEED/53668 JSEED = 40014*(JSEED-K*53668) -K*12211 IF (JSEED .LT. 0) JSEED = JSEED+ICONS ISEEDS(I) = MOD(JSEED,ITWO24) 25 CONTINUE TWOM12 = TWOM24 * 4096. DO 50 I= 1,24 SEEDS(I) = REAL(ISEEDS(I))*TWOM24 NEXT(I) = I-1 50 CONTINUE NEXT(1) = 24 I24 = 24 J24 = 10 CARRY = 0. IF (SEEDS(24) .EQ. 0.) CARRY = TWOM24 ENDIF C C The Generator proper: "Subtract-with-borrow", C as proposed by Marsaglia and Zaman, C Florida State University, March, 1989 C DO 100 IVEC= 1, LENV UNI = SEEDS(J24) - SEEDS(I24) - CARRY IF (UNI .LT. 0.) THEN UNI = UNI + 1.0 CARRY = TWOM24 ELSE CARRY = 0. ENDIF SEEDS(I24) = UNI I24 = NEXT(I24) J24 = NEXT(J24) RVEC(IVEC) = UNI C small numbers (with less than 12 "significant" bits) are "padded". IF (UNI .LT. TWOM12) THEN RVEC(IVEC) = RVEC(IVEC) + TWOM24*SEEDS(J24) C and zero is forbidden in case someone takes a logarithm IF (RVEC(IVEC) .EQ. 0.) RVEC(IVEC) = TWOM24*TWOM24 ENDIF C Skipping to luxury. As proposed by Martin Luscher. IN24 = IN24 + 1 IF (IN24 .EQ. 24) THEN IN24 = 0 KOUNT = KOUNT + NSKIP DO 90 ISK= 1, NSKIP UNI = SEEDS(J24) - SEEDS(I24) - CARRY IF (UNI .LT. 0.) THEN UNI = UNI + 1.0 CARRY = TWOM24 ELSE CARRY = 0. ENDIF SEEDS(I24) = UNI I24 = NEXT(I24) J24 = NEXT(J24) 90 CONTINUE ENDIF 100 CONTINUE KOUNT = KOUNT + LENV IF (KOUNT .GE. IGIGA) THEN MKOUNT = MKOUNT + 1 KOUNT = KOUNT - IGIGA ENDIF RETURN C C Entry to input and float integer seeds from previous run ENTRY RLUXIN(ISDEXT) * IF block added by Phillip Helbig after correpondence with James IF (NOTYET) THEN CC WRITE(6,'(A)') ' PROPER RESULTS ONLY WITH INITIALISATION FROM CC $25 INTEGERS OBTAINED WITH RLUXUT' NOTYET = .FALSE. ENDIF TWOM24 = 1. DO 195 I= 1, 24 NEXT(I) = I-1 195 TWOM24 = TWOM24 * 0.5 NEXT(1) = 24 TWOM12 = TWOM24 * 4096. CC WRITE(6,'(A)') ' FULL INITIALIZATION OF RANLUX WITH 25 INTEGERS:' CC WRITE(6,'(5X,5I12)') ISDEXT DO 200 I= 1, 24 SEEDS(I) = REAL(ISDEXT(I))*TWOM24 200 CONTINUE CARRY = 0. IF (ISDEXT(25) .LT. 0) CARRY = TWOM24 ISD = IABS(ISDEXT(25)) I24 = MOD(ISD,100) ISD = ISD/100 J24 = MOD(ISD,100) ISD = ISD/100 IN24 = MOD(ISD,100) ISD = ISD/100 LUXLEV = ISD IF (LUXLEV .LE. MAXLEV) THEN NSKIP = NDSKIP(LUXLEV) CC WRITE (6,'(A,I2)') ' RANLUX LUXURY LEVEL SET BY RLUXIN TO: ', CC + LUXLEV ELSE IF (LUXLEV .GE. 24) THEN NSKIP = LUXLEV - 24 CC WRITE (6,'(A,I5)') ' RANLUX P-VALUE SET BY RLUXIN TO:',LUXLEV ELSE NSKIP = NDSKIP(MAXLEV) CC WRITE (6,'(A,I5)') ' RANLUX ILLEGAL LUXURY RLUXIN: ',LUXLEV LUXLEV = MAXLEV ENDIF INSEED = -1 RETURN C C Entry to output seeds as integers ENTRY RLUXUT(ISDEXT) DO 300 I= 1, 24 ISDEXT(I) = INT(SEEDS(I)*TWOP12*TWOP12) 300 CONTINUE ISDEXT(25) = I24 + 100*J24 + 10000*IN24 + 1000000*LUXLEV IF (CARRY .GT. 0.) ISDEXT(25) = -ISDEXT(25) RETURN C C Entry to output the "convenient" restart point ENTRY RLUXAT(LOUT,INOUT,K1,K2) LOUT = LUXLEV INOUT = INSEED K1 = KOUNT K2 = MKOUNT RETURN C C Entry to initialize from one or three integers ENTRY RLUXGO(LUX,INS,K1,K2) IF (LUX .LT. 0) THEN LUXLEV = LXDFLT ELSE IF (LUX .LE. MAXLEV) THEN LUXLEV = LUX ELSE IF (LUX .LT. 24 .OR. LUX .GT. 2000) THEN LUXLEV = MAXLEV CC WRITE (6,'(A,I7)') ' RANLUX ILLEGAL LUXURY RLUXGO: ',LUX ELSE LUXLEV = LUX DO 310 ILX= 0, MAXLEV IF (LUX .EQ. NDSKIP(ILX)+24) LUXLEV = ILX 310 CONTINUE ENDIF IF (LUXLEV .LE. MAXLEV) THEN NSKIP = NDSKIP(LUXLEV) CC WRITE(6,'(A,I2,A,I4)') ' RANLUX LUXURY LEVEL SET BY RLUXGO :', CC + LUXLEV,' P=', NSKIP+24 ELSE NSKIP = LUXLEV - 24 CC WRITE (6,'(A,I5)') ' RANLUX P-VALUE SET BY RLUXGO TO:',LUXLEV ENDIF IN24 = 0 CC IF (INS .LT. 0) WRITE (6,'(A)') CC + ' Illegal initialization by RLUXGO, negative input seed' IF (INS .GT. 0) THEN JSEED = INS CC WRITE(6,'(A,3I12)') ' RANLUX INITIALIZED BY RLUXGO FROM SEEDS', CC + JSEED, K1,K2 ELSE JSEED = JSDFLT CC WRITE(6,'(A)')' RANLUX INITIALIZED BY RLUXGO FROM DEFAULT SEED' ENDIF INSEED = JSEED NOTYET = .FALSE. TWOM24 = 1. DO 325 I= 1, 24 TWOM24 = TWOM24 * 0.5 K = JSEED/53668 JSEED = 40014*(JSEED-K*53668) -K*12211 IF (JSEED .LT. 0) JSEED = JSEED+ICONS ISEEDS(I) = MOD(JSEED,ITWO24) 325 CONTINUE TWOM12 = TWOM24 * 4096. DO 350 I= 1,24 SEEDS(I) = REAL(ISEEDS(I))*TWOM24 NEXT(I) = I-1 350 CONTINUE NEXT(1) = 24 I24 = 24 J24 = 10 CARRY = 0. IF (SEEDS(24) .EQ. 0.) CARRY = TWOM24 C If restarting at a break point, skip K1 + IGIGA*K2 C Note that this is the number of numbers delivered to C the user PLUS the number skipped (if luxury .GT. 0). KOUNT = K1 MKOUNT = K2 IF (K1+K2 .NE. 0) THEN DO 500 IOUTER= 1, K2+1 INNER = IGIGA IF (IOUTER .EQ. K2+1) INNER = K1 DO 450 ISK= 1, INNER UNI = SEEDS(J24) - SEEDS(I24) - CARRY IF (UNI .LT. 0.) THEN UNI = UNI + 1.0 CARRY = TWOM24 ELSE CARRY = 0. ENDIF SEEDS(I24) = UNI I24 = NEXT(I24) J24 = NEXT(J24) 450 CONTINUE 500 CONTINUE C Get the right value of IN24 by direct calculation IN24 = MOD(KOUNT, NSKIP+24) IF (MKOUNT .GT. 0) THEN IZIP = MOD(IGIGA, NSKIP+24) IZIP2 = MKOUNT*IZIP + IN24 IN24 = MOD(IZIP2, NSKIP+24) ENDIF C Now IN24 had better be between zero and 23 inclusive IF (IN24 .GT. 23) THEN CC WRITE (6,'(A/A,3I11,A,I5)') CC + ' Error in RESTARTING with RLUXGO:',' The values', INS, CC + K1, K2, ' cannot occur at luxury level', LUXLEV IN24 = 0 ENDIF ENDIF RETURN END C --------------------------------------------------------------------- FUNCTION RANM ( ) C --------------------------------------------------------------------- C C Description: C Function RANM returns a random number, uniformly distributed C in the interval (0,1). It uses the RANLUX random number generator C due to Luscher and James. For details see: C Luscher M., Comput. Phys. Commun. 79 (1994) 100. C James F., Comput. Phys. Commun. 79 (1994) 111. C C --------------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) C Local variables: PARAMETER ( LVEC = 1 ) DIMENSION RVEC(LVEC) C --------------------------------------------------------------------- C CALL RANLUX(RVEC,LVEC) RANM = RVEC(1) C END C --------------------------------------------------------------------- FUNCTION RANDI ( A, B ) C --------------------------------------------------------------------- C C Description: C Function RANDI returns a random number between A and B. C C --------------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) C R = RANM() RANDI = A+(B-A)*R C END C ---------------------------------------------------------------------