Merge branch 'devel' into AFM

[unres.git] / source / unres / src_MD-NEWSC / djacob.f

diff --git a/source/unres/src_MD-NEWSC/djacob.f b/source/unres/src_MD-NEWSC/djacob.f
new file mode 100644 (file)
index 0000000..e3f46bc
--- /dev/null
+++ b/source/unres/src_MD-NEWSC/djacob.f
@@ -0,0 +1,107 @@
+      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