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),IXAT(N),ICODE(4)                         
      DIMENSION BESTP(N,M),BESTV(M)                                     
      CHARACTER CFO*50
      DATA XL/(N)*-10.D0/                                                
      DATA XR/(N)*10.D0/                                                 
      DATA IXAT/(N)*1/                                                  
      DATA ICODE/1,1,1,0/                                               
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            
C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C  PREPARE THE INTERFACE TO MERLIN (LOCAL MINIMIZER).
C            
            NOF = 0                                                     
            ND  = N                                                     
C            
C  INSTRUCT MERLIN TO OPTIMIZE FROM THE GLOBALY BEST POINT FOUND 
C            
            CALL OPTIMA(ND,NOF,XB,VAL,XL,XR,IXAT,ICODE,                 
     >          'in.dat','/dev/null',GRMSO,NF,NOG,NHO,NJO)              
C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C            
            WRITE(*,88)'MERLIN:'                                        
            WRITE(*,CFO) NF,VAL,(XB(J),J=1,N)                            
            WRITE(*,'(57(''=''))')                                      
            WRITE(*,*)'TOTAL FUNCTION CALLS: ',M*NITER+NF               
            WRITE(*,*)'TOTAL GRADIENT CALLS: ',NOG                      
      END                                                               
