spectral.for

Go to the documentation of this file.

      SUBROUTINE spectral(A,D,V)
C>    EIGENVALUES AND EIGENVECTOR OF A 3X3 MATRIX
C     THIS SUBROUTINE CALCULATES THE EIGENVALUES AND EIGENVECTORS OF
C     A SYMMETRIC 3X3 MATRIX A.
C
C     THE OUTPUT CONSISTS OF A VECTOR D CONTAINING THE THREE
C     EIGENVALUES IN ASCENDING ORDER, AND A MATRIX V WHOSE
C     COLUMNS CONTAIN THE CORRESPONDING EIGENVECTORS.
C
      IMPLICIT NONE
C
      INTEGER NP,NROT
      parameter(np=3)
C
      DOUBLE PRECISION D(3),V(3,3),A(3,3),E(3,3)
C
      e = a
C
      CALL jacobi(e,3,np,d,v,nrot)
      CALL eigsrt(d,v,3,np)
C
      RETURN
      END SUBROUTINE spectral

C***********************************************************************

      SUBROUTINE jacobi(A,N,NP,D,V,NROT)
C
C COMPUTES ALL EIGENVALUES AND EIGENVECTORS OF A REAL SYMMETRIC
C  MATRIX A, WHICH IS OF SIZE N BY N, STORED IN A PHYSICAL
C  NP BY NP ARRAY.  ON OUTPUT, ELEMENTS OF A ABOVE THE DIAGONAL
C  ARE DESTROYED, BUT THE DIAGONAL AND SUB-DIAGONAL ARE UNCHANGED
C  AND GIVE FULL INFORMATION ABOUT THE ORIGINAL SYMMETRIC MATRIX.
C  VECTOR D RETURNS THE EIGENVALUES OF A IN ITS FIRST N ELEMENTS.
C  V IS A MATRIX WITH THE SAME LOGICAL AND PHYSICAL DIMENSIONS AS
C  A WHOSE COLUMNS CONTAIN, UPON OUTPUT, THE NORMALIZED
C  EIGENVECTORS OF A.  NROT RETURNS THE NUMBER OF JACOBI ROTATION
C  WHICH WERE REQUIRED.
C
C THIS SUBROUTINE IS TAKEN FROM 'NUMERICAL RECIPES.'
C
      IMPLICIT NONE
      include 'PARAM_UMAT.INC'
C
      INTEGER IP,IQ,N,NMAX,NP,NROT,I,J
      parameter(nmax=100)
C
      DOUBLE PRECISION A(np,np),D(np),V(np,np),B(nmax),Z(nmax),
     +  sm,tresh,g,t,h,theta,s,c,tau


C INITIALIZE V TO THE IDENTITY MATRIX
      DO i=1,3
          v(i,i)=one
        DO j=1,3
          IF (i.NE.j)THEN
           v(i,j)=zero
         ENDIF
       END DO
      END DO
C INITIALIZE B AND D TO THE DIAGONAL OF A, AND Z TO ZERO.
C  THE VECTOR Z WILL ACCUMULATE TERMS OF THE FORM T*A_PQ AS
C  IN EQUATION (11.1.14)
C
      DO ip = 1,n
        b(ip) = a(ip,ip)
        d(ip) = b(ip)
        z(ip) = 0.d0
      END DO


C BEGIN ITERATION
C
      nrot = 0
      DO i=1,50
C
C         SUM OFF-DIAGONAL ELEMENTS
C
          sm = 0.d0
          DO ip=1,n-1
            DO iq=ip+1,n
              sm = sm + dabs(a(ip,iq))
            END DO
          END DO
C
C          IF SM = 0., THEN RETURN.  THIS IS THE NORMAL RETURN,
C          WHICH RELIES ON QUADRATIC CONVERGENCE TO MACHINE
C          UNDERFLOW.
C
          IF (sm.EQ.0.d0) RETURN
C
C          IN THE FIRST THREE SWEEPS CARRY OUT THE PQ ROTATION ONLY IF
C           |A_PQ| > TRESH, WHERE TRESH IS SOME THRESHOLD VALUE,
C           SEE EQUATION (11.1.25).  THEREAFTER TRESH = 0.
C
          IF (i.LT.4) THEN
            tresh = 0.2d0*sm/n**2
          ELSE
            tresh = 0.d0
          END IF
C
          DO ip=1,n-1
            DO iq=ip+1,n
              g = 100.d0*dabs(a(ip,iq))
C
C              AFTER FOUR SWEEPS, SKIP THE ROTATION IF THE
C               OFF-DIAGONAL ELEMENT IS SMALL.
C
              IF ((i.GT.4).AND.(dabs(d(ip))+g.EQ.dabs(d(ip)))
     +            .AND.(dabs(d(iq))+g.EQ.dabs(d(iq)))) THEN
                a(ip,iq) = 0.d0
              ELSE IF (dabs(a(ip,iq)).GT.tresh) THEN
                h = d(iq) - d(ip)
                IF (dabs(h)+g.EQ.dabs(h)) THEN
C
C                  T = 1./(2.*THETA), EQUATION (11.1.10)
C
                  t =a(ip,iq)/h
                ELSE
                  theta = 0.5d0*h/a(ip,iq)
                  t =1.d0/(dabs(theta)+dsqrt(1.d0+theta**2.d0))
                  IF (theta.LT.0.d0) t = -t
                END IF
                c = 1.d0/dsqrt(1.d0 + t**2.d0)
                s = t*c
                tau = s/(1.d0 + c)
                h = t*a(ip,iq)
                z(ip) = z(ip) - h
                z(iq) = z(iq) + h
                d(ip) = d(ip) - h
                d(iq) = d(iq) + h
                a(ip,iq) = 0.d0
C
C               CASE OF ROTATIONS 1 <= J < P
C
                DO j=1,ip-1
                  g = a(j,ip)
                  h = a(j,iq)
                  a(j,ip) = g - s*(h + g*tau)
                  a(j,iq) = h + s*(g - h*tau)
                END DO
C
C                CASE OF ROTATIONS P < J < Q
C
                DO j=ip+1,iq-1
                  g = a(ip,j)
                  h = a(j,iq)
                  a(ip,j) = g - s*(h + g*tau)
                  a(j,iq) = h + s*(g - h*tau)
                END DO
C
C                 CASE OF ROTATIONS Q < J <= N
C
                DO j=iq+1,n
                  g = a(ip,j)
                  h = a(iq,j)
                  a(ip,j) = g - s*(h + g*tau)
                  a(iq,j) = h + s*(g - h*tau)
                END DO
                DO j = 1,n
                  g = v(j,ip)
                  h = v(j,iq)
                  v(j,ip) = g - s*(h + g*tau)
                  v(j,iq) = h + s*(g - h*tau)
                END DO
                nrot = nrot + 1
             END IF
               END DO
             END DO
C
C          UPDATE D WITH THE SUM OF T*A_PQ, AND REINITIALIZE Z
C
       DO ip=1,n
         b(ip) = b(ip) + z(ip)
         d(ip) = b(ip)
         z(ip) = 0.d0
       END DO
      END DO
C
C IF THE ALGORITHM HAS REACHED THIS STAGE, THEN THERE
C  ARE TOO MANY SWEEPS.  PRINT A DIAGNOSTIC AND CUT THE
C  TIME INCREMENT.
C
      WRITE (*,'(/1X,A/)') '50 ITERATIONS IN JACOBI SHOULD NEVER HAPPEN'
C
      RETURN
      END SUBROUTINE jacobi

C**********************************************************************
      SUBROUTINE eigsrt(D,V,N,NP)
C
C     GIVEN THE EIGENVALUES D AND EIGENVECTORS V AS OUTPUT FROM
C     JACOBI, THIS SUBROUTINE SORTS THE EIGENVALUES INTO ASCENDING
C     ORDER AND REARRANGES THE COLMNS OF V ACCORDINGLY.
C
C     THE SUBROUTINE WAS TAKEN FROM 'NUMERICAL RECIPES.'
C
      IMPLICIT NONE
C
      INTEGER N,NP,I,J,K
C
      DOUBLE PRECISION D(np),V(np,np),P
C
      DO i=1,n-1
              k = i
              p = d(i)
              DO j=i+1,n
               IF (d(j).GE.p) THEN
                k = j
                p = d(j)
               END IF
              END DO
              IF (k.NE.i) THEN
               d(k) = d(i)
               d(i) = p
               DO j=1,n
                p = v(j,i)
                v(j,i) = v(j,k)
                v(j,k) = p
               END DO
              END IF
      END DO
C
      RETURN
      END SUBROUTINE eigsrt