PROGRAM NCOR ! -------------------------------------------- C (OBTAINS POSITIVE SEMI-DEFINITE CORRELATION MATRIX (R) C FROM A NEGATIVE DEFINITE PSEUDO-CORRELATION MATRX (Q) IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (NORM=2) DIMENSION A(100,100),V(100,100),W(100,100),P(100),B(100,100) DIMENSION D(100,100),Q(100,100),C(100,100),MM(100) CHARACTER *80 INFIL,OUTFIL COMMON /RNDM/IU,IV C ----------------------------------------------------------------- WRITE(*,*)'THE ORDER OF Q MATRIX ?' READ(*,*) N WRITE(*,*)'INPUT AND OUTPUT FILE NAMES ?' READ(*,*) INFIL,OUTFIL WRITE(*,*)'NUMBER OF ITERATIONS(SAY 100) AND PERTURB(SAY -1.0D4)?' ! PERT IS A NEGATIVE LARGE NUMBER ! THAT IS USED IN PLACE OF NEGATIVE/ZERO EIGENVALUES IN Q ! THE MAGNITUDE MAY BE CHANGED BY THE USER READ(*,*) NIT,PERT C ----------------------------------------------------------------- OPEN(7,FILE=INFIL) ! OPENING THE INUT FILE DO I=1,N READ(7,*)(Q(I,J),J=1,N) ENDDO CLOSE(7) ! CLOSING THE INUT FILE C ----------------------------------------------------------------- OPEN(7,FILE=OUTFIL) !OPENING THE OUTPUT FILE, WILL BE CLOSED LATER WRITE(7,*)'THE ORIGINAL MATRIX (Q) IS GIVEN BELOW' DO I=1,N WRITE(7,1)(Q(I,J),J=1,N) ENDDO 1 FORMAT(10F8.4) C -----------------GENERATE MATRIX FOR EXPERIMENTS ---------------- C WRITE(*,*)'FEED SEED' C READ(*,*) IU C DO I=1,N C DO J=1,N C CALL RANDOM(RAND) C Q(I,J)=0.2D0+RAND*.7 ! OR BY SOME OTHER SPECIFICATION C IF(I.EQ.J) Q(I,J)=1.D0 C ENDDO C ENDDO C C ----------------------------------------------------------------- C ---------- STORE Q IN A: WORK WITH A, PRESERVE Q ---------------- DO I=1,N DO J=1,N A(I,J)=Q(I,J) ENDDO ENDDO CALL EIGEN(A,N,NN,V,W,P,MM) ! FIND EIGENVALUES AND VECTORS OF A WRITE(*,*)'THE EIGENVALUES OF THE ORIGINAL MATRIX ARE AS FOLLOWS' WRITE(7,*)'THE EIGENVALUES OF THE ORIGINAL MATRIX ARE AS FOLLOWS' WRITE(*,*)(A(I,I),I=1,N) WRITE(7,*)(A(I,I),I=1,N) WRITE(*,*)'-----------------------------------------------------' C A IS DISTURBED DUE TO INVOKING EIGEN. RESTORE IT. DO I=1,N DO J=1,N A(I,J)=Q(I,J) ENDDO ENDDO C ----------------------------------------------------------------- C ---------- FIND TRACE OF THE GIVEN MATRIX ------------------------ TRACE=0.D0 DO I=1,N TRACE=TRACE+A(I,I) ENDDO C NORMALIZE A MATRIX AND INITIALIZE D MATRIX AS AN IDENTITY MATRIX C AND INITIALIZE W MATRIX DO I=1,N DO J=1,N W(I,J)=0.D0 D(I,J)=0.D0 IF(I.EQ.J) D(I,J)=1.D0 A(I,J)=A(I,J)/TRACE ENDDO ENDDO C ----- VON NEUMANN DIVERGENCE ----- ITERATION BEGINS ------------- DO IT=1,NIT ! ITERATION BEGINS DO I=1,N DO J=1,N D(I,J)=0.D0 IF(I.EQ.J) D(I,J)=1.D0 ENDDO ENDDO DO I=1,N D(I,I)=DLOG(DSQRT((1.D0-A(I,I))/A(I,I))) ENDDO C FIND LOG(A) -------------------------------------------------- CALL EIGEN(A,N,NN,V,W,P,MM) ! FIND EIGENVALUES AND VECTORS OF A DO I=1,N IF(A(I,I).GT.0.D0) THEN W(I,I)=DLOG(A(I,I))! TO FIND LOG(A) ELSE W(I,I)=PERT ! PERT IS A LARGE (NEGATIVE) NUMBER ENDIF ENDDO DO I=1,N DO J=1,N C(I,J)=0.D0 DO K=1,N C(I,J)=C(I,J)+V(I,K)*W(K,J) ENDDO ENDDO ENDDO DO I=1,N DO J=1,N A(I,J)=0.D0 DO K=1,N A(I,J)=A(I,J)+C(I,K)*V(J,K) ENDDO ENDDO ENDDO C LOG(A) FOUND. ADD D TO A DO I=1,N DO J=1,N A(I,J)=A(I,J)+D(I,J) ENDDO ENDDO C FIND EXP(A+D) CALL EIGEN(A,N,NN,V,W,P,MM) DO I=1,N W(I,I)=DEXP(A(I,I)) ENDDO DO I=1,N DO J=1,N C(I,J)=0.D0 DO K=1,N C(I,J)=C(I,J)+V(I,K)*W(K,J) ENDDO ENDDO ENDDO DO I=1,N DO J=1,N A(I,J)=0.D0 DO K=1,N A(I,J)=A(I,J)+C(I,K)*V(J,K) ENDDO ENDDO ENDDO TRACE=0.D0 DO I=1,N TRACE=TRACE+A(I,I) ENDDO C NORMALIZE A MATRIX AND INITIALIZE D MATRIX AS AN IDENTITY MATRIX DO I=1,N DO J=1,N A(I,J)=A(I,J)/TRACE ENDDO ENDDO DO I=1,N DO J=1,N B(I,J)=A(I,J) ENDDO ENDDO CALL EIGEN(B,N,NN,V,W,P,MM) ENDDO ! ITERATION ENDS NCH=0 ! FOR CHECK IF ALL EIGENVALUES ARE NON-NEGATIVE DO I=1,N IF(B(I,I).LE.0.D0) NCH=1 ENDDO C ----------------------------------------------------------------- C RESTUCTURING THE MATRIX TRACE=0.D0 DO I=1,N TRACE=TRACE+B(I,I) ENDDO DO I=1,N W(I,I)=B(I,I)*N/TRACE ENDDO DO I=1,N DO J=1,N C(I,J)=0.D0 DO K=1,N C(I,J)=C(I,J)+V(I,K)*W(K,J) ENDDO ENDDO ENDDO DO I=1,N DO J=1,N A(I,J)=0.D0 DO K=1,N A(I,J)=A(I,J)+C(I,K)*V(J,K) ENDDO ENDDO ENDDO C COMPUTE NORM S=0.D0 DO I=1,N DO J=1,N S=S+DABS(A(I,J)-Q(I,J))**2 ENDDO ENDDO C ----------------------------------------------------------------- C FINAL RESULTS WRITE(7,*)'FINAL RESULTS - THE OUTPUT MATRIX' DO I=1,N WRITE(7,1)(A(I,J),J=1,N) ENDDO WRITE(7,*)'------------THE FINAL EIGENVALUES ARE --------------' WRITE(7,*)(B(I,I)*N/TRACE,I=1,N) WRITE(7,*)'NORM=',S**(1.D0/NORM) CLOSE(7)! CLOSE THE OUTPUT FILE C ----------------------------------------------------------------- WRITE(*,*)'RESULTS HAVE BEEN STORED IN THE OUTPUT FILE' IF(NCH.NE.0) THEN WRITE(*,*)'WARNING ---- WARNING ------ WARNING -------' WRITE(*,*)'SOME EIGENVALUES ARE NEGATIVE EVEN NOW --------------' WRITE(*,*)'RE-RUN WITH MORE ITERATIONS OR OTHER PERTURBANCE VALUE' WRITE(*,*)'SOME TRIAL AND ERROR MAY BE NEEDED' ENDIF WRITE(*,*)'END' END C ----------------------------------------------------------------- SUBROUTINE EIGEN(A,N,NN,V,W,P,MM) C THIS SUBROUTINE IS ADAPTED FROM KRISHNAMURTHY AND SEN (1976) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION A(100,100),V(100,100),W(100,100),P(100) DIMENSION MM(100) C ------------ INITIALISATION ------------------------- NN=1 DO 50 I=1,N DO 51 J=1,N V(I,J)=0.0 51 W(I,J)=0.0 P(I)=0.0 50 CONTINUE PMAX=0 EPLN=0 TAN=0 SIN=0 COS=0 AI=0 TT=0 EPLN=1.0D-310 C 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 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=DBLE(IV+.0D0) IU=IV RAND=RAND*0.4656613E-09 RETURN END C -----------------------------------------------------------------