Télécharger la plaquette Cast3M


fr en

Télécharger eig35.eso

Retour à la liste

Numérotation des lignes : 

eig35

  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.  
  84.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales