C MATXTEST - CDC - L.P.GERALDO - AP 314 C TESTS A SYMMETRIC MATRIX FOR POSITIVE DEFINITENESS. C THIS PROGRAM ALSO COMPUTATES THE DETERMINANT AND C INVERSE MATRIX C C C OPTIONS FOR MATRIX INPUT (IOP) C C OPTION 1 - INPUT ERRORS AND CORRELATION MATRIX C OPTION 2 - INPUT COVARIANCE MATRIX C OPTION 3 - INPUT CROSS SECTIONS AND RELATIVE COVARIANCE MATRIX C C OPTIONS FOR MATRIX INFORMATION OUTPUT (IND) C TOGETHER WITH POSITIVE DEFINITENESS TEST. C C OPTION - 00 = STOP C OPTION - 10 = ONLY DETERMINANT C OPTION - 01 = ONLY INVERSE MATRIX C OPTION - 11 = BOTH DETERMINANT AND INVERSE MATRIX C C REAL A(50,50),C(50,50),F(50,50),E(50),X(50),DET(3),T,D REAL V(50,50),W(3,3),TS,XS,YS,ZS DATA LDA/50/ DATA LD/3/ DATA IN/'N '/ WRITE(*,5) 5 FORMAT(1X,'DATA INPUT: ') OPEN(3,FILE=' ',STATUS='OLD ') READ(3,10)N 10 FORMAT(I2) READ(3,10)IOP READ(3,10)IND 15 FORMAT(6E12.6) GO TO (20,40,50) IOP 20 CONTINUE READ(3,15) (E(I), I=1,N) DO 25 I=1,N 25 READ(3,15) (C(I,J), J=1,I) DO 30 I=1,N DO 30 J=1,I 30 C(J,I)=C(I,J) DO 35 I=1,N DO 35 J=1,N 35 A(I,J) = E(I)*E(J)*C(I,J) DO 36 I=1,N DO 36 J=1,N 36 V(I,J)=A(I,J) GO TO 70 40 CONTINUE DO 45 I=1,N 45 READ(3,15) (A(I,J), J=1,I) DO 46 I=1,N DO 46 J=1,N 46 A(I,J)=A(J,I) DO 47 I=1,N DO 47 J=1,N 47 C(I,J)=A(I,J)/(A(I,I)*A(J,J))**0.5 DO 48 I=1,N DO 48 J=1,N 48 V(I,J)=A(I,J) GO TO 70 50 CONTINUE READ(3,15) (X(I), I=1,N) DO 55 I=1,N 55 READ(3,15) (F(I,J), J=1,I) DO 60 I=1,N DO 60 J=1,N 60 F(I,J)=F(J,I) DO 65 I=1,N DO 65 J=1,N 65 A(I,J)=X(I)*X(J)*F(I,J) DO 67 I=1,N DO 67 J=1,N 67 C(I,J)=A(I,J)/(A(I,I)*A(J,J))**0.5 DO 68 I=1,N DO 68 J=1,N 68 V(I,J)=A(I,J) 70 CONTINUE C C CORRELATION MATRIX TEST C DO 73 I=1,N DO 73 J=1,I IF(I .EQ. J) GO TO 73 IF(ABS(C(I,J)) .LE. 1) GO TO 73 WRITE(*,71) 71 FORMAT(1X,'... I ....... J ....... C(I,J) ... ') WRITE(*,72)I,J,C(I,J) 72 FORMAT(1X,I5,I10,F18.6) 73 CONTINUE DO 74 I=1,N DO 74 J=1,I IF(I .EQ. J) GO TO 74 IF(ABS(C(I,J)) .GT. 1) GO TO 240 74 CONTINUE C C POSITIVE DEFINITENESS TEST C CALL SPOFA(A,LDA,N,INFO) IF (INFO .NE. 0 ) GO TO 125 C C C OUTPUT DATA C WRITE(*,80) 80 FORMAT(1X,' THIS MATRIX IS POSITIVE DEFINITE ') IF (IND .EQ. 0) GO TO 120 C C DETERMINANT AND INVERSE MATRIX COMPUTATION C CALL SPODI(A,LDA,N,DET,IND) C IF (IND .EQ. 11) GO TO 85 IF (IND .EQ. 1) GO TO 95 85 DT=DET(1)*10.0**DET(2) WRITE(*,90) DT 90 FORMAT(1X, 'DET= ',E12.6) IF (IND .EQ. 10) GO TO 120 91 WRITE(*,92) 92 FORMAT(1X,'Do you want to save the inverse matrix?Y/N ') READ(*,143,ERR=91)IOP IF(IOP .EQ. IN) GO TO 95 WRITE(*,93) 93 FORMAT(1X,'Please type the output file name ') OPEN(5,FILE=' ',STATUS='NEW ') WRITE(5,94) 94 FORMAT(1X,'xxxxxxxxxx INVERSE MATRIX xxxxxxxxxx ') DO 941 I=1,N DO 941 J=1,N 941 A(J,I)=A(I,J) DO 942 J=1,N 942 WRITE(5,115) (A(J,I), I=1,J) GO TO 120 C 95 WRITE(*,100) 100 FORMAT(1X,'xxxxxxxxxx INVERSE MATRIX xxxxxxxxxx') DO 105 I=1,N DO 105 J=1,N 105 A(J,I)=A(I,J) DO 110 J=1,N 110 WRITE(*,115) (A(J,I), I=1,J) 115 FORMAT(1X,50E12.6) 120 STOP 125 WRITE(*,130) 130 FORMAT(1X,'THIS MATRIX IS NOT POSITIVE DEFINITE ') WRITE(*,135) 135 FORMAT(1X,'VERIFY THE FOLLOWING MATRIX LINE: ') I=INFO WRITE(*,140)I 140 FORMAT(1X,'LINE I=',I2) WRITE(*,115) (A(I,J),J=1,I) 141 WRITE(*,142) 142 FORMAT(1X,'Do you want to test all the cross correlations inside t 1his leading principal minor?Y/N ') READ(*,143,ERR=141) IOP 143 FORMAT(A1) IF(IOP .EQ. IN) GO TO 235 WRITE(*,145) 145 FORMAT(1X,'Test of the cross correlation consistency ') N=INFO IF(N .LT. 3) GO TO 238 150 I1=1 160 J1=I1+1 170 K1=J1+1 178 W(1,1)=V(I1,I1) W(2,2)=V(J1,J1) W(3,3)=V(K1,K1) W(1,2)=V(I1,J1) W(1,3)=V(I1,K1) W(2,3)=V(J1,K1) M=3 CALL SPOFA(W,LD,M,INF) IF(INF .EQ. 0) GO TO 188 WRITE(*,185) 185 FORMAT(1X,'This cross correlation is not positive definite. ') WRITE(*,187)I1,J1,K1 187 FORMAT(1X,'The correspondent parameters are: ',3I4) 188 IF(K1 .EQ. N) GO TO 200 K1=K1+1 GO TO 178 200 IF(J1 .EQ. N-1) GO TO 210 J1=J1+1 GO TO 170 210 IF(I1 .EQ. N-2) GO TO 220 I1=I1+1 GO TO 160 220 CONTINUE WRITE(*,230) 230 FORMAT(1X,'All 3x3 principal minors were tested ') 235 STOP 238 WRITE(*,239) 239 FORMAT(1X,'The matrix order is lower than 3. ') STOP 240 WRITE(*,245) 245 FORMAT(1X,'This matrix is unreal because the correlation coefficie 1nts above are higher than 1 ') STOP END C C C SUBROUTINE SPOFA(A,LDA,N,INFO) C SPOFA TESTS IF THE MATRIX IS POSITIVE DEFINITE AND C FACTORS A REAL SYMMETRIC POSITIVE DEFINITE MATRIX C C INTEGER LDA,N,INFO REAL A(LDA,1) REAL SDOT,T REAL S INTEGER J,JM1,K C BEGIN BLOCK WITH...EXITS TO 40 C C DO 30 J = 1, N INFO = J S = 0.0E0 JM1 = J - 1 IF (JM1 .LT. 1) GO TO 20 DO 10 K = 1, JM1 T = A(K,J) - SDOT(K-1,A(1,K),1,A(1,J),1) T = T/A(K,K) A(K,J) = T S = S + T*T 10 CONTINUE 20 CONTINUE S = A(J,J) - S C EXIT IF (S .LE. 0.0E0) GO TO 40 A(J,J) = SQRT(S) 30 CONTINUE INFO = 0 40 CONTINUE RETURN END C C C SUBROUTINE SPODI(A,LDA,N,DET,JOB) C SPODI COMPUTES THE DETERMINANT AND INVERSE OF A CERTAIN C REAL SYMMETRIC POSITIVE DEFINITE MATRIX. C C INTEGER LDA,N,JOB REAL A(LDA,1) REAL DET(2) C INTERNAL VARIABLES REAL T REAL S INTEGER I,J,JM1,K,KP1 C C COMPUTE DETERMINANT C IF (JOB/10 .EQ. 0) GO TO 70 DET(1) = 1.0E0 DET(2) = 0.0E0 S = 10.0E0 DO 50 I = 1,N DET(1) = A(I,I)**2*DET(1) C ...EXIT IF (DET(1) .EQ. 0.0E0) GO TO 60 10 IF (DET(1) .GE. 1.0E0) GO TO 20 DET(1) = S*DET(1) DET(2) = DET(2) - 1.0E0 GO TO 10 20 CONTINUE 30 IF (DET(1) .LT. S) GO TO 40 DET(1) = DET(1)/S DET(2) = DET(2) + 1.0E0 GO TO 30 40 CONTINUE 50 CONTINUE 60 CONTINUE 70 CONTINUE C C COMPUTE INVERSE(R) C IF (MOD(JOB,10) .EQ. 0) GO TO 140 DO 100 K = 1, N A(K,K) = 1.0E0/A(K,K) T = -A(K,K) CALL SSCAL(K-1,T,A(1,K),1) KP1 = K + 1 IF (N .LT. KP1) GO TO 90 DO 80 J = KP1, N T = A(K,J) A(K,J) = 0.0E0 CALL SAXPY(K,T,A(1,K),1,A(1,J),1) 80 CONTINUE 90 CONTINUE 100 CONTINUE C C FORM INVERSE(R) * TRANS(INVERSE(R)) C DO 130 J = 1, N JM1 = J - 1 IF (JM1 .LT. 1) GO TO 120 DO 110 K = 1, JM1 T = A(K,J) CALL SAXPY(K,T,A(1,J),1,A(1,K),1) 110 CONTINUE 120 CONTINUE T = A(J,J) CALL SSCAL(J,T,A(1,J),1) 130 CONTINUE 140 CONTINUE RETURN END C C SUBROUTINE SSCAL(N,SA,SX,INCX) C C SSCAL SCALES A VECTOR BY A CONSTANT C REAL SA,SX(1) INTEGER I,INCX,M,MP1,N,NINCX C M=MOD(N,5) IF(M .EQ. 0) GO TO 40 DO 30 I=1,M SX(I)=SA*SX(I) 30 CONTINUE IF(N .LT. 5) RETURN 40 MP1=M+1 DO 50 I=MP1,N,5 SX(I)=SA*SX(I) SX(I+1)=SA*SX(I+1) SX(I+2)=SA*SX(I+2) SX(I+3)=SA*SX(I+3) SX(I+4)=SA*SX(I+4) 50 CONTINUE RETURN END C C C SUBROUTINE SAXPY(N,SA,SX,INCX,SY,INCY) C C SAXPY COMPUTES THE OPERATION: CONSTANT TIMES A VECTOR C PLUS A VECTOR. C REAL SX(1),SY(1),SA INTEGER I,INCX,INCY,IX,IY,M,MP1,N C M=MOD(N,4) IF(M .EQ. 0) GO TO 40 DO 30 I=1,M SY(I)=SY(I)+SA*SX(I) 30 CONTINUE IF(N .LT. 4) RETURN 40 MP1=M+1 DO 50 I=MP1,N,4 SY(I)=SY(I)+SA*SX(I) SY(I+1)=SY(I+1)+SA*SX(I+1) SY(I+2)=SY(I+2)+SA*SX(I+2) SY(I+3)=SY(I+3)+SA*SX(I+3) 50 CONTINUE RETURN END C C C REAL FUNCTION SDOT(N,SX,INCX,SY,INCY) C C SDOT FORMS THE DOT PRODUCT OF TWO VECTORS C REAL SX(1),SY(1),STEMP INTEGER I,INCX,INCY,IX,IY,M,MP1,N C STEMP=0.0E0 SDOT=0.0E0 IF(N .LE. 0)RETURN 20 M=MOD(N,5) IF(M .EQ. 0) GO TO 40 DO 30 I=1,M STEMP=STEMP+SX(I)*SY(I) 30 CONTINUE IF(N .LT. 5) GO TO 60 40 MP1=M+1 DO 50 I=MP1,N,5 STEMP=STEMP+SX(I)*SY(I)+SX(I+1)*SY(I+1)+SX(I+2)*SY(I+2)+ * SX(I+3)*SY(I+3)+SX(I+4)*SY(I+4) 50 CONTINUE 60 SDOT=STEMP RETURN END