C     ..................................................................
C                                                                       
C        SUBROUTINE CIGMAT                                              
C                                                                       
C        PURPOSE                                                        
C           COMPUTE EIGENVALUES AND EIGENVECTORS OF A COMPLEX MATRIX    
C           COMPLEX*16 VERSION                                          
C                                                                       
C        USAGE                                                          
C           CALL CIGMAT(A,LAMBDA,M,MAX,MULT,INT,INTH,H,VECT,NCAL,HL)    
C                                                                       
C        DESCRIPTION OF PARAMETERS                                      
C           A(MAX,MAX)-      COMPLEX INPUT MATRIX                       
C           LAMBDA(MAX)-     COMPLEX EIGENVALUES IN ORDER BY INCREASING 
C                            REAL PART                                  
C           M-               ORDER OF A                                 
C           MAX-             MAXIMUM ORDER OF A FROM CALLING PROGRAM    
C           MULT(MAX)-       TEMPORARY COMPLEX ARRAY                    
C           INT(MAX)-        TEMPORARY INTEGER ARRAY                    
C           INTH(MAX)-       TEMPORARY INTEGER ARRAY                    
C           H(MAX,MAX)-      TEMPORARY COMPLEX MATRIX                   
C           VECT(MAX,MAX)-   COMPLEX EIGENVECTORS IN SAME SEQUENCE AS   
C                            EIGENVALUES                                
C           NCAL-            NUMBER OF EIGENVALUES(AND EIGENVECTORS)    
C                            DESIRED                                    
C           HL(MAX,MAX)-     TEMPORARY COMPLEX MATRIX                   
C                                                                       
C        REMARKS                                                        
C           THE 11 SUBROUTINE ARGUMENTS INCLUBE 5 COMPLEX 2-DIM         
C           ARRAYS, 1 COMPLEX 1-DIM ARRAY, 2 INTEGER 1-DIM              
C           ARRAYS, AND 3 INTEGER SCALARS.  BE CERTAIN THAT ARGUMENT    
C           TYPES IN MAIN PROGRAM ARE MATCHED WITH THESE IN SUBROUTINE. 
C                                                                       
C        METHOD                                                         
C           QR ALGORITHM OBTAINED FROM M. WERBY                         
C                                                                       
C     ..................................................................
C                                                                       
      SUBROUTINE CIGMAT (A,LAMBDA,M,MAX,MULT,INT,INTH,H,VECT,NCAL,HL)   
      IMPLICIT REAL*8 (A-H,O-Z)
      COMPLEX*16  A(MAX,MAX),LAMBDA(MAX),SHIFT(3),TEMP,TEMP1,TEMP2,SINN,
     1H(MAX,MAX),HL(MAX,MAX),VECT(MAX,MAX),MULT(MAX),COSS               
      LOGICAL TWICE                                                     
      REAL*8 EPS,SUM,ZERO,BIG,RHO,ABIG,ABSSQ,FN,DREAL                   
      INTEGER INT(MAX),INTH(MAX),R,RP1,RP2                              
      ZERO=DBLE(0.0)                                                    
      DO 10 I=1,M                                                       
      DO 10 J=1,M                                                       
10    HL(J,I)=A(J,I)                                                    
C
      NV=NCAL                                                           
      N=M                                                               
      NCAL=N                                                            
      IF (N.NE.1) GO TO 20                                              
      LAMBDA(1)=A(1,1)                                                  
      VECT(1,1)=DCMPLX(1.0,0.0)                                         
      GO TO 640                                                         
20    ICOUNT=0                                                          
      SHIFT(1)=DCMPLX(0.0,0.0)                                          
      IF (N.NE.2) GO TO 50                                              
30    TEMP=(A(1,1)+A(2,2)+CDSQRT((A(1,1)+A(2,2))**2-4.0*(A(2,2)*A(1,1)- 
     1A(2,1)*A(1,2))))/2.0                                              
      IF (DREAL(TEMP).NE.ZERO.OR.DIMAG(TEMP).NE.ZERO) GO TO 40          
      LAMBDA(M)=SHIFT(1)                                                
      LAMBDA(M-1)=A(1,1)+A(2,2)+SHIFT(1)                                
      GO TO 390                                                         
40    LAMBDA(M)=TEMP+SHIFT(1)                                           
      LAMBDA(M-1)=(A(2,2)*A(1,1)-A(2,1)*A(1,2))/(LAMBDA(M)-SHIFT(1))+SHI
     1FT(1)                                                             
      GO TO 390                                                         
C                                                                       
C        REDUCE MATRIX TO HESSENBERG FORM                               
C                                                                       
50    NM2=N-2                                                           
      DO 160 R=1,NM2                                                    
      RP1=R+1                                                           
      RP2=R+2                                                           
      ABIG=ZERO                                                         
      INT(R)=RP1                                                        
      DO 60 I=RP1,N                                                     
      ABSSQ=DREAL(A(I,R))**2+DIMAG(A(I,R))**2                           
      IF (ABSSQ.LE.ABIG) GO TO 60                                       
      INT(R)=I                                                          
      ABIG=ABSSQ                                                        
60    CONTINUE                                                          
      INTER=INT(R)                                                      
      IF (ABIG.EQ.ZERO) GO TO 160                                       
      IF (INTER.EQ.RP1) GO TO 90                                        
      DO 70 I=R,N                                                       
      TEMP=A(RP1,I)                                                     
      A(RP1,I)=A(INTER,I)                                               
70    A(INTER,I)=TEMP                                                   
      DO 80 I=1,N                                                       
      TEMP=A(I,RP1)                                                     
      A(I,RP1)=A(I,INTER)                                               
80    A(I,INTER)=TEMP                                                   
90    DO 100 I=RP2,N                                                    
      MULT(I)=A(I,R)/A(RP1,R)                                           
100   A(I,R)=MULT(I)                                                    
      DO 120 I=1,RP1                                                    
      TEMP=DCMPLX(0.0,0.0)                                              
      DO 110 J=RP2,N                                                    
110   TEMP=TEMP+A(I,J)*MULT(J)                                          
120   A(I,RP1)=A(I,RP1)+TEMP                                            
      DO 140 I=RP2,N                                                    
      TEMP=DCMPLX(0.0,0.0)                                              
      DO 130 J=RP2,N                                                    
130   TEMP=TEMP+A(I,J)*MULT(J)                                          
140   A(I,RP1)=A(I,RP1)+TEMP-MULT(I)*A(RP1,RP1)                         
      DO 150 I=RP2,N                                                    
      DO 150 J=RP2,N                                                    
150   A(I,J)=A(I,J)-MULT(I)*A(RP1,J)                                    
160   CONTINUE                                                          
C                                                                       
C        CALCULATE EPSILON                                              
C                                                                       
      EPS=ZERO                                                          
      DO 170 I=1,N                                                      
170   EPS=EPS+CDABS(A(1,I))                                             
      DO 190 I=2,N                                                      
      SUM=ZERO                                                          
      IM1=I-1                                                           
      DO 180 J=IM1,N                                                    
180   SUM=SUM+CDABS(A(I,J))                                             
      IF (SUM.GT.EPS) EPS=SUM                                           
190   CONTINUE                                                          
CC      FN=DBLE(FLOAT(N))                                               
      FN=DBLE(FLOAT(N))                                                 
      EPS=DSQRT(FN)*EPS*1.D-12                                          
      IF (EPS.EQ.ZERO)EPS=DBLE(1.0E-12)                                 
      DO 200 I=1,N                                                      
      DO 200 J=1,N                                                      
200   H(I,J)=A(I,J)                                                     
210   IF (N.NE.1) GO TO 220                                             
      LAMBDA(M)=A(1,1)+SHIFT(1)                                         
      GO TO 390                                                         
220   IF (N.EQ.2) GO TO 30                                              
230   MN1=M-N+1                                                         
      IF (DREAL(A(N,N)).NE.ZERO.OR.DIMAG(A(N,N)).NE.ZERO) IF(DABS(DREAL 
     1(A(N,N-1)/A(N,N)))+DABS(DIMAG(A(N,N-1)/A(N,N)))-1.D-9)            
     2 250,250,240                                                      
240   IF (DABS(DREAL(A(N,N-1)))+DABS(DIMAG(A(N,N-1))).GE.EPS) GO TO 260 
250   LAMBDA(MN1)=A(N,N)+SHIFT(1)                                       
      ICOUNT=0                                                          
      N=N-1                                                             
      GO TO 220                                                         
C                                                                       
C        DETERMINE SHIFT                                                
C                                                                       
260   SHIFT(2)=(A(N-1,N-1)+A(N,N)+CDSQRT((A(N-1,N-1)+A(N,N))**2-4.0*(A(N
     1,N)*A(N-1,N-1)-A(N,N-1)*A(N-1,N))))/2.0                           
      IF (DREAL(SHIFT(2)).NE.ZERO.OR.DIMAG(SHIFT(2)).NE.ZERO) GO TO 270 
      SHIFT(3)=A(N-1,N-1)+A(N,N)                                        
      GO TO 280                                                         
270   SHIFT(3)=(A(N,N)*A(N-1,N-1)-A(N,N-1)*A(N-1,N))/SHIFT(2)           
280   IF (CDABS(SHIFT(2)-A(N,N)).LT.CDABS(SHIFT(3)-A(N,N))) GO TO 290   
      INDEX=3                                                           
      GO TO 300                                                         
290   INDEX=2                                                           
300   IF (CDABS(A(N-1,N-2)).GE.EPS) GO TO 310                           
      LAMBDA(MN1)=SHIFT(2)+SHIFT(1)                                     
      LAMBDA(MN1+1)=SHIFT(3)+SHIFT(1)                                   
      ICOUNT=0                                                          
      N=N-2                                                             
      GO TO 210                                                         
310   SHIFT(1)=SHIFT(1)+SHIFT(INDEX)                                    
      DO 320 I=1,N                                                      
320   A(I,I)=A(I,I)-SHIFT(INDEX)                                        
C                                                                       
C      PERFORM GIVENS ROTATIONS, QR ITERATES. INCREASED NUMBER
C      OF INTERATIONS TO 13 ON 6/29/94 (R.B. EVANS, SAIC, NEW
C      LONDON, CT) 
C
      IF (ICOUNT.LE.13) GO TO 330                                       
      NCAL=M-N                                                          
      GO TO 390                                                         
330   NM1=N-1                                                           
      TEMP1=A(1,1)                                                      
      TEMP2=A(2,1)                                                      
      DO 380 R=1,NM1                                                    
      RP1=R+1                                                           
      RHO=DSQRT(DREAL(TEMP1)**2+DIMAG(TEMP1)**2+DREAL(TEMP2)**2+        
     1DIMAG(TEMP2)**2)                                                  
      IF (RHO.NE.ZERO) GO TO 340                                        
      TEMP1=A(RP1,RP1)                                                  
      IF (MAX-R.GE.2) TEMP2=A(R+2,R+1)                                  
      GO TO 380                                                         
340   COSS=TEMP1/RHO                                                    
      SINN=TEMP2/RHO                                                    
      INDEX=MAX0(R-1,1)                                                 
      DO 350 I=INDEX,N                                                  
      TEMP=DCONJG(COSS)*A(R,I)+DCONJG(SINN)*A(RP1,I)                    
      A(RP1,I)=-SINN*A(R,I)+COSS*A(RP1,I)                               
350   A(R,I)=TEMP                                                       
      TEMP1=A(RP1,RP1)                                                  
      IF (MAX-R.GE.2) TEMP2=A(R+2,R+1)                                  
      DO 360 I=1,R                                                      
      TEMP=COSS*A(I,R)+SINN*A(I,RP1)                                    
      A(I,RP1)=-DCONJG(SINN)*A(I,R)+DCONJG(COSS)*A(I,RP1)               
360   A(I,R)=TEMP                                                       
      INDEX=MIN0(R+2,N)                                                 
      DO 370 I=RP1,INDEX                                                
      A(I,R)=SINN*A(I,RP1)                                              
370   A(I,RP1)=DCONJG(COSS)*A(I,RP1)                                    
380   CONTINUE                                                          
      ICOUNT=ICOUNT+1                                                   
      GO TO 230                                                         
C                                                                       
C        VALUES IN INCREASING ORDER OF REAL PARTS                       
C
390   CONTINUE                                                          
C
      DO 400 I=1,NCAL                                                   
      K=NCAL+1-I                                                        
      VECT(I,1)=LAMBDA(K)
400   CONTINUE                                                          
      L=NCAL-1                                                          
      DO 420 I=1,L                                                      
      K=I+1                                                             
      DO 410 J=K,NCAL                                                   
      IF (DREAL(VECT(J,1)).GE.DREAL(VECT(I,1))) GO TO 410               
      TEMP=VECT(J,1)                                                    
      VECT(J,1)=VECT(I,1)                                               
      VECT(I,1)=TEMP                                                    
410   CONTINUE                                                          
420   LAMBDA(I)=VECT(I,1)                                               
      LAMBDA(NCAL)=VECT(NCAL,1)                                         
      N=M                                                               
      NM1=N-1                                                           
      IF (N.NE.2) GO TO 430                                             
      EPS=DMAX1(CDABS(LAMBDA(1)),CDABS(LAMBDA(2)))*1.D-8                
      IF (EPS.EQ.ZERO) EPS=DBLE(1.E-12)                                 
      H(1,1)=A(1,1)                                                     
      H(1,2)=A(1,2)                                                     
      H(2,1)=A(2,1)                                                     
      H(2,2)=A(2,2)                                                     
430   DO 440 I=1,M                                                      
      DO 440 J=1,M                                                      
440   A(J,I)=HL(J,I)                                                    
C                                                                       
C        CALCULATE VECTORS                                              
C                                                                       
      NC=MIN0(NCAL,NV)                                                  
      IF (NC.EQ.0) GO TO 640                                            
      DO 630 L=1,NC                                                     
      DO 460 I=1,N                                                      
      DO 450 J=1,N                                                      
450   HL(I,J)=H(I,J)                                                    
460   HL(I,I)=HL(I,I)-LAMBDA(L)                                         
      DO 500 I=1,NM1                                                    
      MULT(I)=DCMPLX(0.0,0.0)                                           
      INTH(I)=0                                                         
      IP1=I+1                                                           
      IF (CDABS(HL(I+1,I)).LE.CDABS(HL(I,I))) GO TO 480                 
      INTH(I)=1                                                         
      DO 470 J=I,N                                                      
      TEMP=HL(I+1,J)                                                    
      HL(I+1,J)=HL(I,J)                                                 
470   HL(I,J)=TEMP                                                      
480   IF (DREAL(HL(I,I)).EQ.ZERO.AND.DIMAG(HL(I,I)).EQ.ZERO) GO TO 500  
      MULT(I)=-HL(I+1,I)/HL(I,I)                                        
      DO 490 J=IP1,N                                                    
490   HL(I+1,J)=HL(I+1,J)+MULT(I)*HL(I,J)                               
500   CONTINUE                                                          
      DO 510 I=1,N                                                      
510   VECT(I,L)=DCMPLX(1.0,0.0)                                         
      TWICE=.FALSE.                                                     
520   IF (DREAL(HL(N,N)).EQ.ZERO.AND.DIMAG(HL(N,N)).EQ.ZERO)            
     1  HL(N,N)=DCMPLX(EPS,ZERO)                                        
      VECT(N,L)=VECT(N,L)/HL(N,N)                                       
      DO 540 I=1,NM1                                                    
      K=N-I                                                             
      DO 530 J=K,NM1                                                    
530   VECT(K,L)=VECT(K,L)-HL(K,J+1)*VECT(J+1,L)                         
      IF (DREAL(HL(K,K)).EQ.ZERO.AND.DIMAG(HL(K,K)).EQ.ZERO) HL(K,K)=   
     1  DCMPLX(EPS,ZERO)                                                
540   VECT(K,L)=VECT(K,L)/HL(K,K)                                       
      BIG=ZERO                                                          
      DO 550 I=1,N                                                      
      SUM=DABS(DREAL(VECT(I,L)))+DABS(DIMAG(VECT(I,L)))                 
      IF (SUM.GT.BIG) BIG=SUM                                           
550   CONTINUE                                                          
      DO 560 I=1,N                                                      
560   VECT(I,L)=VECT(I,L)/BIG                                           
      IF (TWICE) GO TO 580                                              
      DO 570 I=1,NM1                                                    
      IF (INTH(I).EQ.0) GO TO 570                                       
      TEMP=VECT(I,L)                                                    
      VECT(I,L)=VECT(I+1,L)                                             
      VECT(I+1,L)=TEMP                                                  
570   VECT(I+1,L)=VECT(I+1,L)+MULT(I)*VECT(I,L)                         
      TWICE=.TRUE.                                                      
      GO TO 520                                                         
580   IF (N.EQ.2) GO TO 610                                             
      NM2=N-2                                                           
      DO 600 I=1,NM2                                                    
      N1I=N-1-I                                                         
      NI1=N-I+1                                                         
      DO 590 J=NI1,N                                                    
590   VECT(J,L)=H(J,N1I)*VECT(N1I+1,L)+VECT(J,L)                        
      INDEX=INT(N1I)                                                    
      TEMP=VECT(N1I+1,L)                                                
      VECT(N1I+1,L)=VECT(INDEX,L)                                       
600   VECT(INDEX,L)=TEMP                                                
C                                                                       
C        NORMALIZE VECTORS TO UNIT LENGTH                               
C                                                                       
610   CONTINUE                                                          
      TEMP=DCMPLX(0.0,0.0)                                              
      DO 620 I=1,N                                                      
620   TEMP=VECT(I,L)*VECT(I,L)+TEMP                                     
      TEMP=CDSQRT(TEMP)                                                 
      DO 630 I=1,N                                                      
630   VECT(I,L)=VECT(I,L)/TEMP                                          
640   RETURN                                                            
      END                                                               
