C     !----------------- MAIN PROGRAM : ORDCANON ----------------------
C     PROVIDES TO USE REPULSIVE PARTICLE SWARM METHOD TO
C     OBTAIN THE LARGEST CANONICAL CORRELATION & COMPOSITE VARIATE RANKS
C     PRODUCT MOMENT AS WELL AS ABSOLUTE CORRELATION (BRADLEY, 1985) MAY
C     BE USED. PROGRAM BY SK MISHRA, DEPT. OF ECONOMICS, NORTH-EASTERN
C     HILL UNIVERSITY, SHILLONG (INDIA)
C     -----------------------------------------------------------------
C         ADJUST THE PARAMETERS SUITABLY IN SUBROUTINES RPS
C     WHEN THE PROGRAM ASKS FOR PARAMETERS, FEED THEM SUITABLY
C     -----------------------------------------------------------------
      PROGRAM ORDCANON
      PARAMETER(NOB=30,MVAR=9)!CHANGE THE PARAMETERS HERE AS NEEDED.
C     -----------------------------------------------------------------
C     NOB=NO. OF CASES AND MVAR=NO. OF VARIABLES IN ALL M= (M1+M2)
C     TO BE ADJUSTED IN SUBROUTINE CORD(M,X,F) ALSO: STATEMENT 931
      IMPLICIT DOUBLE PRECISION (A-H, O-Z)
      COMMON /KFF/KF,NFCALL,FTIT ! FUNCTION CODE, NO. OF CALLS & TITLE
      CHARACTER *30  METHOD(1)
      CHARACTER *70 FTIT
      CHARACTER *40 INFILE,OUTFILE
      COMMON /CANON/MONE,MTWO
      COMMON /CORDAT/CDAT(NOB,MVAR),QIND1(NOB),QIND2(NOB),R(1),NORM,NCOR
      COMMON /XBASE/XBAS
      COMMON /RNDM/IU,IV ! RANDOM NUMBER GENERATION (IU = 4-DIGIT SEED)
      COMMON /GETRANK/MRNK
      INTEGER IU,IV
      DIMENSION XX(3,50),KKF(3),MM(3),FMINN(3),XBAS(1000,50)
      DIMENSION ZDAT(NOB,MVAR+1),FRANK1(NOB),FRANK2(NOB),RMAT(2,2)
      DIMENSION X(50)! X IS THE DECISION VARIABLE X IN F(X) TO MINIMIZE
C       M = DIMENSION OF THE PROBLEM, KF(=1) = TEST FUNCTION CODE AND
C              FMIN IS THE MIN VALUE OF F(X) OBTAINED FROM RPS
      WRITE(*,*)'====================     WARNING    =============== '
      WRITE(*,*)'ADJUST PARAMETERS IN SUBROUTINES RPS IF NEEDED '
C     ------------------ OPTIMIZATION BY RPS METHOD -------------------
      NORM=2!WORKS WITH THE EUCLIDEAN NORM (IDENTICAL RESULTS IF NORM=1)
      NOPT=1 ! ONLY ONE FUNCTION IS OPTIMIZED
      WRITE(*,*)'=================================================== '
      METHOD(1)=' : REPULSIVE PARTICLE SWARM OPTIMIZATION'
C     INITIALIZE. THIS XBAS WILL BE USED TO INITIALIZE THE POPULATION.
      WRITE(*,*)' '
      WRITE(*,*)'---------- FEED RANDOM NUMBER SEED, AND NCOR ---------'
      WRITE(*,*)' '
      WRITE(*,*)'FEED SEED [ANY 4-DIGIT NUMBER] AND NCOR[0,1]'
      WRITE(*,*)'NCOR(0)=PRODUCT MOMENT; NCOR(1)=ABSOLUTE CORRELATION'
      WRITE(*,*)' '
    1 READ(*,*) IU,NCOR
      IF(NCOR.LT.0.OR.NCOR.GT.1) THEN
      WRITE(*,*)'SORRY. NCOR TAKES ON[0,1] ONLY. FEED SEED & NCOR AGAIN'
      GOTO 1
      ENDIF
      WRITE(*,*)'WANT RANK SCORE OPTIMIZATION? YES(1); NO(OTHER THAN 1)'
      READ(*,*) MRNK
      WRITE(*,*)'INPUT FILE TO READ DATA:YOUR DATA MUST BE IN THIS FILE'
      WRITE(*,*)'CASES (NOB) IN ROWS ; VARIABLES (MVAR) IN COLUMNS'
      READ(*,*) INFILE
      WRITE(*,*)'SPECIFY THE OUTPUT FILE TO STORE THE RESULTS'
      READ(*,*) OUTFILE
      OPEN(9, FILE=OUTFILE)
      OPEN(7,FILE=INFILE)
      DO I=1,NOB
      READ(7,*),CDA,(CDAT(I,J),J=1,MVAR)
      ENDDO
      CLOSE(7)
      DO I=1,NOB
      DO J=1,MVAR
      ZDAT(I,J+1)=CDAT(I,J)
      ENDDO
      ENDDO
      WRITE(*,*)'DATA HAS BEEN READ. WOULD YOU UNITIZE VARIABLES? [YES=1
     & ELSE NO UNITIZATION]'
      WRITE(*,*)'UNITIZE MEANS TRANSFORMATION FROM X(I,J) TO UNITIZED X'
      WRITE(*,*)'[X(I,J)-MIN(X(.,J))]/[MAX(X(.,J))-MIN(X(.,J))]'
      READ(*,*) NUN
      IF(NUN.EQ.1) THEN
      DO J=1,MVAR
      CMIN=CDAT(1,J)
      CMAX=CDAT(1,J)
      DO I=2,NOB
      IF(CMIN.GT.CDAT(I,J)) CMIN=CDAT(I,J)
      IF(CMAX.LT.CDAT(I,J)) CMAX=CDAT(I,J)
      ENDDO
      DO I=1,NOB
      CDAT(I,J)=(CDAT(I,J)-CMIN)/(CMAX-CMIN)
      ENDDO
      ENDDO
      ENDIF
C     ------------------------------------------------------------------
C     THIS XBAS WILL BE USED AS INITIAL X
      DO I=1,1000
      DO J=1,50
      CALL RANDOM(RAND)
      XBAS(I,J)=RAND ! RANDOM NUMBER BETWEEN (0, 1)
      ENDDO
      ENDDO
C     ------------------------------------------------------------------
      WRITE(*,*)' *****************************************************'
C     ------------------------------------------------------------------
      K=1
      WRITE(*,*)'PARTICLE SWARM PROGRAM TO OBTAIN CANONICAL CORRELATION'
      CALL RPS(M,X,FMINRPS,Q1) !CALLS RPS AND RETURNS OPTIMAL X AND FMIN
      WRITE(*,*)'RPS BRINGS THE FOLLOWING VALUES TO THE MAIN PROGRAM'
      WRITE(*,*)(X(JOPT),JOPT=1,M),' OPTIMUM FUNCTION=',FMINRPS
      IF(KF.EQ.1) THEN
      WRITE(9,*)'REPULSIVE PARTICLE SWARM OPTIMIZATION RESULTS'
      WRITE(9,*)'THE LARGEST CANONICAL R BETWEEN THE SETS OF VARIABLES'
      WRITE(9,*)' ABS(R)=',DABS(R(1)),'; SQUARE(R)=',R(1)**2
      IF(NCOR.EQ.0) THEN
      WRITE(*,*)'NOTE: THESE ARE KARL PEARSON TYPE CORRELATION (NCOR=0)'
      WRITE(9,*)'NOTE: THESE ARE KARL PEARSON TYPE CORRELATION (NCOR=0)'
      ELSE
      WRITE(*,*)'NOTE: THESE ARE BRADLEY TYPE CORRELATION (NCOR=1)'
      WRITE(9,*)'NOTE: THESE ARE BRADLEY TYPE CORRELATION (NCOR=1)'
      ENDIF
      WRITE(*,*)'______________________________________________________'
      WRITE(9,*)'______________________________________________________'

      DO II=1,NOB
      FRANK1(II)=QIND1(II)
      FRANK2(II)=QIND2(II)
      ENDDO
      ENDIF
      FMIN=FMINRPS
C     ------------------------------------------------------------------
      DO J=1,M
      XX(K,J)=X(J)
      ENDDO
      KKF(K)=KF
      MM(K)=M
      FMINN(K)=FMIN
      WRITE(*,*)' '
      WRITE(*,*)' '
      WRITE(*,*)'---------------------- FINAL RESULTS=================='
      WRITE(*,*)'FUNCT CODE=',KKF(K),'  FMIN=',FMINN(K),' : DIM=',MM(K)
      WRITE(*,*)'OPTIMAL DECISION VARIABLES : ',METHOD(K)
      WRITE(*,*)'FOR THE FIRST SET OF VARIABLES WEIGHTS ARE AS FOLLOWS'
      WRITE(9,*)'FOR THE FIRST SET OF VARIABLES WEIGHTS ARE AS FOLLOWS'
      WRITE(9,*)(XX(K,J),J=1,MONE)
      WRITE(*,*)(XX(K,J),J=1,MONE)
      WRITE(*,*)'FOR THE SECOND SET OF VARIABLES WEIGHTS ARE AS FOLLOWS'
      WRITE(9,*)'FOR THE SECOND SET OF VARIABLES WEIGHTS ARE AS FOLLOWS'
      WRITE(9,*)(XX(K,J),J=MONE+1,M)
      WRITE(*,*)(XX(K,J),J=MONE+1,M)
      WRITE(*,*)'/////////////////////////////////////////////////////'
      WRITE(*,*)'OPTIMIZATION PROGRAM ENDED'
      WRITE(*,*)'******************************************************'
      WRITE(*,*)'MEASURE OF EQUALITY/INEQUALITY'
      WRITE(*,*)'RPS: BEFORE AND AFTER OPTIMIZATION = ',Q0,Q1
      WRITE(*,*)' '
      WRITE(*,*)'RESULTS STORED IN FILE= ',OUTFILE
      WRITE(*,*)'OPEN BY MSWORD OR EDIT OR ANY OTHER EDITOR'
      WRITE(*,*)' '
      WRITE(*,*)'NOTE:VECTORS OF CORRELATIONS & INDEX(BOTH TOGETHER) ARE
     & IDETERMINATE FOR SIGN &   MAY BE MULTIPLED BY (-1) IF NEEDED'
      WRITE(*,*)'THAT IS IF R(J) IS TRANSFORMED TO -R(J) FOR ALL J THEN
     &THE INDEX(I) TOO IS         TRANSFORMED TO -INDEX(I) FOR ALL I'
      WRITE(9,*)' '
      WRITE(9,*)'NOTE: VECTORS OF CORRELATIONS AND INDEX (BOTH TOGETHER)
     & ARE IDETERMINATE FOR SIGN AND MAY BE MULTIPLED BY (-1) IF NEEDED'
      WRITE(9,*)'THAT IS IF R(J) IS TRANSFORMED TO -R(J) FOR ALL J THEN
     &THE INDEX(I) TOO IS TRANSFORMED TO -INDEX(I) FOR ALL I'
      CALL DORANK(FRANK1,NOB)
      CALL DORANK(FRANK2,NOB)
      DO I=1,NOB
      ZDAT(I,1)=FRANK1(I)
      ZDAT(I,2)=FRANK2(I)
      ENDDO
      IF(NCOR.EQ.0) THEN
      CALL CORREL(ZDAT,NOB,2,RMAT)
      ELSE
      CALL DOCORA(ZDAT,NOB,2,RMAT)
      ENDIF
      WRITE(9,*)'=================================================== '
      WRITE(*,*)'=================================================== '
      WRITE(9,*)'1ST 2 ARE CANONICAL SCORES AND LAST 2 ARE THEIR RANK'
      WRITE(*,*)'1ST 2 ARE CANONICAL SCORES AND LAST 2 ARE THEIR RANK'
      WRITE(9,*)'=================================================== '
      WRITE(*,*)'=================================================== '
      DO I=1,NOB
      IF(MRNK.EQ.1) THEN
      QIND1(I)=0.D0
      QIND2(I)=0.D0
      DO J=1,MONE
      QIND1(I)=QIND1(I)+ZDAT(I,J+1)*XX(NOPT,J)
      ENDDO
      DO J=MONE+1,MVAR
      QIND2(I)=QIND2(I)+ZDAT(I,J+1)*XX(NOPT,J)
      ENDDO
      ENDIF
      WRITE(9,2)I,QIND1(I),QIND2(I),(ZDAT(I,J),J=1,2)
      WRITE(*,2)I,QIND1(I),QIND2(I),(ZDAT(I,J),J=1,2)
      ENDDO
    2 FORMAT(I5,2F15.6,2F10.3)
      WRITE(9,*)'SQUARE OF CANONICAL CORRELATION =',RMAT(1,2)**2
      WRITE(*,*)'SQUARE OF CANONICAL CORRELATION =',RMAT(1,2)**2
      WRITE(9,*)'ABSOLUTE OF CANONICAL CORRELATION =',DABS(RMAT(1,2))
      WRITE(*,*)'ABSOLUTE OF CANONICAL CORRELATION =',DABS(RMAT(1,2))
      IF(NCOR.EQ.0) THEN
      WRITE(*,*)'NOTE: THESE ARE KARL PEARSON TYPE CORRELATION (NCOR=0)'
      WRITE(9,*)'NOTE: THESE ARE KARL PEARSON TYPE CORRELATION (NCOR=0)'
      ELSE
      WRITE(*,*)'NOTE: THESE ARE BRADLEY TYPE CORRELATION (NCOR=1)'
      WRITE(9,*)'NOTE: THESE ARE BRADLEY TYPE CORRELATION (NCOR=1)'
      ENDIF
      CLOSE(9)
      WRITE(*,*) 'THE JOB IS OVER'
      END
C     -----------------------------------------------------------------
      SUBROUTINE RPS(M,ABEST,FBEST,G1)
C     PROGRAM TO FIND GLOBAL MINIMUM BY REPULSIVE PARTICLE SWARM METHOD
C     WRITTEN BY SK MISHRA, DEPT. OF ECONOMICS, NEHU, SHILLONG (INDIA)
C     -----------------------------------------------------------------
      PARAMETER (N=50,NN=10,MX=100,NSTEP=7,ITRN=10000,NSIGMA=1,ITOP=1)
      PARAMETER (NPRN=50) ! DISPLAYS RESULTS AT EVERY 500 TH ITERATION
C     PARAMETER(N=50,NN=25,MX=100,NSTEP=9,ITRN=10000,NSIGMA=1,ITOP=3)
C     PARAMETER (N=100,NN=15,MX=100,NSTEP=9,ITRN=10000,NSIGMA=1,ITOP=3)
C     IN CERTAIN CASES THE ONE OR THE OTHER SPECIFICATION WORKS BETTER
C     DIFFERENT SPECIFICATIONS OF PARAMETERS MAY SUIT DIFFERENT TYPES
C     OF FUNCTIONS OR DIMENSIONS - ONE HAS TO DO SOME TRIAL AND ERROR
C     -----------------------------------------------------------------
C     N = POPULATION SIZE. IN MOST OF THE CASES N=30 IS OK. ITS VALUE
C     MAY BE INCREASED TO 50 OR 100 TOO. THE PARAMETER NN IS THE SIZE OF
C     RANDOMLY CHOSEN NEIGHBOURS. 15 TO 25 (BUT SUFFICIENTLY LESS THAN
C     N) IS A GOOD CHOICE. MX IS THE MAXIMAL SIZE OF DECISION VARIABLES.
C     IN F(X1, X2,...,XM) M SHOULD BE LESS THAN OR EQUAL TO MX. ITRN IS
C     THE NO. OF ITERATIONS. IT MAY DEPEND ON THE PROBLEM. 200(AT LEAST)
C     TO 500 ITERATIONS MAY BE GOOD ENOUGH. BUT FOR FUNCTIONS LIKE
C     ROSENBROCKOR GRIEWANK OF LARGE SIZE (SAY M=30) IT IS NEEDED THAT
C     ITRN IS LARGE, SAY 5000 OR EVEN 10000.
C     SIGMA INTRODUCES PERTURBATION & HELPS THE SEARCH JUMP OUT OF LOCAL
C     OPTIMA. FOR EXAMPLE : RASTRIGIN FUNCTION OF DMENSION 3O OR LARGER
C     NSTEP DOES LOCAL SEARCH BY TUNNELLING AND WORKS WELL BETWEEN 5 AND
C     15, WHICH IS MUCH ON THE HIGHER SIDE.
C     ITOP <=1 (RING); ITOP=2 (RING AND RANDOM); ITOP=>3 (RANDOM)
C     NSIGMA=0 (NO CHAOTIC PERTURBATION);NSIGMA=1 (CHAOTIC PERTURBATION)
C     NOTE THAT NSIGMA=1 NEED NOT ALWAYS WORK BETTER (OR WORSE)
C     SUBROUTINE FUNC( ) DEFINES OR CALLS THE FUNCTION TO BE OPTIMIZED.
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /RNDM/IU,IV
      COMMON /KFF/KF,NFCALL,FTIT
      INTEGER IU,IV
      CHARACTER *70 FTIT
      DIMENSION X(N,MX),V(N,MX),A(MX),VI(MX),TIT(50),ABEST(*)
      DIMENSION XX(N,MX),F(N),V1(MX),V2(MX),V3(MX),V4(MX),BST(MX)
C     A1 A2 AND A3 ARE CONSTANTS AND W IS THE INERTIA WEIGHT.
C     OCCASIONALLY, TINKERING WITH THESE VALUES, ESPECIALLY A3, MAY BE
C     NEEDED.
      DATA A1,A2,A3,W,SIGMA /.5D00,.5D00,.0005D00,.5D00,1.D-03/
      EPSILON=1.D-12 ! ACCURACY NEEDED FOR TERMINATON
C     --------------------CHOOSING THE TEST FUNCTION ------------------'
      CALL FSELECT(KF,M,FTIT)
C     -----------------------------------------------------------------
      FFMIN=1.D30
      LCOUNT=0
      NFCALL=0
      WRITE(*,*)'4-DIGITS SEED FOR RANDOM NUMBER GENERATION'
      READ(*,*) IU
      DATA FMIN /1.0E30/
C     GENERATE N-SIZE POPULATION OF M-TUPLE PARAMETERS X(I,J) RANDOMLY
      DO I=1,N
        DO J=1,M
        CALL RANDOM(RAND)
         X(I,J)=RAND
C     WE GENERATE RANDOM(-5,5). HERE MULTIPLIER IS 10. TINKERING IN SOME
C     CASES MAY BE NEEDED
         ENDDO
        F(I)=1.0D30
      ENDDO
C     INITIALISE VELOCITIES V(I) FOR EACH INDIVIDUAL IN THE POPULATION
      DO I=1,N
      DO J=1,M
      CALL RANDOM(RAND)
       V(I,J)=(RAND-0.5D+00)
C       V(I,J)=RAND
      ENDDO
      ENDDO
      DO 100 ITER=1,ITRN
C     WRITE(*,*)'ITERATION=',ITER
C     LET EACH INDIVIDUAL SEARCH FOR THE BEST IN ITS NEIGHBOURHOOD
        DO I=1,N
           DO J=1,M
           A(J)=X(I,J)
           VI(J)=V(I,J)
           ENDDO
           CALL LSRCH(A,M,VI,NSTEP,FI)
           IF(FI.LT.F(I)) THEN
            F(I)=FI
            DO IN=1,M
            BST(IN)=A(IN)
            ENDDO
C     F(I) CONTAINS THE LOCAL BEST VALUE OF FUNCTION FOR ITH INDIVIDUAL
C     XX(I,J) IS THE M-TUPLE VALUE OF X ASSOCIATED WITH LOCAL BEST F(I)
             DO J=1,M
             XX(I,J)=A(J)
             ENDDO
             ENDIF
        ENDDO
C      NOW LET EVERY INDIVIDUAL RANDOMLY COSULT NN(<<N) COLLEAGUES AND
C      FIND THE BEST AMONG THEM
      DO I=1,N
C     ------------------------------------------------------------------
      IF(ITOP.GE.3) THEN
C     RANDOM TOPOLOGY ******************************************
C     CHOOSE NN COLLEAGUES RANDOMLY AND FIND THE BEST AMONG THEM
          BEST=1.0D30
           DO II=1,NN
                 CALL RANDOM(RAND)
               NF=INT(RAND*N)+1
                IF(BEST.GT.F(NF)) THEN
                 BEST=F(NF)
                NFBEST=NF
                 ENDIF
            ENDDO
      ENDIF
C----------------------------------------------------------------------
      IF(ITOP.EQ.2) THEN
C     RING + RANDOM TOPOLOGY ******************************************
C     REQUIRES THAT THE SUBROUTINE NEIGHBOR IS TURNED ALIVE
       BEST=1.0D30
          CALL NEIGHBOR(I,N,I1,I3)
          DO II=1,NN
                IF(II.EQ.1) NF=I1
                 IF(II.EQ.2) NF=I
                  IF(II.EQ.3) NF=I3
                      IF(II.GT.3) THEN
                     CALL RANDOM(RAND)
                      NF=INT(RAND*N)+1
                     ENDIF
                  IF(BEST.GT.F(NF)) THEN
                  BEST=F(NF)
                  NFBEST=NF
                 ENDIF
             ENDDO
       ENDIF
C---------------------------------------------------------------------
      IF(ITOP.LE.1) THEN
C     RING TOPOLOGY **************************************************
C     REQUIRES THAT THE SUBROUTINE NEIGHBOR IS TURNED ALIVE
        BEST=1.0D30
           CALL NEIGHBOR(I,N,I1,I3)
              DO II=1,3
              IF (II.NE.I) THEN
             IF(II.EQ.1) NF=I1
              IF(II.EQ.3) NF=I3
                  IF(BEST.GT.F(NF)) THEN
                   BEST=F(NF)
                   NFBEST=NF
                  ENDIF
                  ENDIF
            ENDDO
       ENDIF
C---------------------------------------------------------------------
C     IN THE LIGHT OF HIS OWN AND HIS BEST COLLEAGUES EXPERIENCE, THE
C     INDIVIDUAL I WILL MODIFY HIS MOVE AS PER THE FOLLOWING CRITERION
C     FIRST, ADJUSTMENT BASED ON ONES OWN EXPERIENCE
C     AND OWN BEST EXPERIENCE IN THE PAST (XX(I))
           DO J=1,M
           CALL RANDOM(RAND)
           V1(J)=A1*RAND*(XX(I,J)-X(I,J))

C     THEN BASED ON THE OTHER COLLEAGUES BEST EXPERIENCE WITH WEIGHT W
C     HERE W IS CALLED AN INERTIA WEIGHT 0.01< W < 0.7
C     A2 IS THE CONSTANT NEAR BUT LESS THAN UNITY
           CALL RANDOM(RAND)
           V2(J)=V(I,J)
           IF(F(NFBEST).LT.F(I)) THEN
           V2(J)=A2*W*RAND*(XX(NFBEST,J)-X(I,J))
           ENDIF
C     THEN SOME RANDOMNESS AND A CONSTANT A3 CLOSE TO BUT LESS THAN UNITY
           CALL RANDOM(RAND)
           RND1=RAND
           CALL RANDOM(RAND)
            V3(J)=A3*RAND*W*RND1
C            V3(J)=A3*RAND*W
C     THEN ON PAST VELOCITY WITH INERTIA WEIGHT W
           V4(J)=W*V(I,J)
C     FINALLY A SUM OF THEM
           V(I,J)= V1(J)+V2(J)+V3(J)+V4(J)
           ENDDO
      ENDDO
C     CHANGE X
      DO I=1,N
      DO J=1,M
      RANDS=0.D00
C     ------------------------------------------------------------------
      IF(NSIGMA.EQ.1) THEN
       CALL RANDOM(RAND) ! FOR CHAOTIC PERTURBATION
       IF(DABS(RAND-.5D00).LT.SIGMA) RANDS=RAND-0.5D00
C     SIGMA CONDITIONED RANDS INTRODUCES CHAOTIC ELEMENT IN TO LOCATION
C     IN SOME CASES THIS PERTURBATION HAS WORKED VERY EFFECTIVELY WITH
C     PARAMETER (N=100,NN=15,MX=100,NSTEP=9,ITRN=100000,NSIGMA=1,ITOP=2)
      ENDIF
C     -----------------------------------------------------------------
      X(I,J)=X(I,J)+V(I,J)*(1.D00+RANDS)
      ENDDO
      ENDDO
       DO I=1,N
         IF(F(I).LT.FMIN) THEN
         FMIN=F(I)
         II=I
         DO J=1,M
         BST(J)=XX(II,J)
         ENDDO
         ENDIF
         ENDDO

      IF(LCOUNT.EQ.NPRN) THEN
      LCOUNT=0
      WRITE(*,*)'OPTIMAL SOLUTION UPTO THIS (FUNCTION CALLS=',NFCALL,')'
      WRITE(*,*)'X = ',(BST(J),J=1,M),' MIN F = ',FMIN
C      WRITE(*,*)'NO. OF FUNCTION CALLS = ',NFCALL
      DO J=1,M
      ABEST(J)=BST(J)
      ENDDO
      IF(DABS(FFMIN-FMIN).LT.EPSILON) GOTO 999
      FFMIN=FMIN
      ENDIF
      LCOUNT=LCOUNT+1
  100 CONTINUE
  999 WRITE(*,*)'------------------------------------------------------'
      DO I=1,N
      IF(F(I).LT.FBEST) THEN
      FBEST=F(I)
      DO J=1,M
      ABEST(J)=XX(I,J)
      ENDDO
      ENDIF
      ENDDO
      CALL FUNC(ABEST,M,FBEST)
      CALL GINI(F,N,G1)
      WRITE(*,*)'FINAL X = ',(BST(J),J=1,M),' FINAL MIN F = ',FMIN
      WRITE(*,*)'COMPUTATION OVER:FOR ',FTIT
      WRITE(*,*)'NO. OF VARIABLES=',M,'  END.'
      RETURN
      END
C     ----------------------------------------------------------------
      SUBROUTINE LSRCH(A,M,VI,NSTEP,FI)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      CHARACTER *70 FTIT
      COMMON /KFF/KF,NFCALL,FTIT
      COMMON /RNDM/IU,IV
      INTEGER IU,IV
      DIMENSION A(*),B(100),VI(*)
      AMN=1.0D30
      DO J=1,NSTEP
         DO JJ=1,M
         B(JJ)=A(JJ)+(J-(NSTEP/2)-1)*VI(JJ)
         ENDDO
      CALL FUNC(B,M,FI)
        IF(FI.LT.AMN) THEN
        AMN=FI
        DO JJ=1,M
        A(JJ)=B(JJ)
        ENDDO
        ENDIF
      ENDDO
      FI=AMN
      RETURN
      END
C     -----------------------------------------------------------------
C     THIS SUBROUTINE IS NEEDED IF THE NEIGHBOURHOOD HAS RING TOPOLOGY
C     EITHER PURE OR HYBRIDIZED
       SUBROUTINE NEIGHBOR(I,N,J,K)
       IF(I-1.GE.1 .AND. I.LT.N) THEN
       J=I-1
       K=I+1
       ELSE
       IF(I-1.LT.1) THEN
       J=N-I+1
       K=I+1
       ENDIF
       IF(I.EQ.N) THEN
       J=I-1
       K=1
       ENDIF
       ENDIF
       RETURN
       END
C     -----------------------------------------------------------------
C     RANDOM NUMBER GENERATOR (UNIFORM BETWEEN 0 AND 1 - BOTH EXCLUSIVE)
      SUBROUTINE RANDOM(RAND1)
       DOUBLE PRECISION  RAND1
       COMMON /RNDM/IU,IV
      INTEGER IU,IV
       IV=IU*65539
       IF(IV.LT.0) THEN
       IV=IV+2147483647+1
       ENDIF
       RAND=IV
       IU=IV
       RAND=RAND*0.4656613E-09
       RAND1= DBLE(RAND)
       RETURN
       END
C     -----------------------------------------------------------------
      SUBROUTINE GINI(F,N,G)
      PARAMETER (K=1) !K=1 GINI COEFFICENT; K=2 COEFFICIENT OF VARIATION
C     THIS PROGRAM COMPUTES MEASURE OF INEQUALITY
C     IF K =1 GET THE GINI COEFFICIENT. IF K=2 GET COEFF OF VARIATIONE
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION F(*)
      S=0.D0
      DO I=1,N
      S=S+F(I)
      ENDDO
      S=S/N
      H=0.D00
      DO I=1,N-1
      DO J=I+1,N
      H=H+(DABS(F(I)-F(J)))**K
      ENDDO
      ENDDO
      H=(H/(N**2))**(1.D0/K)! FOR K=1 H IS MEAN DEVIATION;
C                             FOR K=2 H IS STANDARD DEVIATION
      WRITE(*,*)'MEASURES OF DISPERSION AND CENTRAL TENDENCY = ',G,S
      G=DEXP(-H)! G IS THE MEASURE OF EQUALITY (NOT GINI OR CV)
C     G=H/DABS(S) !IF S NOT ZERO, K=1 THEN G=GINI, K=2 G=COEFF VARIATION
      RETURN
      END
C     -----------------------------------------------------------------
      SUBROUTINE FSELECT(KF,M,FTIT)
      COMMON /CANON/MONE,MTWO
C      THE PROGRAM REQUIRES INPUTS FROM THE USER ON THE FOLLOWING ------
C     (1) FUNCTION CODE (KF), (2) NO. OF VARIABLES IN THE FUNCTION (M);
      CHARACTER *70 TIT(100),FTIT
      NFN=1
      KF=1
      WRITE(*,*)'----------------------------------------------------'
      DATA TIT(1)/'COMPUTE CANONICAL CORRELATION FROM 2 DATA SETS'/
C     -----------------------------------------------------------------
      DO I=1,NFN
      WRITE(*,*)TIT(I)
      ENDDO
      WRITE(*,*)'----------------------------------------------------'
      WRITE(*,*)'SPECIFY NO. OF VARIABLES IN SET-1[=M1] AND SET-2[=M2]'
      READ(*,*) MONE, MTWO
      M=MONE+MTWO
      FTIT=TIT(KF) ! STORE THE NAME OF THE CHOSEN FUNCTION IN FTIT
      RETURN
      END
C     -----------------------------------------------------------------
      SUBROUTINE FUNC(X,M,F)
C     TEST FUNCTIONS FOR GLOBAL OPTIMIZATION PROGRAM
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /RNDM/IU,IV
      COMMON /KFF/KF,NFCALL,FTIT
      INTEGER IU,IV
      DIMENSION X(*)
      CHARACTER *70 FTIT
      NFCALL=NFCALL+1 ! INCREMENT TO NUMBER OF FUNCTION CALLS
C     KF IS THE CODE OF THE TEST FUNCTION
      IF(KF.EQ.1) THEN
      CALL CORD(M,X,F)
      RETURN
      ENDIF
C     =================================================================
      WRITE(*,*)'FUNCTION NOT DEFINED. PROGRAM ABORTED'
      STOP
      END
C     ------------------------------------------------------------------
      SUBROUTINE CORD(M,X,F)
      PARAMETER (NOB=30,MVAR=9)! CHANGE THE PARAMETERS HERE AS NEEDED.
C     ------------------------------------------------------------------
C     NOB=NO. OF OBSERVATIONS (CASES) & MVAR= NO. OF VARIABLES
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /CANON/MONE,MTWO
      COMMON /RNDM/IU,IV
      COMMON /CORDAT/CDAT(NOB,MVAR),QIND1(NOB),QIND2(NOB),R(1),NORM,NCOR
      COMMON /GETRANK/MRNK
      INTEGER IU,IV
      DIMENSION X(*),Z(NOB,2)
      DO I=1,M
      IF(X(I).LT.-1.0D0.OR.X(I).GT.1.0D0) THEN
      CALL RANDOM(RAND)
      X(I)=(RAND-0.5D0)*2
      ENDIF
      ENDDO
      XNORM=0.D0
      DO J=1,M
      XNORM=XNORM+X(J)**2
      ENDDO
      XNORM=DSQRT(XNORM)
      DO J=1,M
      X(J)=X(J)/XNORM
      ENDDO
C     CONSTRUCT INDEX
      DO I=1,NOB
      QIND1(I)=0.D0
      QIND2(I)=0.D0
      DO J=1,MONE
      QIND1(I)=QIND1(I)+CDAT(I,J)*X(J)
      ENDDO
      DO J=MONE+1,M
      QIND2(I)=QIND2(I)+CDAT(I,J)*X(J)
      ENDDO
      ENDDO

C     ------------------------------------------------------------------
C      !FIND THE RANK OF QIND
      IF(MRNK.EQ.1) THEN
      CALL DORANK(QIND1,NOB)
      CALL DORANK(QIND2,NOB)
      ENDIF
C     ------------------------------------------------------------------
C     COMPUTE CORRELATIONS
      DO I=1,NOB
      Z(I,1)=QIND1(I)
      Z(I,2)=QIND2(I)
      ENDDO

      IF(NCOR.EQ.0) THEN
      CALL CORLN(Z,NOB,RHO)
      ELSE
      CALL CORA(Z,NOB,RHO)
      ENDIF
      R(1)=RHO
      F= DABS(R(1))**NORM
C     ------------------------------------------------------------------
      F=-F
      RETURN
      END
      SUBROUTINE CORLN(Z,NOB,RHO)
C     NOB = NO. OF CASES
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION Z(NOB,2),AV(2),SD(2)
      DO J=1,2
      AV(J)=0.D0
      SD(J)=0.D0
      DO I=1,NOB
      AV(J)=AV(J)+Z(I,J)
      SD(J)=SD(J)+Z(I,J)**2
      ENDDO
      ENDDO
      DO J=1,2
      AV(J)=AV(J)/NOB
      SD(J)=DSQRT(SD(J)/NOB-AV(J)**2)
      ENDDO
C      WRITE(*,*)'AV AND SD ', AV(1),AV(2),SD(1),SD(2)
      RHO=0.D0
      DO I=1,NOB
      RHO=RHO+(Z(I,1)-AV(1))*(Z(I,2)-AV(2))
      ENDDO
      RHO=(RHO/NOB)/(SD(1)*SD(2))
      RETURN
      END
C     -----------------------------------------------------------------
      SUBROUTINE CORA(Z,N,R)
C     COMPUTING BRADLEY'S ABSOLUTE CORRELATION MATRIX
C     BRADLEY, C. (1985) "THE ABSOLUTE CORRELATION", THE MATHEMATICAL
C     GAZETTE, 69(447): 12-17.
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION Z(N,2),X(N),Y(N)
C     ------------------------------------------------------------------
C     PUT Z INTO X AND Y
      DO I=1,N
      X(I)=Z(I,1)
      Y(I)=Z(I,2)
      ENDDO
C     ARRANGE X ANY IN AN ASCENDING ORDER
      DO I=1,N-1
      DO II=I+1,N
      IF(X(I).GT.X(II)) THEN
      TEMP=X(I)
      X(I)=X(II)
      X(II)=TEMP
      ENDIF
      IF(Y(I).GT.Y(II)) THEN
      TEMP=Y(I)
      Y(I)=Y(II)
      Y(II)=TEMP
      ENDIF
      ENDDO
      ENDDO
C     FIND MEDIAN
      IF(INT(N/2).EQ.N/2.D0) THEN
      XMED=(X(N/2)+X(N/2+1))/2.D0
      YMED=(Y(N/2)+Y(N/2+1))/2.D0
      ENDIF
      IF(INT(N/2).NE.N/2.D0) THEN
      XMED=X(N/2+1)
      YMED=Y(N/2+1)
      ENDIF
C     SUBTRACT RESPECTIVE MEDIANS FROM X AND Y AND FIND ABS DEVIATIONS
      VX=0.D0
      VY=0.D0
      DO I=1,N
      X(I)=X(I)-XMED
      Y(I)=Y(I)-YMED
      VX=VX+DABS(X(I))
      VY=VY+DABS(Y(I))
      ENDDO
C     SCALE THE VARIABLES X AND Y SUCH THAT VX=VY
      IF(VX.EQ.0.D0.OR.VY.EQ.0.D0) THEN
      R=0.D0
      RETURN
      ENDIF
      DO I=1,N
      X(I)=X(I)*VY/VX
      ENDDO
C     COMPUTE CORRELATION COEFFICIENT
      VZ=0.D0
      R=0.D0
      DO I=1,N
      VZ=VZ+DABS(X(I))+DABS(Y(I))
      R=R+DABS(X(I)+Y(I))-DABS(X(I)-Y(I))
      ENDDO
      R=R/VZ
      RETURN
      END
C     ------------------------------------------------------------------
      SUBROUTINE DORANK(X,N)! N IS THE NUMBER OF OBSERVATIONS
      PARAMETER (NRL=0) ! THIS VALUE IS TO BE SET BY THE USER
C                      !THE VALUE OF NRL DECIDES THE SCHEME OF RANKINGS
C     !THIS PROGRAM RETURNS RANK-ORDER OF A GIVEN VECTOR
      PARAMETER (MXD=1000)! MXD IS MAX DIMENSION FOR TEMPORARY VARIABLES
      ! THAT ARE LOCAL AND DO NOT GO TO THE INVOKING PROGRAM
      ! X IS THE VARIABLE TO BE SUBSTITUTED BY ITS RANK VALUES
C     NRULE=0 FOR ORDINAL RANKING (1-2-3-4 RULE);
C     NRULE=1 FOR DENSE RANKING (1-2-2-3 RULE);
C     NRULE=2 FOR STANDARD COMPETITION RANKING (1-2-2-4 RULE);
C     NRULE=3 FOR MODIFIED COMPETITION RANKING (1-3-3-4 RULE);
C     NRULE=4 FOR FRACTIONAL RANKING (1-2.5-2.5-4 RULE);
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION X(N),NF(MXD),NCF(MXD),RANK(MXD),ID(MXD),XX(MXD)
C     GENERATE ID(I),I=1,N
      DO I=1,N
      ID(I)=I
      NF(I)=0
      ENDDO
C     ARRANGE DATA (X) AND THE IDS IN ASCENDING ORDER
      DO I=1,N-1
      DO II=I,N
      IF(X(II).LT.X(I)) THEN
      TEMP=X(I)
      X(I)=X(II)
      X(II)=TEMP
      ITEMP=ID(I)
      ID(I)=ID(II)
      ID(II)=ITEMP
      ENDIF
      ENDDO
      ENDDO
C     MAKE DISCRETE UNGROUPED FREQUENCY TABLE
      K=0
      J=1
    1 K=K+1
      XX(K)=X(J)
      NF(K)=0
      DO I=J,N
      IF(XX(K).EQ.X(I)) THEN
      NF(K)=NF(K)+1
      ELSE
      J=I
      IF(J.LE.N) THEN
      GOTO 1
      ELSE
      GOTO 2
      ENDIF
      ENDIF
      ENDDO
    2 KK=K
      DO K=1,KK
      IF(K.EQ.1) THEN
      NCF(K)=NF(K)
      ELSE
      NCF(K)=NCF(K-1)+NF(K)
      ENDIF
      ENDDO
      DO I=1,N
      RANK(I)=1.D0
      ENDDO

      IF(NRL.GT.4) THEN
      WRITE(*,*)'RANKING RULE CODE GREATER THAN 4 NOT PERMITTED',NRL
      STOP
      ENDIF

      IF(NRL.LT.0) THEN
      WRITE(*,*)'RANKING RULE CODE LESS THAN 0 NOT PERMITTED',NRL
      STOP
      ENDIF

      IF(NRL.EQ.0) THEN
      DO I=1,N
      RANK(I)=I
      ENDDO
      ENDIF
C     ------------------------------------------------------------------
      IF(NRL.GT.0) THEN
       DO K=1,KK
       IF(K.EQ.1) THEN
       K1=1
       ELSE
       K1=NCF(K-1)+1
       ENDIF
       K2=NCF(K)
       DO I=K1,K2
       SUM=0.D0
       DO II=K1,K2
       SUM=SUM+II
       ENDDO
       KX=(K2-K1+1)
       IF(NRL.EQ.1)RANK(I)=K ! DENSE RANKING (1-2-2-3 RULE)
       IF(NRL.EQ.2)RANK(I)=K1!STANDARD COMPETITION RANKING(1-2-2-4 RULE)
       IF(NRL.EQ.3)RANK(I)=K2!MODIFIED COMPETITION RANKING(1-3-3-4 RULE)
       IF(NRL.EQ.4)RANK(I)=SUM/KX !FRACTIONAL RANKING (1-2.5-2.5-4 RULE)
       ENDDO
       ENDDO
      ENDIF
C     ------------------------------------------------------------------
      DO I=1,N
      X(ID(I))=RANK(I) ! BRINGS THE DATA TO ORIGINAL SEQUENCE
      ENDDO
      RETURN
      END
C     ----------------------------------------------------------------
      SUBROUTINE CORREL(X,N,M,RMAT)
      PARAMETER (NMX=30)!DO NOT CHANGE UNLESS NO. OF VARIABLES EXCEED 30
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION X(N,M),RMAT(2,2),AV(NMX),SD(NMX)
      DO J=1,2
      AV(J)=0.D0
      SD(J)=0.D0
      DO I=1,N
      AV(J)=AV(J)+X(I,J)
      SD(J)=SD(J)+X(I,J)**2
      ENDDO
      AV(J)=AV(J)/N
      SD(J)=DSQRT(SD(J)/N-AV(J)**2)
      ENDDO
      DO J=1,2
      DO JJ=1,2
      RMAT(J,JJ)=0.D0
      DO I=1,N
      RMAT(J,JJ)=RMAT(J,JJ)+X(I,J)*X(I,JJ)
      ENDDO
      ENDDO
      ENDDO
      DO J=1,2
      DO JJ=1,2
      RMAT(J,JJ)=RMAT(J,JJ)/N-AV(J)*AV(JJ)
      RMAT(J,JJ)=RMAT(J,JJ)/(SD(J)*SD(JJ))
      ENDDO
      ENDDO
      RETURN
      END
C     ------------------------------------------------------------------
      SUBROUTINE DOCORA(ZDAT,N,M,RMAT)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION ZDAT(N,M),RMAT(2,2),Z(N,2)
      DO I=1,N
      Z(I,1)=ZDAT(I,1)
      Z(I,2)=ZDAT(I,2)
      ENDDO
      CALL CORA(Z,N,R)
      RMAT(1,2)=R
      RMAT(2,1)=R
      DO J=1,2
      RMAT(J,J)=1.D0
      ENDDO
      RETURN
      END