C*********************************************************************** C PROGRAM "RANK OF G AND F MATRICES" BY ROBERT FRACZKIEWICZ * C DEPARTMENT OF CHEMISTRY, UNIVERSITY OF HOUSTON, 1992 * C CHEM86@JETSON.UH.EDU * C*********************************************************************** IMPLICIT REAL*8(A-H,O-Z) DIMENSION NR(2000),NC(2000),NFO(2000),Z(2000),AL(100,100), 1 AK(100,100),NROW(4),NCOL(4),NPO(4),DATIN(4),FI(150) CHARACTER*60 NAME PRINT * PRINT *,' RANK OF REAL SYMMETRIC MATRIX ' PRINT *,' RANKFG VERSION 4.0 (1992) ' PRINT * PRINT *,'NF=NUMBER OF INDEPENDENT FORCE CONSTANTS. IF NF=0 THEN' PRINT *,'PROGRAM CALCULATES RANK OF G-MATRIX, OTHERWISE THAT OF ' PRINT *,'F-MATRIX. FORCE CONSTANTS FILE IS NEEDED IN THE LATTER ' PRINT *,'CASE (EXTENSION ".FI")' PRINT *,'INPUT THE TOTAL NUMBER OF FORCE CONSTANTS:' READ(*,*) NF PRINT *,'INPUT THE name OF THE FILE (WITHOUT EXTENSION !!!) ' PRINT *,'WHERE THE Z- OR G-MATRIX IS STORED AND PRESS . ' PRINT *,'THE ".ZMAT" OR ".GMAT" EXTENSION WILL BE AUTOMATICALLY ' PRINT *,'ASSUMED. THE LENGTH OF THE name CANNOT EXCEED 14 ' PRINT *,'CHARACTERS. ' READ(*,'(A60)') NAME NAMLEN = LSTNBL(NAME) PRINT *,'INPUT THE DIMENSION OF THE MATRIX' READ(*,*) NQ PRINT *,'THE RESULTS OF THE PROGRAM WILL BE STORED IN THE' PRINT *,'OUTPUT FILE WITH EXTENSION ".RAOUTF" OR ".RAOUTG".' NRANK=NQ DO 180 I=1,2000 180 NFO(I)=0 IF(NF.GT.0) THEN NOZ=0 OPEN(UNIT=16, FILE=NAME(1:NAMLEN)//'.RAOUTF', STATUS='NEW') OPEN(UNIT=15, FILE=NAME(1:NAMLEN)//'.ZMAT', 1 ERR=9991, STATUS='OLD') 191 READ(15,18) (NROW(L),NCOL(L),NPO(L),DATIN(L),L=1,4) 18 FORMAT(4(3I3,F9.6)) DO 196 L=1,4 IF(NROW(L))198,196,193 193 NOZ=NOZ+1 NR(NOZ)=NROW(L) NC(NOZ)=NCOL(L) NFO(NOZ)=NPO(L) Z(NOZ)=DATIN(L) IF(NROW(L).NE.NCOL(L)) THEN NOZ=NOZ+1 NR(NOZ)=NCOL(L) NC(NOZ)=NROW(L) NFO(NOZ)=NPO(L) Z(NOZ)=DATIN(L) ENDIF 196 CONTINUE GO TO 191 198 CLOSE(UNIT=15) PRINT *,'INPUT THE name OF THE FILE (WITHOUT EXTENSION !!!) ' PRINT *,'WHERE THE FORCE CONSTATNTS ARE STORED AND PRESS ' PRINT *,'. THE ".FI" EXTENSION WILL BE AUTOMATICALLY ' PRINT *,'ASSUMED. THE LENGTH OF THE name CANNOT EXCEED 14 ' PRINT *,'CHARACTERS. ' READ(*,'(A60)') NAME NAMLEN = LSTNBL(NAME) OPEN(UNIT=15, FILE=NAME(1:NAMLEN)//'.FI', 1 ERR=9991, STATUS='OLD') READ(15,*) (FI(I), I=1,NF) CLOSE(UNIT=15) DO 195 I=1,NOZ 195 Z(I)=Z(I)*FI(NFO(I)) ELSE NOZ=0 OPEN(UNIT=16, FILE=NAME(1:NAMLEN)//'.RAOUTG', STATUS='NEW') OPEN(UNIT=15, FILE=NAME(1:NAMLEN)//'.GMAT', 1 ERR=9991, STATUS='OLD') 199 READ(15,19) (NROW(L),NCOL(L),DATIN(L),L=1,4) 19 FORMAT(4(2I3,F12.6)) DO 201 L=1,4 IF(NROW(L))202,201,200 200 NOZ=NOZ+1 NR(NOZ)=NROW(L) NC(NOZ)=NCOL(L) Z(NOZ)=DATIN(L) IF(NROW(L).NE.NCOL(L)) THEN NOZ=NOZ+1 NR(NOZ)=NCOL(L) NC(NOZ)=NROW(L) Z(NOZ)=DATIN(L) ENDIF 201 CONTINUE GO TO 199 202 CLOSE(UNIT=15) ENDIF DO 100 I=1,100 DO 100 J=1,100 AL(I,J)=0.0 100 AK(I,J)=0.0 DO 110 L=1,NOZ I=NR(L) J=NC(L) 110 AK(I,J)=AK(I,J)+Z(L) NW=1 IF(ABS(AK(NW,NW)).LT.1.0E-05) THEN 20 AX=0.0 DO 10 I=NW+1,NQ AA=ABS(AK(NW,I)) IF(AA.GT.AX) THEN AX=AA KOL=I ENDIF 10 CONTINUE IF(AX.LT.9.0E-06) THEN NRANK=NRANK-1 WRITE(16,41) NW,AX 41 FORMAT(1X,I3,' - TH ROW IS ALMOST EQUAL TO ZERO, MAX = ',E13.6) Z(NW)=0.0 NW=NW+1 GO TO 20 ENDIF DO 30 L=NW,NQ POM=AK(L,NW) AK(L,NW)=AK(L,KOL) 30 AK(L,KOL)=POM ENDIF Z(NW)=AK(NW,NW) DO 40 I=NW+1,NQ DO 50 J=NW,I-1 IF(ABS(Z(J)).LT.1.0E-06) THEN AL(I,J)=0.0 ELSE AL(I,J)=AK(I,J) DO 60 K=1,J-1 60 AL(I,J)=AL(I,J)-AL(I,K)*AL(J,K)*Z(K) AL(I,J)=AL(I,J)/Z(J) ENDIF 50 CONTINUE Z(I)=AK(I,I) DO 70 K=1,I-1 70 Z(I)=Z(I)-AL(I,K)*AL(I,K)*Z(K) IF(ABS(Z(I)).LT.1.0E-05) NRANK=NRANK-1 40 CONTINUE WRITE(16,45) NAME(1:NAMLEN),NRANK 45 FORMAT(1X,'FILE : ',A14,/1X,'RANK OF THE MATRIX = ',i3) WRITE(16,43) 43 FORMAT(/1x,'THE MATRIX AFTER SINGULAR VALUE DECOMPOSITION'//) DO 80 K=1,NQ 80 WRITE(16,44) K,K,Z(K) 44 FORMAT(1X,'A(',I2,',',I2,')=',E13.6) CLOSE(UNIT=16) STOP 1 C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 9991 PRINT * PRINT *,'INPUT OPERATION UNSUCCESFUL. CANNOT FIND INPUT FILE' STOP 3 C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ END C ------------------- LSTNBL ----- C Routine which finds last nonblank character in a character variable C INTEGER FUNCTION LSTNBL(CHRVAR) CHARACTER*(*) CHRVAR INTEGER I, L L = LEN(CHRVAR) DO 100 I = L, 1, -1 IF(CHRVAR(I:I) .NE. ' ')GOTO 200 100 CONTINUE 200 LSTNBL = I RETURN END