C     MAIN PROGRAM : GENERATE A SEMIPOSITIVE CORRELATION MATRIX FROM
C     A GIVEN CORRELATION MATRIX WITH SOME KNOWN ELEMENTS
C     -----------------------------------------------------------------
C                      METHOD:DIFFERENTIAL EVOLUTION
C     -----------------------------------------------------------------
C     ADJUST THE PARAMETERS SUITABLY IN SUBROUTINES DE
C     WHEN THE PROGRAM ASKS FOR PARAMETERS, FEED THEM SUITABLY
C     ================ MAIN PROGRAM ===================================

      IMPLICIT DOUBLE PRECISION (A-H, O-Z)
      CHARACTER *70 INFILE,OUTFIL
      COMMON /XBASE/XBAS
      COMMON /RNDM/IU,IV ! RANDOM NUMBER GENERATION (IU = 4-DIGIT SEED)
      COMMON /NEAREST/Z,MORDER
      DIMENSION Z(10,10) ! THE INPUT CORRELATION MATRIX
      INTEGER IU,IV
      DIMENSION XBAS(500,50)
      DIMENSION X(50)! X IS THE DECISION VARIABLE X IN F(X) TO MINIMIZE
C     M IS THE DIMENSION OF THE PROBLEM, KF IS TEST FUNCTION CODE AND
C     FMIN IS THE MIN VALUE OF F(X) OBTAINED FROM DE
C     ------------------------------------------------------------------
      WRITE(*,*)'====== CONSTRUCTION OF VALID CORRELATION MATRIX ======'
      WRITE(*,*)'====== OPTIMIZATION BY DIFFERENTIAL EVOLUTION ========'
C     ------------------------------------------------------------------
      WRITE(*,*)'ORDER OF INPUT MATRIX (MORDER)& NAME OF INPUT FILE ?'
      READ(*,*) MORDER,INFILE
      WRITE(*,*)'SPECIFY THE OUTPUT FILE TO STORE VALID OUTPUT MATRICES'
      READ(*,*) OUTFIL
C     READ THE GIVEN CORRELATION MATRIX (UPPER DIAGONAL ONLY)
      OPEN(7,FILE=INFILE) ! OPEN INPUT FILE  AND READ THE MATRIX
      DO I=1,MORDER
      READ(7,*)(Z(I,J),J=I,MORDER)
C      write(*,*)(Z(I,J),J=I,MORDER)
      ENDDO
      CLOSE(7)
C     -----------------------------------------------------------------
C      WRITE(*,*)'====================     WARNING    =============== '
C      WRITE(*,*)'ADJUST PARAMETERS IN SUBROUTINES DE SUBROUTINE'
C      WRITE(*,*)'====================     WARNING    =============== '

C     INITIALIZATION. THIS XBAS WILL BE USED IN PROGRAMS TO
C     INITIALIZE THE POPULATION.
      WRITE(*,*)' '
      WRITE(*,*)'FEED RANDOM NUMBER SEED [4-DIGIT ODD INTEGER] TO BEGIN'
      READ(*,*) IU
C     THIS XBAS WILL BE USED IN ALL THE THREE METHODS AS INITIAL X
      DO I=1,500
      DO J=1,50
      CALL RANDOM(RAND)
      XBAS(I,J)=(RAND-0.5D00)*2 ! RANDOM NUMBER BETWEEN (-1, 1)
      ENDDO
      ENDDO
      WRITE(*,*)' *****************************************************'
C     ------------------------------------------------------------------

C      WRITE(*,*)'TO PROCEED TYPE ANY CHARACTER AND STRIKE ENTER'
C      READ(*,*) PROCEED
      CALL DE(M,X,FMINDE) ! CALLS DE AND RETURNS OPTIMAL X AND FMIN
C     ------------------------------------------------------------------
      CALL NCORX(X,M,OUTFIL)
      WRITE(*,*)'PROGRAM ENDED, FOR RESULTS OPEN OUTPUT FILE ',OUTFIL
      WRITE(*,*)'******************************************************'
      END
C     -----------------------------------------------------------------
      SUBROUTINE DE(M,A,FBEST)
C     PROGRAM: "DIFFERENTIAL EVOLUTION ALGORITHM" OF GLOBAL OPTIMIZATION
C     THIS METHOD WAS PROPOSED BY R. STORN AND K. PRICE IN 1995. REF --
C     "DIFFERENTIAL EVOLUTION - A SIMPLE AND EFFICIENT ADAPTIVE SCHEME
C     FOR GLOBAL OPTIMIZATION OVER CONTINUOUS SPACES" : TECHNICAL REPORT
C     INTERNATIONAL COMPUTER SCIENCE INSTITUTE, BERKLEY, 1995.
C     PROGRAM BY SK MISHRA, DEPT. OF ECONOMICS, NEHU, SHILLONG (INDIA)
C     -----------------------------------------------------------------
C     PROGRAM DE
      IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! TYPE DECLARATION
      PARAMETER(NMAX=500,MMAX=50) ! MAXIMUM DIMENSION PARAMETERS
      PARAMETER (RX1=1.d0, RX2=0.d0) ! TO BE ADJUSTED SUITABLY, IF NEEDED
C     RX1 AND RX2 CONTROL THE SCHEME OF CROSSOVER. (0 <= RX1 <= RX2) <=1
C     RX1 DETERMINES THE UPPER LIMIT OF SCHEME 1 (AND LOWER LIMIT OF
C     SCHEME 2; RX2 IS THE UPPER LIMIT OF SCHEME 2 AND LOWER LIMIT OF
C     SCHEME 3. THUS RX1 = .2 AND RX2 = .8 MEANS 0-20% SCHEME1, 20 TO 80
C     PERCENT SCHEME 2 AND THE REST (80 TO 100 %) SCHEME 3.
C     PARAMETER(NCROSS=2) ! CROSS-OVER SCHEME (NCROSS <=0 OR =1 OR =>2)
      PARAMETER(IPRINT=500,EPS=1.D-08)!FOR WATCHING INTERMEDIATE RESULTS
C     IT PRINTS THE INTERMEDIATE RESULTS AFTER EACH IPRINT ITERATION AND
C     EPS DETERMINES ACCURACY FOR TERMINATION. IF EPS= 0, ALL ITERATIONS
C     WOULD BE UNDERGONE EVEN IF NO IMPROVEMENT IN RESULTS IS THERE.
C     ULTIMATELY "DID NOT CONVERGE" IS REOPORTED.
      COMMON /RNDM/IU,IV ! RANDOM NUMBER GENERATION (IU = 4-DIGIT SEED)
      INTEGER IU,IV      ! FOR RANDOM NUMBER GENERATION
      COMMON /XBASE/XBAS
      common /nfcal/nfcall
      CHARACTER *70 FTIT ! TITLE OF THE FUNCTION
      CHARACTER *15 CFIL !OUTPUT FILE
C     -----------------------------------------------------------------
C     THE PROGRAM REQUIRES INPUTS FROM THE USER ON THE FOLLOWING ------
C     (1) FUNCTION CODE (KF), (2) NO. OF VARIABLES IN THE FUNCTION (M);
C     (3) N=POPULATION SIZE (SUGGESTED 10 TIMES OF NO. OF VARIABLES, M,
C         FOR SMALLER PROBLEMS N=100 WORKS VERY WELL);
C     (4) PCROS = PROB. OF CROSS-OVER (SUGGESTED : ABOUT 0.85 TO .99);
C     (5) FACT = SCALE (SUGGESTED 0.5 TO .95 OR  1, ETC);
C     (6) ITER = MAXIMUM NUMBER OF ITERATIONS PERMITTED (5000 OR MORE)
C     (7) RANDOM NUMBER SEED (4 DIGITS INTEGER)
C     ----------------------------------------------------------------
      DIMENSION X(NMAX,MMAX),Y(NMAX,MMAX),A(MMAX),FV(NMAX)
      DIMENSION IR(3),XBAS(500,50)
C     ----------------------------------------------------------------
C     ------- SELECT THE FUNCTION TO MINIMIZE AND ITS DIMENSION -------
      CALL FSELECT(KF,M,FTIT)
      CFIL='CORRESULTS'! IT IS AN INTERMEDIATE FILE
C     SPECIFY OTHER PARAMETERS ---------------------------------------
      WRITE(*,*)'POPULATION SIZE [N] AND NO. OF ITERATIONS [ITER] ?'
      WRITE(*,*)'SUGGESTED: N=>100 OR =>10.M; ITERATION 500 OR LARGER'
      READ(*,*) N,ITER
C      WRITE(*,*)'CROSSOVER PROBABILITY [PCROS] AND SCALE [FACT] ?'
C      WRITE(*,*)'SUGGESTED : PCROS ABOUT 0.9; FACT=.5 OR LARGER BUT <=1'
C      READ(*,*) PCROS,FACT
      PCROS=0.9d0
      FACT=0.5d0
      WRITE(*,*)'RANDOM NUMBER SEED ?'
      WRITE(*,*)'A FOUR-DIGIT POSITIVE ODD INTEGER, SAY, 1171'
      READ(*,*) IU

      NFCALL=0 ! INITIALIZE COUNTER FOR FUNCTION CALLS
      GBEST=1.D30 ! TO BE USED FOR TERMINATION CRITERION
C     INITIALIZATION : GENERATE X(N,M) RANDOMLY
      DO I=1,N
      DO J=1,M
C      CALL RANDOM(RAND) ! GENERATES INITION X WITHIN
C      X(I,J)=(RAND-.5D00)*2000 ! GENERATES INITION X WITHIN
C     RANDOM NUMBERS BETWEEN -RRANGE AND +RRANGE (BOTH EXCLUSIVE)
      X(I,J)=XBAS(I,J)! TAKES THESE NUMBERS FROM THE MAIN PROGRAM
      ENDDO
      ENDDO
      WRITE(*,*)'COMPUTING --- PLEASE WAIT '
      IPCOUNT=0
      DO 100 ITR=1,ITER  ! ITERATION BEGINS
C     -----------------------------------------------------------------
C     EVALUATE ALL X FOR THE GIVEN FUNCTION
      DO I=1,N
      DO J=1,M
      A(J)=X(I,J)
      ENDDO
      CALL FUNC(A,M,F)
C     STORE FUNCTION VALUES IN FV VECTOR
      FV(I)=F
      ENDDO

C     ----------------------------------------------------------------
C     FIND THE FITTEST (BEST) INDIVIDUAL AT THIS ITERATION
                FBEST=FV(1)
                KB=1
                DO IB=2,N
                     IF(FV(IB).LT.FBEST) THEN
                     FBEST=FV(IB)
                     KB=IB
                     ENDIF
                ENDDO
C     BEST FITNESS VALUE = FBEST : INDIVIDUAL X(KB)
C     -----------------------------------------------------------------
C     GENERATE OFFSPRINGS
      DO I=1,N    ! I LOOP BEGINS
C     INITIALIZE CHILDREN IDENTICAL TO PARENTS; THEY WILL CHANGE LATER
           DO J=1,M
           Y(I,J)=X(I,J)
           ENDDO
C     SELECT RANDOMLY THREE OTHER INDIVIDUALS
   20     DO IRI=1,3  ! IRI LOOP BEGINS
          IR(IRI)=0

          CALL RANDOM(RAND)
           IRJ=INT(RAND*N)+1
C     CHECK THAT THESE THREE INDIVIDUALS ARE DISTICT AND OTHER THAN I
           IF(IRI.EQ.1.AND.IRJ.NE.I) THEN
           IR(IRI)=IRJ
           ENDIF
           IF(IRI.EQ.2.AND.IRJ.NE.I.AND.IRJ.NE.IR(1)) THEN
           IR(IRI)=IRJ
           ENDIF
       IF(IRI.EQ.3.AND.IRJ.NE.I.AND.IRJ.NE.IR(1).AND.IRJ.NE.IR(2)) THEN
            IR(IRI)=IRJ
            ENDIF
            ENDDO   ! IRI LOOP ENDS
C     CHECK IF ALL THE THREE IR ARE POSITIVE (INTEGERS)
          DO IX=1,3
          IF(IR(IX).LE.0) THEN
          GOTO 20  ! IF NOT THEN REGENERATE
          ENDIF
          ENDDO
C     THREE RANDOMLY CHOSEN INDIVIDUALS DIFFERENT FROM I AND DIFFERENT
C     FROM EACH OTHER ARE IR(1),IR(2) AND IR(3)
C     ===================== RANDOMIZATION OF NCROSS ===================
C     RANDOMIZES NCROSS
      NCROSS=0
      CALL RANDOM(RAND)
      IF(RAND.GT.RX1) NCROSS=1 ! IF RX1=>1, SCHEME 2 NEVER IMPLEMENTED
      IF(RAND.GT.RX2) NCROSS=2 ! IF RX2=>1, SCHEME 3 NEVER IMPLEMENTED

C     ---------------------- SCHEME 1 ----------------------------------
C      NO CROSS OVER, ONLY REPLACEMENT THAT IS PROBABILISTIC
          IF(NCROSS.LE.0) THEN
          DO J=1,M      ! J LOOP BEGINS
          CALL RANDOM(RAND)
          IF(RAND.LE.PCROS) THEN ! REPLACE IF RAND < PCROS
          A(J)=X(IR(1),J)+(X(IR(2),J)-X(IR(3),J))*FACT ! CANDIDATE CHILD
          ENDIF
          ENDDO  ! J LOOP ENDS
          ENDIF

C     ----------------------- SCHEME 2 ---------------------------------
C     THE STANDARD CROSSOVER SCHEME
C     CROSSOVER SCHEME (EXPONENTIAL) SUGGESTED BY KENNETH PRICE IN HIS
C     PERSONAL LETTER TO THE AUTHOR (DATED SEPTEMBER 29, 2006)
      IF(NCROSS.EQ.1) THEN
        CALL RANDOM(RAND)
        JR=INT(RAND*M)+1
        J=JR
    2   A(J)=X(IR(1),J)+FACT*(X(IR(2),J)-X(IR(3),J))
        J=J+1
        IF(J.GT.M) J=1
        IF(J.EQ.JR) GOTO 10
        CALL RANDOM(RAND)
        IF(PCROS.LE.RAND) GOTO 2
    6   A(J)=X(I,J)
        J=J+1
        IF(J.GT.M) J=1
        IF (J.EQ.JR) GOTO 10
        GOTO 6
   10   CONTINUE
      ENDIF
C     ------------------------ SCHEME 3 --------------------------------
C     ESPECIALLY SUITABLE TO NON-DECOMPOSABLE (NON-SEPERABLE) FUNCTIONS
C     CROSSOVER SCHEME (NEW) SUGGESTED BY KENNETH PRICE IN HIS
C     PERSONAL LETTER TO THE AUTHOR (DATED OCTOBER 18, 2006)
      IF(NCROSS.GE.2) THEN
          CALL RANDOM(RAND)
          IF(RAND.LE.PCROS) THEN
             CALL NORMAL(RN)
             DO J=1,M
             A(J)=X(I,J)+(X(IR(1),J)+ X(IR(2),J)-2*X(I,J))*RN
             ENDDO
           ELSE
            DO J=1,M
            A(J)=X(I,J)+(X(IR(1),J)- X(IR(2),J))! FACT ASSUMED TO BE 1
            ENDDO
           ENDIF
      ENDIF
C     -----------------------------------------------------------------
          CALL FUNC(A,M,F) ! EVALUATE THE OFFSPRING
          IF(F.LT.FV(I)) THEN ! IF BETTER, REPLACE PARENTS BY THE CHILD
          FV(I)=F
          DO J=1,M
          Y(I,J)=A(J)
          ENDDO
          ENDIF
      ENDDO   ! I LOOP ENDS
      DO I=1,N
      DO J=1,M
      X(I,J)=Y(I,J) ! NEW GENERATION IS A MIX OF BETTER PARENTS AND
C                     BETTER CHILDREN
      ENDDO
      ENDDO
      IPCOUNT=IPCOUNT+1
      IF(IPCOUNT.EQ.IPRINT) THEN
      DO J=1,M
      A(J)=X(KB,J)
      ENDDO
      WRITE(*,*)(X(KB,J),J=1,M),'  FBEST UPTO NOW = ',FBEST
      WRITE(*,*)'TOTAL NUMBER OF FUNCTION CALLS =',NFCALL
           IF(DABS(FBEST-GBEST).LT.EPS) THEN
           WRITE(*,*) FTIT
           WRITE(*,*)'COMPUTATION OVER'
           GOTO 999
           ELSE
           GBEST=FBEST
           ENDIF
      IPCOUNT=0
      ENDIF
C     ----------------------------------------------------------------
  100 ENDDO   ! ITERATION ENDS : GO FOR NEXT ITERATION, IF APPLICABLE
C     ----------------------------------------------------------------
      WRITE(*,*)'DID NOT CONVERGE. REDUCE EPS OR RAISE ITER OR DO BOTH'
      WRITE(*,*)'INCREASE N, PCROS, OR SCALE FACTOR (FACT)'
  999 OPEN(7,FILE=CFIL)
      WRITE(7,*) N
      DO I=1,N
      WRITE(7,*)(X(I,J),J=1,M),FV(I)
      ENDDO
      CLOSE(7)
      RETURN
      END
C     -----------------------------------------------------------------
      SUBROUTINE NORMAL(R)
C     PROGRAM TO GENERATE N(0,1) FROM RECTANGULAR RANDOM NUMBERS
C     IT USES BOX-MULLER VARIATE TRANSFORMATION FOR THIS PURPOSE.
C     -----------------------------------------------------------------
C     ----- BOX-MULLER METHOD BY GEP BOX AND ME MULLER (1958) ---------
C     BOX, G. E. P. AND MULLER, M. E. "A NOTE ON THE GENERATION OF
C     RANDOM NORMAL DEVIATES." ANN. MATH. STAT. 29, 610-611, 1958.
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     ALSO,  X=[(-2*LN(U1))**.5]*(SIN(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
      INTEGER 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))
      R=R*DCOS(U2*6.283185307179586476925286766559D00)
C     R=R*DCOS(U2*6.28318530718D00)
      RETURN
      END
C     -----------------------------------------------------------------
C     RANDOM NUMBER GENERATOR (UNIFORM BETWEEN 0 AND 1 - BOTH EXCLUSIVE)
      SUBROUTINE RANDOM(RAND1)
       DOUBLE PRECISION  RAND1
       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
       RAND1= dble(RAND)
       RETURN
       END
C     -----------------------------------------------------------------C     -----------------------------------------------------------------
      SUBROUTINE FSELECT(KF,M,FTIT)
C     PARAMETER (NFUNCT=1) ! NO. OF FUNCTIONS IN THE LIST
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
C      WRITE(*,*)'----------------------------------------------------'
      KF=1 ! NO. OF FUNCTIONS
      DATA TIT(1)/'KF=1 CONSTRUCT CORRELATION MATRIX:M-VARIABLES M=?'/
C     -----------------------------------------------------------------
C      DO I=1,NFUNCT
C      WRITE(*,*)TIT(I)
C      ENDDO
C      WRITE(*,*)'*****************************************************'
      WRITE(*,*)'NO. OF VARIABLES=UNKNOWN CORRELATION COEFFICIENTS [M]?'
      READ(*,*) M
      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 /nfcal/nfcall
      DIMENSION X(*)
      NFCALL=NFCALL+1 ! INCREMENT TO NUMBER OF FUNCTION CALLS
C     KF IS THE CODE OF THE TEST FUNCTION
C     -----------------------------------------------------------------
C     CONSTRUCT CORRELATION MATRIX
      CALL CONCOR(M,X,F)
      RETURN
C     =================================================================
      END
C     ------------------------------------------------------------------
      SUBROUTINE EIGEN(A,N,V,W)
C     COMPUTES EIGENVALUES AND VECTORS OF A REAL SYMMETRIC MATRIX
C     A(N,N) =GIVEN REAL SYMMETRIC MATRIX WHOSE EIGENVALUES AND VECTORS
C     ARE BE FOUND. ITS ORDER IS N X N
C     W(N,N) CONTAINS EIGENVALUES IN ITS MAIN DIAGONAL. OTHER ELEMENTS=0
C     V(N,N) CONTAINS EIGENVECTORS
C     PROGRAM BY KRISNAMURTHY,EV & SEN (1976) COMPUTER-BASED NUMERICAL
C     ALGORITHMS. AFFILIATED EAST-WEST PRESS, NEW DELHI

      DOUBLE PRECISION A(10,10),V(10,10),W(10,10),P(10)
      DOUBLE PRECISION PMAX,EPLN,TAN,SIN,COS,AI,TT,TA,TB
      DIMENSION MM(10)
C     ------------ INITIALISATION --------------------------------------
      DO I=1,N
      DO J=1,N
      V(I,J)=0.d0
      W(I,J)=A(I,J)
      ENDDO
      P(I)=0.d0
      ENDDO
      PMAX=0.d0
      EPLN=0.d0
      TAN=0.d0
      SIN=0.d0
      COS=0.d0
      AI=0.d0
      TT=0
      NN=1
      EPLN=1.0D-100
C     ------------------------------------------------------------------
      IF(NN.NE.0) THEN
      DO I=1,N
      DO J=1,N
      V(I,J)=0.d0
      IF(I.EQ.J) V(I,J)=1.d0
      ENDDO
      ENDDO
      ENDIF
      NR=0
      MI=N-1
      DO I=1,MI
      P(I)=0.d0
      MJ=I+1
      DO J=MJ,N
      IF(P(I).LE.DABS(A(I,J))) THEN
      P(I)=DABS(A(I,J))
      MM(I)=J
      ENDIF
      ENDDO
      ENDDO

    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
C     EPLN=DABS(PMAX)*1.0D-09
      IF (PMAX.LE.EPLN) THEN
C     WRITE(*,*)'PMAX EPLN',PMAX, EPLN
C     PAUSE'CONVERGENCE CRITERION IS MET'
      GO TO 12
      ENDIF
      NR=NR+1
      TA=2.d0*A(IP,JP)
      TB=(DABS(A(IP,IP)-A(JP,JP))+
     1DSQRT((A(IP,IP)-A(JP,JP))**2+4.d0*A(IP,JP)**2))
      TAN=TA/TB
      IF(A(IP,IP).LT.A(JP,JP)) TAN=-TAN
      COS=1.d0/DSQRT(1.d0+TAN**2)
      SIN=TAN*COS
      AI=A(IP,IP)
      A(IP,IP)=(COS**2)*(AI+TAN*(2.d0*A(IP,JP)+TAN*A(JP,JP)))
      A(JP,JP)=(COS**2)*(A(JP,JP)-TAN*(2.d0*A(IP,JP)-TAN*AI))
      A(IP,JP)=0.d0
      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.d0) 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.d0
      MJ=I+1
      P(I)=0.d0
      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.d0
      P(JP)=0.d0
      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
      P(I)=A(I,I)
      ENDDO
      DO I=1,N
      DO J=1,N
      A(I,J)=W(I,J)
      W(I,J)=0.D0
      ENDDO
      W(I,I)=P(I)
      ENDDO
      RETURN
      END



C     ------------------------------------------------------------------
      SUBROUTINE CONCOR(M,X,F)
C     CONSTRUCTING VALID CORRELATION MATRICES
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /NEAREST/Z,MORDER
      COMMON /RNDM/IU,IV
      DIMENSION Z(10,10),A(10,10)
      DIMENSION X(*),V(10,10),W(10,10),P(10)
C     ==================================================================
C     CHECK THE NUMBER OF INVALID ELEMENTS
      MINVAL=0
      DO I=1,MORDER
      DO J=I,MORDER
      IF(DABS(Z(I,J)).GT.1.D0) MINVAL=MINVAL+1
      ENDDO
      ENDDO
      IF(M.NE.MINVAL) THEN
      WRITE(*,*)' '
      WRITE(*,*)'??????????????????????????????????????????????????????'
      WRITE(*,*)'PARAMETER DOES NOT MATCH.'
      WRITE(*,*)'THE VALUE OF M SHOULD BE=',MINVAL
      WRITE(*,*)'RERUN THE PROGRAM WITH M =', MINVAL
      STOP
      ENDIF
C     ==================================================================
        DO I=1,M
        IF(X(I).LT.-1.D0 .OR. X(I).GT.1.D0) THEN
        CALL RANDOM(RAND)
        X(I)=(RAND-.5d0)*2
        ENDIF
        ENDDO
C     CONSTRUCT THE MATRIX(MM,MM)
      ICOUNT=0
      DO I=1,MORDER
      A(I,I)=1.d0
         DO J=I+1,MORDER
            IF(DABS(Z(I,J)).GT.1.D0) THEN
            ICOUNT=ICOUNT+1
            A(I,J)=X(ICOUNT)
            ELSE
            A(I,J)=Z(I,J)
            ENDIF
         ENDDO
      ENDDO
C     FILLING THE LOWER DIAGONAL
      DO I=1,MORDER
      DO J=1,I
      A(I,J)=A(J,I)
      ENDDO
      ENDDO
C     FIND EIGENVALUES AND EIGENVECTORS OF MATRIX A
      CALL EIGEN(A,MORDER,V,W)
C     STORE EIGENVALUES (DIAGONAL OF RETURNING W) INTO P
      F=0.D0
      PSUM=0.D0 ! SUM OF MAGNITUDE OF EIGENVALUES
      DO I=1,MORDER
      P(I)=W(I,I)
      IF(P(I).LT.0.D0) PSUM=PSUM+DABS(P(I))
      ENDDO
      PROD=1.D0
      DO I=1,MORDER
      IF(P(I).LT.0.D0) THEN
      F=F+P(I)**2
      PROD=PROD*P(I)
      ENDIF
      ENDDO
      IF(PROD.LT.0.D0.OR.PROD.GT.1.D0) F=(F+PSUM+PROD**2)**2
      RETURN
      END
C     ------------------------------------------------------------------
      SUBROUTINE NCORX(X,M,OUTFIL)
C     NEAREST CORRELATION MATRIX PROBLEM
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /NEAREST/Z,MORDER
      DIMENSION Z(10,10),X(*)
      DIMENSION A(10,10),V(10,10),W(10,10)
      CHARACTER *70 OUTFIL
      CHARACTER *15 CFIL
C     -----------------------------------------------------------------
      CFIL='CORRESULTS'! IT IS AN INTERMEDIATE FILE
C     -----------------------------------------------------------------
      OPEN (7, FILE=CFIL) ! OPENS INTERMEDIATE FILE FOR INPUT
      OPEN(8, FILE=OUTFIL) ! OPENS OUTPUT FILE TO STORE VALID MATRICES
C     CONSTRUCT THE CORRELATION MATRIX
      READ(7,*) N ! READS POPULATION SIZE FROM INTERMEDIATE FILE
      DO IC=1,N
      READ(7,*)(X(J),J=1,M) ! READS VECTOR FROM INTERMEDIATE FILE

      ICOUNT=0
      DO I=1,MORDER
      A(I,I)=1.d0
         DO J=I+1,MORDER
            IF(DABS(Z(I,J)).GT.1.D0) THEN
            ICOUNT=ICOUNT+1
            A(I,J)=X(ICOUNT)
            ELSE
            A(I,J)=Z(I,J)
            ENDIF
         ENDDO
      ENDDO
C     FILLING THE LOWER DIAGONAL
      DO I=1,MORDER
      DO J=1,I
      A(I,J)=A(J,I)
      ENDDO
      ENDDO
      WRITE(*,*)' '
      WRITE(8,*)'******************************************************'
      WRITE(8,*)'A VALID CORRELATION MATRIX'
      DO I=1,MORDER
      WRITE(8,1)(A(I,J),J=1,MORDER)
      ENDDO
      WRITE(*,*) ' '
      CALL EIGEN(A,MORDER,V,W)
      MSIGN=0
      SUMW=0.D0
      PROD=1.D0
      DO I=1,MORDER
      SUMW=SUMW+W(I,I)
      PROD=PROD*W(I,I)
      IF(W(I,I).LT.0) MSIGN=1
      ENDDO
      WRITE(8,*)'EIGENVALUES, SUM AND PRODUCT OF EIGENVALUES'
      WRITE(8,1)(W(I,I),I=1,MORDER),SUMW,PROD
      WRITE(8,*)'EIGENVECTORS '
      DO I=1,MORDER
      WRITE(8,1)(V(I,J),J=1,MORDER)
      ENDDO
    1 FORMAT(8F10.7)
      IF(MSIGN.EQ.1) THEN
      WRITE(8,*)'FAILURE OF THE METHOD'
      ELSE
      WRITE(8,*)'SUCCESS OF THE METHOD'
      ENDIF
      ENDDO
      CLOSE(7)
      CLOSE(8)
      RETURN
      END