Fixed a_max expression (diag routine). Alternative formulation of jro…

…tate (conceptionally without local copy of one of the input vectors).
cp2k · Jun 26, 2019 · 92e01a7 · 92e01a7
1 parent 5a5c2c7
commit 92e01a7
Showing 1 changed file with 9 additions and 7 deletions.
diff --git a/src/common/mathlib.F b/src/common/mathlib.F
@@ -1575,7 +1575,7 @@ SUBROUTINE diag(n, a, d, v)
                                                             t, tau, theta, tresh
       REAL(KIND=dp), DIMENSION(n)                        :: b, z
 
-      a_max = a(n-1, n)
+      a_max = ABS(a(n-1, n))
       DO ip = 1, n-2
          a_max = MAX(a_max, MAXVAL(ABS(a(ip, ip+1:n))))
          b(ip) = a(ip, ip) ! get_diag(a)
@@ -1628,7 +1628,7 @@ SUBROUTINE diag(n, a, d, v)
             END DO
          END DO
          b = b+z
-         a_max = a(n-1, n)
+         a_max = ABS(a(n-1, n))
          DO ip = 1, n-2
             a_max = MAX(a_max, MAXVAL(ABS(a(ip, ip+1:n))))
          END DO
@@ -1651,12 +1651,13 @@ SUBROUTINE jrotate(a, b, ss, tt)
       REAL(KIND=dp), DIMENSION(:), INTENT(INOUT)         :: a, b
       REAL(KIND=dp), INTENT(IN)                          :: ss, tt
 
-      REAL(KIND=dp)                                      :: d
-      REAL(KIND=dp), DIMENSION(SIZE(a))                  :: c
+      REAL(KIND=dp)                                      :: u, v
 
-      c(:) = a(:); d = 1-ss*tt
-      a(:) = a(:)*d-b(:)*ss
-      b(:) = b(:)*d+c(:)*ss
+      u = 1.0_dp-ss*tt
+      v = ss/u
+
+      a = a*u-b*ss
+      b = b*(u+ss*v)+a*v
 
    END SUBROUTINE jrotate
 
@@ -1737,3 +1738,4 @@ SUBROUTINE swap_vector(a, b)
    END SUBROUTINE swap_vector
 
 END MODULE mathlib
+