SUBROUTINE DJACOB(N,NMAX,MAXJAC,E,A,C,AII)                        
      IMPLICIT REAL*8 (A-H,O-Z) 
C     THE JACOBI DIAGONALIZATION PROCEDURE
      COMMON INP,IOUT,IPN                              
      DIMENSION A(NMAX,N),C(NMAX,N),AII(150),AJJ(150)
      SIN45 = .70710678                                                 
      COS45 = .70710678                                                 
      S45SQ = 0.50                                                      
      C45SQ = 0.50                                                      
C     UNIT EIGENVECTOR MATRIX                                           
      DO 70 I = 1,N                                                     
      DO 7 J = I,N                                                      
      A(J,I)=A(I,J)                                                     
      C(I,J) = 0.0                                                      
    7 C(J,I) = 0.0                                                      
   70  C(I,I) = 1.0                                                     
C     DETERMINATION OF SEARCH ARGUMENT, TEST                            
      AMAX = 0.0                                                        
      DO 1 I = 1,N                                                      
      DO 1 J = 1,I                                                      
      TEMPA=DABS(A(I,J))                                                 
      IF (AMAX-TEMPA) 2,1,1                                             
    2 AMAX = TEMPA                                                      
    1 CONTINUE                                                          
      TEST = AMAX*E                                                     
C     SEARCH FOR LARGEST OFF DIAGONAL ELEMENT                           
      DO 72 IJAC=1,MAXJAC                                               
      AIJMAX = 0.0                                                      
      DO 3 I = 2,N                                                      
      LIM = I-1                                                         
      DO 3 J = 1,LIM                                                    
      TAIJ=DABS(A(I,J))                                                  
       IF (AIJMAX-TAIJ) 4,3,3                                           
    4 AIJMAX = TAIJ                                                     
      IPIV = I                                                          
      JPIV = J                                                          
    3 CONTINUE                                                          
      IF(AIJMAX-TEST)300,300,5                                          
C     PARAMETERS FOR ROTATION                                           
    5 TAII = A(IPIV,IPIV)                                               
      TAJJ = A(JPIV,JPIV)                                               
      TAIJ = A(IPIV,JPIV)                                               
      TMT = TAII-TAJJ                                                   
      IF(DABS(TMT/TAIJ)-1.0D-12) 60,60,6                                 
   60 IF(TAIJ) 10,10,11                                                 
    6 ZAMMA=TAIJ/(2.0*TMT)                                              
   90 IF(DABS(ZAMMA)-0.38268)8,8,9                                       
    9 IF(ZAMMA)10,10,11                                                 
   10 SINT = -SIN45                                                     
      GO TO 12                                                          
   11 SINT = SIN45                                                      
   12 COST = COS45                                                      
      SINSQ = S45SQ                                                     
      COSSQ = C45SQ                                                     
      GO TO 120                                                         
    8 GAMSQ=ZAMMA*ZAMMA                                                 
      SINT=2.0*ZAMMA/(1.0+GAMSQ)                                        
      COST = (1.0-GAMSQ)/(1.0+GAMSQ)                                    
      SINSQ=SINT*SINT                                                   
      COSSQ=COST*COST                                                   
C     ROTATION                                                          
  120 DO 13 K = 1,N                                                     
      TAIK = A(IPIV,K)                                                  
      TAJK = A(JPIV,K)                                                  
      A(IPIV,K) = TAIK*COST+TAJK*SINT                                   
      A(JPIV,K) = TAJK*COST-TAIK*SINT                                   
      TCIK = C(IPIV,K)                                                  
      TCJK = C(JPIV,K)                                                  
      C(IPIV,K) = TCIK*COST+TCJK*SINT                                   
   13 C(JPIV,K) = TCJK*COST-TCIK*SINT                                   
      A(IPIV,IPIV) = TAII*COSSQ+TAJJ*SINSQ+2.0*TAIJ*SINT*COST           
      A(JPIV,JPIV) = TAII*SINSQ+TAJJ*COSSQ-2.0*TAIJ*SINT*COST           
      A(IPIV,JPIV) = TAIJ*(COSSQ-SINSQ)-SINT*COST*TMT                   
      A(JPIV,IPIV) = A(IPIV,JPIV)                                       
      DO 30 K = 1,N                                                     
      A(K,IPIV) = A(IPIV,K)                                             
   30 A(K,JPIV) = A(JPIV,K)                                             
   72 CONTINUE                                                          
      WRITE (IOUT,1000) AIJMAX                                             
 1000 FORMAT (/1X,'NONCONVERGENT JACOBI. LARGEST OFF-DIAGONAL ELE',     
     1 'MENT = ',1PE14.7)                                               
C     ARRANGEMENT OF EIGENVALUES IN ASCENDING ORDER                     
  300 DO 14 I=1,N                                                       
   14 AJJ(I)=A(I,I)                                                     
      LT=N+1                                                            
      DO15 L=1,N                                                        
      LT=LT-1                                                           
      AIIMIN=1.0E+30                                                    
      DO16 I=1,N                                                        
      IF(AJJ(I)-AIIMIN)17,16,16                                         
   17 AIIMIN=AJJ(I)                                                     
      IT=I                                                              
   16 CONTINUE                                                          
      IN=L                                                              
      AII(IN)=AIIMIN                                                    
      AJJ(IT)=1.0E+30                                                   
      DO15 K=1,N                                                        
   15 A(IN,K)=C(IT,K)                                                   
      DO 18 I=1,N                                                       
      IF(A(I,1))19,22,22                                                
   19 T=-1.0                                                            
      GO TO 91                                                          
   22 T=1.0                                                             
   91 DO 18 J=1,N                                                       
   18 C(J,I)=T*A(I,J)                                                   
      RETURN                                                            
      END