Merge branch 'adasko' of mmka.chem.univ.gda.pl:unres into adasko
[unres.git] / source / unres / src_CSA / djacob.f
1       SUBROUTINE DJACOB(N,NMAX,MAXJAC,E,A,C,AII)                        
2       IMPLICIT REAL*8 (A-H,O-Z) 
3 C     THE JACOBI DIAGONALIZATION PROCEDURE
4       COMMON INP,IOUT,IPN                              
5       DIMENSION A(NMAX,N),C(NMAX,N),AII(150),AJJ(150)
6       SIN45 = .70710678                                                 
7       COS45 = .70710678                                                 
8       S45SQ = 0.50                                                      
9       C45SQ = 0.50                                                      
10 C     UNIT EIGENVECTOR MATRIX                                           
11       DO 70 I = 1,N                                                     
12       DO 7 J = I,N                                                      
13       A(J,I)=A(I,J)                                                     
14       C(I,J) = 0.0                                                      
15     7 C(J,I) = 0.0                                                      
16    70  C(I,I) = 1.0                                                     
17 C     DETERMINATION OF SEARCH ARGUMENT, TEST                            
18       AMAX = 0.0                                                        
19       DO 1 I = 1,N                                                      
20       DO 1 J = 1,I                                                      
21       TEMPA=DABS(A(I,J))                                                 
22       IF (AMAX-TEMPA) 2,1,1                                             
23     2 AMAX = TEMPA                                                      
24     1 CONTINUE                                                          
25       TEST = AMAX*E                                                     
26 C     SEARCH FOR LARGEST OFF DIAGONAL ELEMENT                           
27       DO 72 IJAC=1,MAXJAC                                               
28       AIJMAX = 0.0                                                      
29       DO 3 I = 2,N                                                      
30       LIM = I-1                                                         
31       DO 3 J = 1,LIM                                                    
32       TAIJ=DABS(A(I,J))                                                  
33        IF (AIJMAX-TAIJ) 4,3,3                                           
34     4 AIJMAX = TAIJ                                                     
35       IPIV = I                                                          
36       JPIV = J                                                          
37     3 CONTINUE                                                          
38       IF(AIJMAX-TEST)300,300,5                                          
39 C     PARAMETERS FOR ROTATION                                           
40     5 TAII = A(IPIV,IPIV)                                               
41       TAJJ = A(JPIV,JPIV)                                               
42       TAIJ = A(IPIV,JPIV)                                               
43       TMT = TAII-TAJJ                                                   
44       IF(DABS(TMT/TAIJ)-1.0D-12) 60,60,6                                 
45    60 IF(TAIJ) 10,10,11                                                 
46     6 ZAMMA=TAIJ/(2.0*TMT)                                              
47    90 IF(DABS(ZAMMA)-0.38268)8,8,9                                       
48     9 IF(ZAMMA)10,10,11                                                 
49    10 SINT = -SIN45                                                     
50       GO TO 12                                                          
51    11 SINT = SIN45                                                      
52    12 COST = COS45                                                      
53       SINSQ = S45SQ                                                     
54       COSSQ = C45SQ                                                     
55       GO TO 120                                                         
56     8 GAMSQ=ZAMMA*ZAMMA                                                 
57       SINT=2.0*ZAMMA/(1.0+GAMSQ)                                        
58       COST = (1.0-GAMSQ)/(1.0+GAMSQ)                                    
59       SINSQ=SINT*SINT                                                   
60       COSSQ=COST*COST                                                   
61 C     ROTATION                                                          
62   120 DO 13 K = 1,N                                                     
63       TAIK = A(IPIV,K)                                                  
64       TAJK = A(JPIV,K)                                                  
65       A(IPIV,K) = TAIK*COST+TAJK*SINT                                   
66       A(JPIV,K) = TAJK*COST-TAIK*SINT                                   
67       TCIK = C(IPIV,K)                                                  
68       TCJK = C(JPIV,K)                                                  
69       C(IPIV,K) = TCIK*COST+TCJK*SINT                                   
70    13 C(JPIV,K) = TCJK*COST-TCIK*SINT                                   
71       A(IPIV,IPIV) = TAII*COSSQ+TAJJ*SINSQ+2.0*TAIJ*SINT*COST           
72       A(JPIV,JPIV) = TAII*SINSQ+TAJJ*COSSQ-2.0*TAIJ*SINT*COST           
73       A(IPIV,JPIV) = TAIJ*(COSSQ-SINSQ)-SINT*COST*TMT                   
74       A(JPIV,IPIV) = A(IPIV,JPIV)                                       
75       DO 30 K = 1,N                                                     
76       A(K,IPIV) = A(IPIV,K)                                             
77    30 A(K,JPIV) = A(JPIV,K)                                             
78    72 CONTINUE                                                          
79       WRITE (IOUT,1000) AIJMAX                                             
80  1000 FORMAT (/1X,'NONCONVERGENT JACOBI. LARGEST OFF-DIAGONAL ELE',     
81      1 'MENT = ',1PE14.7)                                               
82 C     ARRANGEMENT OF EIGENVALUES IN ASCENDING ORDER                     
83   300 DO 14 I=1,N                                                       
84    14 AJJ(I)=A(I,I)                                                     
85       LT=N+1                                                            
86       DO15 L=1,N                                                        
87       LT=LT-1                                                           
88       AIIMIN=1.0E+30                                                    
89       DO16 I=1,N                                                        
90       IF(AJJ(I)-AIIMIN)17,16,16                                         
91    17 AIIMIN=AJJ(I)                                                     
92       IT=I                                                              
93    16 CONTINUE                                                          
94       IN=L                                                              
95       AII(IN)=AIIMIN                                                    
96       AJJ(IT)=1.0E+30                                                   
97       DO15 K=1,N                                                        
98    15 A(IN,K)=C(IT,K)                                                   
99       DO 18 I=1,N                                                       
100       IF(A(I,1))19,22,22                                                
101    19 T=-1.0                                                            
102       GO TO 91                                                          
103    22 T=1.0                                                             
104    91 DO 18 J=1,N                                                       
105    18 C(J,I)=T*A(I,J)                                                   
106       RETURN                                                            
107       END