PROGRAM NCORMAT ! -------------------------------------------- C (OBTAINS POSITIVE SEMI-DEFINITE CORRELATION MATRIX (R) C FROM A NON-POSITIVE-SEMI-DEFINITE or PSEUDO-CORRELATION MATRX (Q) C Non-Positive-Semi-definite matrix has at least one - BUT NOT ALL - C of its eigenvalues negative IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (NORM=2) DIMENSION A(200,200),V(200,200),W(200,200),P(200),R(200,200) DIMENSION D(200,200),Q(200,200),C(200,200),MM(200),QQ(200,200) CHARACTER *80 INFIL,OUTFIL COMMON /RNDM/IU,IV C ----------------------------------------------------------------- WRITE(*,*)'THE ORDER OF Q MATRIX ?' READ(*,*) N WRITE(*,*)'MATRIX Q TO BE READ OR GENERATED? TO BE READ FEED(0)' WRITE(*,*)'IF Q TO BE GENERATED, THEN FEED(ANY NON-ZERO NUMBER)' READ(*,*) IGEN WRITE(*,*)'OUTPUT FILE NAME (TO STORE THE RESULTS)?' READ(*,*) OUTFIL WRITE(*,*)'NUMBER OF ITERATIONS(SAY 100) AND ACC (SAY 1.0D-10)?' WRITE(*,*)'ACC IS A SMALL NUMBER FOR ACCURACY' READ(*,*) NIT,ACC WRITE(*,*)'FEED A PERTURBATION VALUE TO REPLACE NONPOSITIVE EIGEN' WRITE(*,*)'VALUE. IT WOULD BE A SMALL POSITIVE VALUE, SAY, 1.D-10' READ(*,*) EPS C ----------------------------------------------------------------- IF(IGEN.EQ.0) THEN WRITE(*,*)'INPUT FILE NAME (TO READ THE Q MATRIX)?' READ(*,*) INFIL 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 -----------------GENERATE MATRIX FOR EXPERIMENTS ---------------- ELSE WRITE(*,*)'FEED SEED' READ(*,*) IU DO I=1,N DO J=1,N CALL RANDOM(RAND) Q(I,J)=0.2D0+RAND*.7 ! OR BY SOME OTHER SPECIFICATION IF(I.EQ.J) Q(I,J)=1.D0 ENDDO ENDDO ENDIF C ----------------------------------------------------------------- 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 ---------- 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,INITIALIZE OTHERS DO I=1,N DO J=1,N W(I,J)=0.D0 !INITIALISE AND KEEP NONNEGATIVE EIGENVALUES IN W(I,I) IF(A(I,J).GT.0.D0 .AND. I.EQ.J) W(I,J)=A(I,J) IF(I.EQ.J.AND.W(I,J).EQ.0.D0) W(I,J)=EPS A(I,J)=Q(I,J) ! RESTORE THE ORIGINAL MATRIX IN A D(I,J)=0.D0 ! INTIALIZE D MATRIX BY ZERO ENDDO ENDDO C ----------------------------------------------------------------- S=1.0D30 ! NORM VERY LARGE TO BEGIN WITH DO IT=1,NIT ! ITERATION BEGINS DO I=1,N DO J=1,N R(I,J)=A(I,J)-D(I,J) A(I,J)=R(I,J) ENDDO ENDDO CALL EIGEN(A,N,NN,V,W,P,MM) ! FIND EIGENVALUES AND VECTORS OF A DO I=1,N DO J=1,N W(I,J)=0.D0 IF(I.EQ.J .AND. A(I,J).GT.0.D0) W(I,J)=A(I,J) IF(I.EQ.J.AND.W(I,J).EQ.0.D0) W(I,J)=EPS ENDDO ENDDO C CONSTRUCTING THE MATRIX FROM EIGENVECTORS AND NON-NEG EIGENVALUES 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 MATRIX HAS BEEN RECONSTRUCTED DO I=1,N DO J=1,N D(I,J)=A(I,J)-R(I,J) IF(I.EQ.J) A(I,J)=1.D0 ENDDO ENDDO C COMPUTE NORM SS=0.D0 DO I=1,N DO J=1,N SS=SS+DABS(A(I,J)-Q(I,J))**2 ENDDO ENDDO IF(DABS(SS-S).LT.ACC) THEN S=SS GO TO 10 ELSE S=SS ENDIF ENDDO ! ITERATIONS END C ----------------------------------------------------------------- C FINAL RESULTS 10 WRITE(*,*)'NO. OF ITERATIONS PERFORMED = ',IT WRITE(7,*)'FINAL RESULTS - THE OUTPUT MATRIX' DO I=1,N WRITE(7,1)(A(I,J),J=1,N) DO J=1,N R(I,J)=A(I,J) ENDDO ENDDO WRITE(7,*)'------THE FINAL EIGENVALUES ARE --------' CALL EIGEN(R,N,NN,V,W,P,MM) WRITE(7,*)(R(I,I),I=1,N) WRITE(7,*)'NORM=',SS**(1.D0/NORM) WRITE(7,*)'NO. OF ITERATIONS PERFORMED = ',IT CLOSE(7)! CLOSE THE OUTPUT FILE C ----------------------------------------------------------------- WRITE(*,*)'RESULTS HAVE BEEN STORED IN THE OUTPUT FILE' 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(200,200),V(200,200),W(200,200),P(200) DIMENSION MM(200) 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 -----------------------------------------------------------------