Télécharger eig35.eso

Retour à la liste

Numérotation des lignes :

  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