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