C     =================================================================
C     PROGRAM FOR COMPUTING NEAREST CORRELATION MATRIX BY REPULSIVE
C     PARTICLE SWARM METHOD OF GLOBAL OPTIMIZATION  WRITTEN BY:
C       PROF. SK MISHRA, DEPARTMENT OF ECONOMICS, NEHU, SHILLONG INDIA
C                     APRIL 2007 : AVAILABLE AT:
C                   http//www.webng.com/economics
C                                  AND
C                 http://freewebs.com/nehu_economics
C     =================================================================
      PARAMETER (MX=10)! PRESENTLY SET AT 10. ONE MAY CHANGE IT.
C      ACCORDINGLY IT HAS TO BE CHANGED IN OTHER SUBROUTINES ALSO.
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /RNDM/IU,IV
      COMMON /KFF/KF,NFCALL,Q
      COMMON /EIGENS/VO(MX,MX),WO(MX,MX),WEIGHT(MX,MX),PPROD,PROD,NORM
      CHARACTER *70 FTIT,INFIL,INWFIL,OUTFIL
      DIMENSION Q(MX,MX),QQ(MX,MX),V(MX,MX),W(MX,MX),P(MX),R(MX,MX)
      CALL FSELECT(KF,M,FTIT)
C     READ THE Q (INVALID CORRELATION MATRIX) FROM INPUT FILE
      WRITE(*,*)'FILE THAT STORES Q (INVALID CORRELATION MATRIX) ?'
      READ(*,*) INFIL
      OPEN(7,FILE=INFIL)
      DO I=1,M
      READ(7,*) (Q(I,J),J=1,M) ! Q IS THE INVALID CORRELATION MATRIX
      DO J=1,M
      QQ(I,J)=Q(I,J) ! Q IS STORED IN QQ TO PRESERVE THE ORIGINAL MATRIX
      ENDDO
      ENDDO
      CLOSE(7)
      WRITE(*,*)'ANY WEIGHT MATRIX ? YES (1/ NON-ZERO), NO (0)'
      READ(*,*) NW
      IF(NW.NE.0) THEN
      WRITE(*,*)'FILE IN WHICH WEIGHT MATRIX HAS BEEN STORED ?'
      READ(*,*) INWFIL
      OPEN(7, FILE=INWFIL)
      DO I=1,M
      READ(7,*)(WEIGHT(I,J),J=1,M)
      ENDDO
      CLOSE(7)
      ELSE
      DO I=1,M
      DO J=1,M
      WEIGHT(I,J)=1.D0
      ENDDO
      ENDDO
      ENDIF
      WRITE(*,*)'CHOOSE THE NORM [ABS(1),FROBENIUS(2),CHEBYSHEV(=>99)]?'
      READ(*,*) NORM
      WRITE(*,*)'CHOOSE A SMALL(SUGGESTED 0.01)GUARD DETERMINANT[>ZERO]'
      WRITE(*,*)'[THIS GUARD DOES NOT ALLOW THE RESULTING NEAREST MATRIX
     & TO BECOME NEARSINGULAR]'
      READ(*,*) PPROD
      WRITE(*,*) 'GIVE AN OUTPUT FILE NAME TO STORE RESULTS'
      READ(*,*) OUTFIL
      OPEN(11,FILE=OUTFIL)! OPENS OUTPUT FILE
C     FIND EIGENVALUES AND EIGENVECTORS OF Q MATRIX
      CALL EIGEN(Q,VO,WO,M)
C-----------------------------------------------------------------------
      WRITE(*,*)'ORIGINAL MATRIX'
      WRITE(11,*)'ORIGINAL MATRIX'
      DO I=1,M
      WRITE(*,1)(Q(I,J),J=1,M)
      WRITE(11,1)(Q(I,J),J=1,M)
      ENDDO
    1 FORMAT(15F8.4) ! OUTPUT FORMAT. ONE MAY CHANGE IT AS NEEDED.
      WRITE(*,*)'EIGENVALUES'
      WRITE(11,*)'EIGENVALUES'
      DO I=1,M
      WRITE(*,2)(WO(I,J),J=1,M)
      WRITE(11,2)(WO(I,J),J=1,M)
      ENDDO
    2 FORMAT(10F12.7)! OUTPUT FORMAT. ONE MAY CHANGE IT AS NEEDED.
C-----------------------------------------------------------------------
C     CHECK IF ANY EIGENVALUE IS NEGATIVE
      LNEG=0
      DO I=1,M
      IF(WO(I,I).LT.0) LNEG=1 ! IF ANY EIGENVALUE NEGETIVE LNEG=1
      ENDDO
      IF(LNEG.EQ.1) THEN
      CALL RPS(M,P,FMIN,FTIT)! P IS THE VECTOR OF EIGENVALUES
      ENDIF
      CALL EIGEN(QQ,V,W,M)
C-----------------------------------------------------------------------
      LNEG=0
      DO I=1,M
      IF(W(I,I).LT.0) LNEG=1
      ENDDO
      DO I=1,M
      W(I,I)=P(I)
      ENDDO
      DO  J=1,M
      DO  JJ=1,M
      R(J,JJ)=0.0
      DO  I=1,M
      R(J,JJ)=R(J,JJ)+V(J,I)*W(I,JJ)
      ENDDO
      ENDDO
      ENDDO
      DO  J=1,M
      DO  JJ=1,M
      W(J,JJ)=0.0
      DO  I=1,M
      W(J,JJ)=W(J,JJ)+R(J,I)*V(JJ,I)
      ENDDO
      ENDDO
      ENDDO
      DO  I=1,M
      DO  J=1,M
      R(I,J)=W(I,J)
      ENDDO
      R(I,I)=1.0
      ENDDO
      CALL EIGEN(R,V,W,M)
      F=0
      DO I=1,M
      DO J=1,M
      DIF=DABS(Q(I,J)-R(I,J))
      IF(NORM.LT.99) THEN
      F=F+WEIGHT(I,J)*DIF**NORM
      ELSE
      IF(DIF*WEIGHT(I,J).GT.F) F=DIF*WEIGHT(I,J)
      ENDIF
      ENDDO
      ENDDO
      IF(LNEG.NE.0) THEN
      WRITE(*,*)' '
      WRITE(*,*)'-----------------------------------------------------'
      WRITE(*,*)'ESTIMATED MATRIX ------------------------------------'
      WRITE(*,*)' '
      WRITE(11,*)' '
      WRITE(11,*)'-----------------------------------------------------'
      WRITE(11,*)'ESTIMATED MATRIX ------------------------------------'
      WRITE(11,*)' '
      DO I=1,M
      WRITE(*,3)(R(I,J),J=1,M)
      WRITE(11,3)(R(I,J),J=1,M)
      ENDDO
    3 FORMAT(12F12.8) ! OUTPUT FORMAT. ONE MAY CHANGE IT AS NEEDED.
      WRITE(*,*)' '
      WRITE(*,*)'MIN NORM WITH PENALTIES (NOT TO BE INTERPRETED)=',F
      WRITE(11,*)' '
      WRITE(11,*)'MIN NORM WITH PENALTIES (NOT TO BE INTERPRETED)=',F
      SUM=0
      DET=1.D0
      DO I=1,M
      SUM=SUM+W(I,I)
      DET=DET*W(I,I)
      ENDDO
      WRITE(*,*)' '
      WRITE(*,*)'EIGENVALUES'
      WRITE(*,2)(W(I,I),I=1,M)
      WRITE(*,*)' '
      WRITE(11,*)' '
      WRITE(11,*)'EIGENVALUES'
      WRITE(11,2)(W(I,I),I=1,M)
      WRITE(11,*)' '
      WRITE(*,*)'SUM AND PRODUCT OF EIGENVALUES=',SUM,DET
      WRITE(11,*)'SUM AND PRODUCT OF EIGENVALUES=',SUM,DET
      ENDIF
      IF(LNEG.EQ.0) THEN
      WRITE(*,*)'THE GIVEN MATRIX IS VALID (POSITIVE SEMI-DEFINITE).'
      WRITE(11,*)'THE GIVEN MATRIX IS VALID (POSITIVE SEMI-DEFINITE).'
      ELSE
      WRITE(*,*)'THE RESULT MATRIX IS VALID (POSITIVE SEMI-DEFINITE).'
      WRITE(11,*)'THE RESULT MATRIX IS VALID (POSITIVE SEMI-DEFINITE).'
      ENDIF
      CLOSE(11)
      END
C     -----------------------------------------------------------------
      SUBROUTINE FUNC(X,M,F)
C     FINDING THE NEAREST CORRELATION MATRIX BY DIFFERENTIAL
C     EVOLUTION METHOD OF GLOBAL OPTIMIZATION  WRITTEN BY
C     PROF. SK MISHRA, DEPARTMENT OF ECONOMICS, NEHU, SHILLONG INDIA
      PARAMETER (MX=10)! PRESENTLY SET AT 10. ONE MAY CHANGE IT.
C      ACCORDINGLY IT HAS TO BE CHANGED IN MAIN/OTHER SUBROUTINES ALSO.
C      THIS SUBROUTNE IS USED BY THE GLOBAL OPTIMIZATION PROGRAM
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /RNDM/IU,IV
      COMMON /KFF/KF,NFCALL,Q
      COMMON /EIGENS/VO(MX,MX),WO(MX,MX),WEIGHT(MX,MX),PPROD,PROD,NORM
      DIMENSION X(*),Q(MX,MX),W(MX,MX),V(MX,MX),R(MX,MX)
      NFCALL=NFCALL+1 ! INCREMENT TO NUMBER OF FUNCTION CALLS
C     KF IS THE CODE OF THE TEST FUNCTION
C     -----------------------------------------------------------------
      LNEG=0
      DO I=1,M
      IF(WO(I,I).LT.0) LNEG=1
      ENDDO
      IF(LNEG.EQ.1) THEN
      SUM=0
      PROD=1.D0
      DO I=1,M
      IF(X(I).LT.0) THEN
      CALL RANDOM(RAND)
      X(I)=RAND
      ENDIF
      SUM=SUM+X(I)
      PROD=PROD*X(I)
      ENDDO
      DO I=1,M
      X(I)=X(I)*M/SUM
      ENDDO
      DO I=1,M
      W(I,I)=X(I)
      ENDDO

      DO  J=1,M
      DO  JJ=1,M
      R(J,JJ)=0.0
      DO  I=1,M
      R(J,JJ)=R(J,JJ)+VO(J,I)*W(I,JJ)
      ENDDO
      ENDDO
      ENDDO
      DO  J=1,M
      DO  JJ=1,M
      W(J,JJ)=0.0
      DO  I=1,M
      W(J,JJ)=W(J,JJ)+R(J,I)*VO(JJ,I)
      ENDDO
      ENDDO
      ENDDO
      DO  I=1,M
      DO  J=1,M
      R(I,J)=W(I,J)
      ENDDO
      R(I,I)=1.0
      ENDDO

      CALL EIGEN(R,V,W,M)
      PF=0
      SUM=0
      DO I=1,M
      SUM=SUM+W(I,I)
      IF(W(I,I).LT.0) PF=PF+100+(100*(1-W(I,I)))**2
      ENDDO
      F=0
      DO I=1,M
      F=F+100*(1.D0-R(I,I))**2
      ENDDO
      FN=0.0
      FX=0
      DO I=1,M-1
      DO J=I+1,M
      DIF=DABS(Q(I,J)-R(I,J))
      WDIF=DIF*WEIGHT(I,J)
      IF(WDIF.GT.FX) FX=WDIF
      IF(NORM.LT.99)FN=FN+WEIGHT(I,J)*DIF**NORM
      ENDDO
      ENDDO
      FNN=2*FN
      FN=0
      DO I=1,M
      DIF=DABS(Q(I,I)-R(I,I))
      WDIF=DIF*WEIGHT(I,I)
      IF(WDIF.GT.FX) FX=WDIF
      IF(NORM.LT.99)FN=FN+WEIGHT(I,I)*DIF**NORM
      ENDDO

         IF(NORM.LT.99) THEN
         F=F+FNN+FN+PF+DABS(SUM-M)*100
         ELSE
         F=F+FX+PF+DABS(SUM-M)*100
         ENDIF
      ENDIF
      IF(PROD.LT.PPROD) F=F+(1000*(PPROD-PROD))**2

      RETURN
      END
C     =================================================================
C     ==============  REPULSIVE PARTICLE SWARM METHOD  ================
      SUBROUTINE RPS(M,BST,FMINIM,FTIT)
C     PROGRAM TO FIND GLOBAL MINIMUM BY REPULSIVE PARTICLE SWARM METHOD
C     WRITTEN BY SK MISHRA, DEPT. OF ECONOMICS, NEHU, SHILLONG (INDIA)
C     -----------------------------------------------------------------
      PARAMETER (MZ=10)! THIS CORRESPONDS TO MX IN THE MAIN PROGRAM
C     NOT TO BE CONFUSED WITH THE MX THAT APPEARS IN THE PARAMETER BELOW
C     PRESENTLY SET AT 10. ONE MAY CHANGE IT.
C      ACCORDINGLY IT HAS TO BE CHANGED IN MAIN/OTHER SUBROUTINES ALSO.
C     -----------------------------------------------------------------
      PARAMETER (N=100,NN=50,MX=50,NSTEP=11,ITRN=10000,NSIGMA=1,ITOP=3)
      PARAMETER (NPRN=500) ! ECHOS 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,Q
      COMMON /EIGENS/VO(MZ,MZ),WO(MZ,MZ),WEIGHT(MZ,MZ),PPROD,PROD,NORM
      INTEGER IU,IV
      CHARACTER *70 FTIT
      DIMENSION X(N,MX),V(N,MX),A(MX),VI(MX),Q(MZ,MZ)
      DIMENSION XX(N,MX),F(N),V1(MX),V2(MX),V3(MX),V4(MX),BST(MZ)
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,EPSI /.5D0,.5D0,5.D-04,.5D00,1.D-03,1.D-08/
C     -----------------------------------------------------------------
C     -----------------------------------------------------------------
      GGBEST=1.D30 ! TO BE USED FOR TERMINATION CRITERION
      LCOUNT=0
      NFCALL=0
      WRITE(*,*)'4-DIGITS SEED FOR RANDOM NUMBER GENERATION'
      WRITE(*,*)'A FOUR-DIGIT POSITIVE ODD INTEGER, SAY, 1171'
      READ(*,*) IU
      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-0.5D00)*2000
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     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 ******************************************
       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 **************************************************
        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
      IF(DABS(FMIN-GGBEST).LT.EPSI) THEN
           WRITE(*,*)'COMPUTATION OVER'
           FMINIM=FMIN
           RETURN
           ELSE
           GGBEST=FMIN
           ENDIF
      ENDIF
      LCOUNT=LCOUNT+1
  100 CONTINUE
      WRITE(*,*)'COMPUTATION OVER:',FTIT
      FMINIM=FMIN
      RETURN
      END
C     ----------------------------------------------------------------
      SUBROUTINE LSRCH(A,M,VI,NSTEP,FI)
      PARAMETER (MX=10)! PRESENTLY SET AT 10. ONE MAY CHANGE IT.
C      ACCORDINGLY IT HAS TO BE CHANGED IN MAIN/OTHER SUBROUTINES ALSO.
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /KFF/KF,NFCALL,Q
      COMMON /RNDM/IU,IV
      DIMENSION A(*),B(100),VI(*),Q(MX,MX) ! Q NOT USED HERE
      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     =================================================================
      SUBROUTINE NORMAL(R)
C     PROGRAM TO GENERATE N(0,1) FROM RECTANGULAR RANDOM NUMBERS
C     IT USES VARIATE TRANSFORMATION FOR THIS PURPOSE.
C     -----------------------------------------------------------------
C     IF U1 AND U2 ARE UNIFORMLY DISTRIBUTED RANDOM NUMBERS (0,1),
C     THEN X=[(-2*LN(U1))**.5]*(COS(2*PI*U2) IS N(0,1)
C     PI = 4*ARCTAN(1.0)= 3.1415926535897932384626433832795
C     2*PI = 6.283185307179586476925286766559
C     -----------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /RNDM/IU,IV
C     -----------------------------------------------------------------
      CALL RANDOM(RAND) ! INVOKES RANDOM TO GENERATE UNIFORM RAND [0, 1]
      U1=RAND ! U1 IS UNIFORMLY DISTRIBUTED [0, 1]
      CALL RANDOM(RAND) ! INVOKES RANDOM TO GENERATE UNIFORM RAND [0, 1]
      U2=RAND ! U1 IS UNIFORMLY DISTRIBUTED [0, 1]
      R=DSQRT(-2.D0*DLOG(U1))
C      R=R*DCOS(U2*6.283185307179586476925286766559D00)
      R=R*DCOS(U2*6.28318530718D00)
      RETURN
      END
C     ------------------------------------------------------------------
C     RANDOM NUMBER GENERATOR (UNIFORM BETWEEN 0 AND 1 - BOTH EXCLUSIVE)
      SUBROUTINE RANDOM(RAND)
       IMPLICIT DOUBLE PRECISION  (A-H, O-Z)
       COMMON /RNDM/IU,IV
       IV=IU*65539
       IF(IV.LT.0) THEN
       IV=IV+2147483647+1
       ENDIF
       RAND=IV
       IU=IV
       RAND=RAND*0.4656613E-09
       RETURN
       END
C     -----------------------------------------------------------------
      SUBROUTINE FSELECT(KF,M,FTIT)
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(1),FTIT
      WRITE(*,*)'----------------------------------------------------'
      DATA TIT(1)/'NEAREST CORRELATION MATRIX PROBLEM'/
C     -----------------------------------------------------------------
      KF=1
      WRITE(*,*)'      THE NEAREST CORRELATION MATRIX PROBLEM'
      WRITE(*,*)'THE METHOD AND PROGRAM DEVELOPED BY PROF. SK MISHRA'
      WRITE(*,*)' DEPT. OF ECONOMICS, NORTH-EASTERN HILL UNIVERSITY'
      WRITE(*,*)'                SHILLONG (INDIA)'
      WRITE(*,*)'----------------------------------------------------'
      WRITE(*,*)' '
      WRITE(*,*)'PLEASE FEED THE INFORMATION AS ASKED FOR'
      WRITE(*,*)' '
      WRITE(*,*)'DIMENSION OF THE INVALID CORRELATION MATRIX [M] ?'
      READ(*,*) M
      FTIT=TIT(KF) ! STORE THE NAME OF THE CHOSEN FUNCTION IN FTIT
      RETURN
      END
C ----------------------------------------------------
      SUBROUTINE EIGEN(A,V,W,N)
      PARAMETER (MX=10)! PRESENTLY SET AT 10. ONE MAY CHANGE IT.
C      ACCORDINGLY IT HAS TO BE CHANGED IN MAIN/OTHER SUBROUTINES ALSO.
C     KRISHNAMURTHY, EV AND SK SEN (1976). "COMPUTER-BASED NUMERICAL
C     ALGORITHMS", AFFILIATED EAST-WEST PRESS, NEW DELHI.
C     COMPUTES EIGENVALUES AND VECTORS OF A REAL SYMMETRIX MATRIX
      DOUBLE PRECISION A(MX,MX),V(MX,MX),W(MX,MX),P(MX)
      DOUBLE PRECISION PMAX,EPLN,TAN,SIN,COS,AI,TT,TA,TB
      DIMENSION MM(MX)
C ------------ INITIALISATION -------------------------
      DO 50 I=1,N
      DO 51 J=1,N
      V(I,J)=0.0
51    W(I,J)=A(I,J)
      P(I)=0.0
50    CONTINUE
      PMAX=0
      EPLN=0
      TAN=0
      SIN=0
      COS=0
      AI=0
      TT=0
      NN=1
      EPLN=1.0D-99
C ------------------------------------------------------
      IF(NN.NE.0) THEN
      DO 3 I=1,N
      DO 3 J=1,N
      V(I,J)=0.0
      IF(I.EQ.J) V(I,J)=1.0
3     CONTINUE
      ENDIF
2     NR=0
5     MI=N-1
      DO 6 I=1,MI
      P(I)=0.0
      MJ=I+1
      DO 6 J=MJ,N
      IF(P(I).GT.DABS(A(I,J))) GO TO 6
      P(I)=DABS(A(I,J))
      MM(I)=J
6      CONTINUE
7       DO 8 I=1,MI
        IF(I.LE.1) GOTO 10
        IF(PMAX.GT.P(I)) GOTO 8
10       PMAX=P(I)
         IP=I
         JP=MM(I)
8     CONTINUE
      IF (PMAX.LE.EPLN) THEN
      GO TO 12
      ENDIF
      NR=NR+1
13    TA=2.0*A(IP,JP)
      TB=(DABS(A(IP,IP)-A(JP,JP))+
     &DSQRT((A(IP,IP)-A(JP,JP))**2+4.0*A(IP,JP)**2))
      TAN=TA/TB
      IF(A(IP,IP).LT.A(JP,JP)) TAN=-TAN
14    COS=1.0/DSQRT(1.0+TAN**2)
      SIN=TAN*COS
      AI=A(IP,IP)
      A(IP,IP)=(COS**2)*(AI+TAN*(2.0*A(IP,JP)+TAN*A(JP,JP)))
      A(JP,JP)=(COS**2)*(A(JP,JP)-TAN*(2.0*A(IP,JP)-TAN*AI))
      A(IP,JP)=0.0
      IF(A(IP,IP).GE.A(JP,JP)) GO TO 15
      TT=A(IP,IP)
      A(IP,IP)=A(JP,JP)
      A(JP,JP)=TT
      IF(SIN.GE.0) GO TO 16
      TT=COS
      GO TO 17
16    TT=-COS
17    COS=DABS(SIN)
      SIN=TT
15    DO 18 I=1,MI
      IF(I-IP) 19, 18, 20
20    IF(I.EQ.JP)GO TO 18
19    IF(MM(I).EQ.IP) GO TO 21
      IF(MM(I).NE.JP) GO TO 18
21    K=MM(I)
      TT=A(I,K)
      A(I,K)=0.0
      MJ=I+1
      P(I)=0.0
      DO 22 J=MJ,N
      IF(P(I).GT.DABS(A(I,J))) GO TO 22
      P(I)=DABS(A(I,J))
      MM(I)=J
22    CONTINUE
      A(I,K)=TT
18    CONTINUE
      P(IP)=0.0
      P(JP)=0.0
      DO 23 I=1,N
      IF(I-IP) 24, 23, 25
24    TT=A(I,IP)
      A(I,IP)=COS*TT+SIN*A(I,JP)
      IF(P(I).GE.DABS(A(I,IP))) GO TO 26
      P(I)=DABS(A(I,IP))
      MM(I)=IP
26    A(I,JP)=-SIN*TT+COS*A(I,JP)
      IF(P(I).GE.DABS(A(I,JP))) GO TO 23
30    P(I)=DABS(A(I,JP))
      MM(I)=JP
      GO TO 23
25    IF(I.LT.JP) GO TO 27
      IF(I.GT.JP) GO TO 28
      IF(I.EQ.JP) GO TO 23
27    TT=A(IP,I)
      A(IP,I)=COS*TT+SIN*A(I,JP)
      IF(P(IP).GE.DABS(A(IP,I))) GO TO 29
      P(IP)=DABS(A(IP,I))
      MM(IP)=I
29    A(I,JP)=-TT*SIN+COS*A(I,JP)
      IF(P(I).GE.DABS(A(I,JP))) GO TO 23
      GO TO 30
28    TT=A(IP,I)
      A(IP,I)=TT*COS+SIN*A(JP,I)
      IF(P(IP).GE.DABS(A(IP,I))) GO TO 31
      P(IP)=DABS(A(IP,I))
      MM(IP)=I
31    A(JP,I)=-TT*SIN+COS*A(JP,I)
      IF(P(JP).GE.DABS(A(JP,I))) GO TO 23
      P(JP)=DABS(A(JP,I))
      MM(JP)=I
23    CONTINUE
      IF(NN.EQ.0) GOTO 7
      DO 32 I=1,N
      TT=V(I,IP)
      V(I,IP)=TT*COS+SIN*V(I,JP)
      V(I,JP)=-TT*SIN+COS*V(I,JP)
32    CONTINUE
      GO TO 7
12    DO I=1,N
      DO J=1,N
      TEMP=A(I,J)
      A(I,J)=W(I,J)
      W(I,J)=TEMP
      ENDDO
      ENDDO
      DO I=1,N
      DO J=1,N
      IF(I.NE.J) W(I,J)=0.0
      ENDDO
      ENDDO
C     NOW A STORES THE ORIGINAL MATRIX, W(I,I) STORES WIGENVALUES AND
C     V STORES EIGENVECTORS
      RETURN
      END