      SUBROUTINE GALRKN(HW,M,ZZ,EGVL,FNORM,
     1              BEGVL,UMATL,VMATL,MX,MULT,A,B,H,HL,
     2              INT,INTH,DEPW,GRADW,FINTW,NPW,AVGW,LHW,
     3              DEPB,RHOB,GRHOB,CGRADB,CINTB,NPB,AVGB,LHB)
C
C     THIS SUBROUTINE FINDS COMPLEX EIGENVALUES EGVL AND EIGEN-
C     FUNCTIONS USING THE GALERKIN-MATRIX METHOD. THE GALERKIN 
C     MATRIX METHOD REDUCES THE PROBLEM TO A GENERALIZED MATRIX 
C     EIGENVALUE PROBLEM OF THE FORM: A UMATL = B UMATL DIAG(EGVL).
C     THE EIGENVALUES ARE THE SAME FOR BOTH PROBLEMS AND THE EIGEN-
C     FUNCTIONS MAY BE CONSTRUCTED FROM THE MATRIX EIGENVECTORS IN
C     UMATL. THE GENERALIZED MATRIX EIGENVALUE PROBLEM IS SYMMETRIZED 
C     USING THE CHOLESKY DECOMPOSITION OF B TO ASSURE OTHOGONALITY 
C     OF THE EIGENFUNCTIONS CONSTRUCTED FROM THE MATRIX EIGENVECTORS
C     IN UMATL.
C 
C     THIS VERSION ALLOWS LINEAR VARIATONS IN THE SQUARE OF THE 
C     INDEX OF REFRACTION IN BOTH WATER AND BOTTOM. THE SQUARE OF
C     THE INDEX OF REFRACTON IN THE BOTTOM CAN BE COMPLEX DUE TO 
C     ATTENUATION. IT ALSO ALLOWS A PIECEWISE EXPONENTIAL DENSITY 
C     VARIATION IN THE BOTTOM.
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      COMMON /BLKEVN/ HB,CW,CB,FKW,FKB,ROHW,ROHB,ATEN
C
      DIMENSION INT(MX),INTH(MX)
C
      REAL*8 ZZ(MX),ZERO,ESP
C 
      DIMENSION DEPW(MX),GRADW(MX),FINTW(MX)
      DIMENSION DEPB(MX),RHOB(MX),GRHOB(MX)
      REAL*8 B(MX,MX),VMATL(MX,MX)
      COMPLEX*16 CGRADB(MX),CINTB(MX)
      COMPLEX*16 EGVL(MX),FNORM(MX),BEGVL(MX)
      COMPLEX*16 A(MX,MX),H(MX,MX),HL(MX,MX)
      COMPLEX*16 UMATL(MX,MX)
      COMPLEX*16 CDSRT,GAMA0
      COMPLEX*16 MULT(MX),CZERO,ONE
      COMPLEX*16 ALPHA1,BETA1
C
      ESP=.1D-9
      ZERO=DBLE(0.0)     
      ONE=DCMPLX(1.0,0.0)
      CZERO=DCMPLX(0.0,0.0)
C
C     FIND THE BASIC EIGENVALUES
C 
      CALL EGNVAL(HW,M,ZZ,MX)
C
C     FIND THE NORMALIZATION FACTORS FOR THE BASIC EIGENFUNCTIONS
C
      DO 100 I=1,M
C
      BEGVL(I)=DCMPLX(ZZ(I),ZERO)
      ALPHA1=CDSRT(FKW**2+FKB**2*(BEGVL(I)-ONE))
      BETA1=CDSRT(FKB**2*BEGVL(I))
      FNORM(I)=GAMA0(HW,ALPHA1,BETA1)
C
C      WRITE(*,50) I,BEGVL(I),FNORM(I),ALPHA1,BETA1
C      WRITE(9,50) I,BEGVL(I),FNORM(I),ALPHA1,BETA1
C   50 FORMAT(1X,I3,2X,4(2E12.4,1X))
C
  100 CONTINUE
C
C     COMPUTE THE REAL POSITIVE DEFINITE SYMMETRIC MATRIX B THAT
C     APPEARS ON THE RIGHT HAND SIDE OF THE GENERALIZED MATRIX EIGEN-
C     VALUE PROBLEM A UMATL = B UMATL DIAG(EGVL)
C
      CALL BMAT(HW,M,BEGVL,FNORM,B,MX,DEPB,RHOB,GRHOB,NPB,LHB)
C
C     COMPUTE THE CHOLESKY DECOMPOSTION OF B
C
C                             T
C          B  =   VMATL  VMATL
C
C     INTO THE LOWER TRIANGULAR MATRIX VMATL AND ITS TRANSPOSE
C
      CALL CHOLDC(B,VMATL,M,MX)
C
C     CONSTRUCT THE COMPLEX SYMMETRIC MATRIX A THAT APPEARS ON THE
C     LEFT HAND SIDE OF THE GENERALIZED MATRIX EIGENVALUE PROBLEM
C     A UMATL = B UMATL DIAG(EGVL) 
C
      CALL AMAT(HW,M,BEGVL,FNORM,A,MX,
     1          DEPW,GRADW,FINTW,NPW,AVGW,LHW,
     2          DEPB,RHOB,GRHOB,CGRADB,CINTB,NPB,AVGB,LHB)
C
C     PRINT THE MATRIX A (COMMENTED OUT)
C
C      CALL PRTMAT(M,M,'A   ',A,MX)
C
C     USE FORWARD SUBSTITUTION TWICE TO COMPUTE THE COMPLEX SYMMETRIC 
C     MATRIX
C                -1         T -1
C           VMATL   A  VMATL     ,    
C
C     WHICH IS STORED IN A, AND WHOSE EIGENVALUES ARE THE SAME AS THE
C     GENERALIZED EIGENVALUE PROBLEM A UMATL = B UMATL DIAG(EGVL) 
C
      CALL FORESUB(VMATL,HL,A,M,MX)
C
C     DEFINE THE TRANSPOSE OF HL
C
      DO 300 I=1,M
      DO 200 II=1,M
      H(II,I)=HL(I,II)
  200 CONTINUE
  300 CONTINUE
C
C     SOLVE FOR THE SYMETRIIC MATIX
C
      CALL FORESUB(VMATL,A,H,M,MX)
C
C     FIND THE EIGENVALUES EGVL AND EIGENVECTORS UMATL BY THE
C     QR ALGORITHM
C
      NCAL=M
      CALL CIGMAT(A,EGVL,M,MX,MULT,INT,INTH,H,UMATL,NCAL,HL)
C
      IF(NCAL .NE. M) THEN
      WRITE(*,*) 'NCAL,M= ',NCAL,M
      WRITE(9,*) 'NCAL,M= ',NCAL,M
      STOP 'SPECIFIED NUMBER OF EIGNEVALUES NOT FOUND'
      END IF
C
C     IN SOME CASES WHEN THE EIGENVALUES HAVE A SMALL IMAGINARY
C     PART, CIGMAT FINDS THE COMPLEX CONJUGATE OF THE PHYSICAL 
C     EIGENVALUE. IN THESE CASES THE COMPLEX CONJUGATE IS TAKEN.
C
      DO 600 I=1,M
      ICONJG=0
C
      IF((DIMAG(EGVL(I)) .GT. ZERO) .AND.
     1   (DIMAG(EGVL(I)) .LT. ESP )) THEN
      ICONJG=1
      EGVL(I)=DCONJG(EGVL(I))
      DO 500 II=1,M
      UMATL(II,I)=DCONJG(UMATL(II,I))
  500 CONTINUE
      END IF
C
      IF(ICONJG .NE. 0) THEN
      WRITE(9,*) '*** WARNING CONJUGATE EIGENVLAUE FOUND BY CIGMAT',
     1           ' MODE = ',I,' IN GALRKN ***'
      WRITE(*,*) '*** WARNING CONJUGATE EIGENVLAUE FOUND BY CIGMAT',
     1           ' MODE = ',I,' IN GALRKN ***'
      END IF
C
  600 CONTINUE
C
C     DEFINE THE TRANSPOSE OF VMATL
C
      DO 800 I=1,M
      DO 700 II=1,M
      H(II,I)=DCMPLX(VMATL(I,II),0.0D0)
  700 CONTINUE
  800 CONTINUE
C
C     RECOVER THE EIGENVECTORS OF THE GENERALIZED MATRIX EIGENVALUE
C     PROBLEM BY BACKSUBSTITUTION AND STORE IN HL.
C
      CALL BACKSUB(H,HL,UMATL,M,MX)
C
C     RETURN THE EIGENVECTORS IN UMATL AND STORE THE REAL SYMMETRIC
C     POSITIVE DEFINITE MATRIX B IN THE COMPLEX MATRIX A FOR ORTHO-
C     GONALITY CHECK
C
      DO 1000 I=1,M
      DO 900 II=1,M
      UMATL(II,I)=HL(II,I)
      A(II,I)=DCMPLX(B(II,I),0.0D0)
  900 CONTINUE
 1000 CONTINUE
C
C     CHECK ORTHOGONALITY (COMMENTED OUT)
C
C      CALL TMMULT(M,UMATL,A,UMATL,MX,H)
C      CALL PRTMAT(M,M,'ORTH',A,MX)
C
      RETURN 
      END      
