Esope source files

   1 : C EIG35     SOURCE    CHAT      05/01/12    23:28:43     5004           
   2 : c**************************************************************************
   3 :       SUBROUTINE EIG35(V,D,ROT)
   4 : c**************************************************************************
   5 :       IMPLICIT INTEGER(I-N)
   6 :       integer rot, its, i,j,k
   7 :       real*8  g,h,aij,sm,thresh,t,c,s,tau,v(3,3),d(3),a(3),b(3),z(3)
   8 :       a(1) = v(1,2)
   9 :       a(2) = v(2,3)
  10 :       a(3) = v(3,1)
  11 :       do i = 1,3
  12 :         d(i) = v(i,i)
  13 :         b(i) = d(i)
  14 :         z(i) = 0.0d0
  15 :         do j = 1,3
  16 :           v(i,j) = 0.0d0
  17 :         end do
  18 :         v(i,i) = 1.d0
  19 :       end do
  20 :       rot = 0
  21 :       do its = 1,50
  22 : c.... set convergence test and threshold
  23 :         sm = abs(a(1)) + abs(a(2)) + abs(a(3))
  24 :         if (sm.eq.0.0d0) return
  25 :         if(its.lt.4) then
  26 :           thresh = 0.011d0*sm
  27 :         else
  28 :           thresh = 0.0d0
  29 :         end if
  30 : c.... perform sweeps for rotations
  31 :         do i = 1,3
  32 :           j = mod(i,3) + 1
  33 :           k = mod(j,3) + 1
  34 :           aij  = a(i)
  35 :           g    = 100.d0*abs(aij)
  36 :           if(abs(d(i))+g.ne.abs(d(i)) .or.
  37 :      1       abs(d(j))+g.ne.abs(d(j))) then
  38 :             if(abs(aij).gt.thresh) then
  39 :               a(i) = 0.0d0
  40 :               h    = d(j) - d(i)
  41 :               if(abs(h)+g.eq.abs(h)) then
  42 :                 t = aij/h
  43 :               else
  44 :                 t = sign(2.d0,h/aij)/(abs(h/aij)+sqrt(4.d0+(h/aij)**2))
  45 :               endif
  46 : c.... set rotation parameters
  47 :               c    = 1.d0/sqrt(1.d0+t*t)
  48 :               s    = t*c
  49 :               tau  = s/(1.d0+c)
  50 : c.... rotate diagonal terms
  51 :               h    = t*aij
  52 :               z(i) = z(i) - h
  53 :               z(j) = z(j) + h
  54 :               d(i) = d(i) - h
  55 :               d(j) = d(j) + h
  56 : c.... rotate off-diagonal terms
  57 :               h    = a(j)
  58 :               g    = a(k)
  59 :               a(j) = h + s*(g - h*tau)
  60 :               a(k) = g - s*(h + g*tau)
  61 : c.... rotate eigenvectors
  62 :               do k = 1,3
  63 :                 g      = v(k,i)
  64 :                 h      = v(k,j)
  65 :                 v(k,i) = g - s*(h + g*tau)
  66 :                 v(k,j) = h + s*(g - h*tau)
  67 :               end do 
  68 :               rot = rot + 1
  69 :             end if
  70 :           else
  71 :             a(i) = 0.0d0
  72 :           end if
  73 :         end do 
  74 : c.... update the diagonal terms
  75 :         do i = 1,3
  76 :           b(i) = b(i) + z(i)
  77 :           d(i) = b(i)
  78 :           z(i) = 0.0d0
  79 :         end do 
  80 :       end do 
  81 :       end
  82 :  
  83 :