Sources Esope | Cast3M

b3_jac
C B3_JAC    SOURCE    PV090527  23/01/27    21:15:07     11574          
* ----------------------------------------------------------------------------
* Numerical diagonalization of 3x3 matrcies
* Copyright (C) 2006  Joachim Kopp
* ----------------------------------------------------------------------------
* This library is free software; you can redistribute it and/or
* modify it under the terms of the GNU Lesser General Public
* License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
*
* This library is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
* Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with this library; if not, write to the Free Software
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
* ----------------------------------------------------------------------------
 
 
* ----------------------------------------------------------------------------
      subroutine b3_jac(A,Q,W)
* ----------------------------------------------------------------------------
* Calculates the eigenvalues and normalized eigenvectors of a symmetric 3x3
* matrix A using the Jacobi algorithm.
* The upper triangular part of A is destroyed during the calculation,
* the diagonal elements are read but not destroyed, and the lower
* triangular elements are not referenced at all.
* ----------------------------------------------------------------------------
* Parameters:
*   A: The symmetric input matrix
*   Q: Storage buffer for eigenvectors
*   W: Storage buffer for eigenvalues
* ----------------------------------------------------------------------
*     .. Arguments ..
*      implicit none
      real*8 A(3,3),W(3),Q(3,3)
      real*8 B(3,3)
*     .. Parameters ..
      integer N
      parameter (N=3)
*     .. Local Variables ..
      real*8 SD, SO
      real*8 S, C, T
      real*8 G, H, Z, THETA
      real*8 THRESH
      integer I, X, Y, R
*     determinant
      real*8 det 
 
*     on verifie que le determinant ne soit pas nul
      det = a(1,1) * a(2,2) * a(3,3) - a(1,1) * a(2,3) * a(3,2) - a(1,2)
     # * a(2,1) * a(3,3) + a(1,2) * a(2,3) * a(3,1) + a(1,3) * a(2,1) * 
     #a(3,2) - a(1,3) * a(2,2) * a(3,1)
c      print*,'Dans b3_jac det=',det
      if(det.eq.0.d0) then
c        on ne calcul pas, on assigne
         do i=1,3
             W(i)=0.d0
             do j=1,3
                if(i.eq.j) then
                    Q(j,i)=1.d0
                else
                    Q(i,j)=0.d0
                end if
             end do
         end do
         return
      end if
 
 
*     Initialize Q to the identitity matrix
*     --- This loop can be omitted if only the eigenvalues are desired -
*      A=B
      do x=1,n
       do y=1,n
        b(x,y)=a(x,y)
       end do
      end do
      DO 10 X = 1, N
        Q(X,X) = 1.0D0
        DO 11, Y = 1, X-1
          Q(X, Y) = 0.0D0
          Q(Y, X) = 0.0D0
   11   CONTINUE
   10 CONTINUE
 
*     Initialize W to diag(A)
      DO 20 X = 1, N
        W(X) = A(X, X)
   20 CONTINUE
 
*     Calculate SQR(tr(A))  
      SD = 0.0D0
      DO 30 X = 1, N
        SD = SD + ABS(W(X))
   30 CONTINUE
      SD = SD**2
 
*     Main iteration loop
      DO 40 I = 1, 50
*       Test for convergence
        SO = 0.0D0
        DO 50 X = 1, N
          DO 51 Y = X+1, N
            SO = SO + ABS(A(X, Y))
   51     CONTINUE
   50   CONTINUE
        IF (SO .EQ. 0.0D0) THEN
c         print*,'convergence jacobi en ',I,' iterations'
          RETURN
        END IF
 
        IF (I .LT. 4) THEN
          THRESH = 0.2D0 * SO / N**2
        ELSE
          THRESH = 0.0D0
        END IF
 
*       Do sweep
        DO 60 X = 1, N
          DO 61 Y = X+1, N
            G = 100.0D0 * ( ABS(A(X, Y)) )
            IF ( I .GT. 4 .AND. (ABS(W(X)) + G) .EQ. ABS(W(X))
     &                    .AND. ABS(W(Y)) + G .EQ. ABS(W(Y)) ) THEN
              A(X, Y) = 0.0D0
            ELSE IF (ABS(A(X, Y)) .GT. THRESH) THEN
*             Calculate Jacobi transformation
              H = W(Y) - W(X)
              IF ( ABS(H) + G .EQ. ABS(H) ) THEN
                T = A(X, Y) / H
              ELSE
                THETA = 0.5D0 * H / A(X, Y)
                IF (THETA .LT. 0.0D0) THEN
                  T = -1.0D0 / (SQRT(1.0D0 + THETA**2) - THETA)
                ELSE
                  T = 1.0D0 / (SQRT(1.0D0 + THETA**2) + THETA)
                END IF
              END IF
 
              C = 1.0D0 / SQRT( 1.0D0 + T**2 )
              S = T * C
              Z = T * A(X, Y)
 
*             Apply Jacobi transformation
              A(X, Y) = 0.0D0
              W(X)    = W(X) - Z
              W(Y)    = W(Y) + Z
              DO 70 R = 1, X-1
                T       = A(R, X)
                A(R, X) = C * T - S * A(R, Y)
                A(R, Y) = S * T + C * A(R, Y)
   70         CONTINUE
              DO 80, R = X+1, Y-1
                T       = A(X, R)
                A(X, R) = C * T - S * A(R, Y)
                A(R, Y) = S * T + C * A(R, Y)
   80         CONTINUE
              DO 90, R = Y+1, N
                T       = A(X, R)
                A(X, R) = C * T - S * A(Y, R)
                A(Y, R) = S * T + C * A(Y, R)
   90         CONTINUE
 
*             Update eigenvectors
*             --- This loop can be omitted if only the eigenvalues are desired ---
              DO 100, R = 1, N
                T       = Q(R, X)
                Q(R, X) = C * T - S * Q(R, Y)
                Q(R, Y) = S * T + C * Q(R, Y)
  100         CONTINUE
            END IF
   61     CONTINUE
   60   CONTINUE
   40 CONTINUE
 
      PRINT *, "B3_JAC: No convergence."
      print*,'A'
      call afic33(B)
      print*,'W'
      do i=1,3
        print*,W(i)
      end do
      print*,'Q'
      call afic33(Q)
      return      
      END 
* End of subroutine DSYEVJ3