Télécharger dvdiag.eso

Retour à la liste

Numérotation des lignes :

  1. C DVDIAG SOURCE KICH 19/09/26 21:15:05 10311
  2. subroutine dvdiag(a,d)
  3. c
  4. c ***************
  5. c diagonalisation
  6. c ***************
  7. c
  8. implicit real*8 (a-h,o-z)
  9. implicit integer (i-n)
  10. dimension sig(6),a(3,3),d(3),e(3)
  11. n = 3
  12. np = 3
  13. do 18 i=n,2,-1
  14. l=i-1
  15. h=0.
  16. scale=0.
  17. if(l.gt.1)then
  18. do 11 k=1,l
  19. scale=scale+abs(a(i,k))
  20. 11 continue
  21. if(scale.eq.0.)then
  22. e(i)=a(i,l)
  23. else
  24. do 12 k=1,l
  25. a(i,k)=a(i,k)/scale
  26. h=h+a(i,k)**2
  27. 12 continue
  28. f=a(i,l)
  29. g=-sign(sqrt(h),f)
  30. e(i)=scale*g
  31. h=h-f*g
  32. a(i,l)=f-g
  33. f=0.
  34. do 15 j=1,l
  35. a(j,i)=a(i,j)/h
  36. g=0.
  37. do 13 k=1,j
  38. g=g+a(j,k)*a(i,k)
  39. 13 continue
  40. if(l.gt.j)then
  41. do 14 k=j+1,l
  42. g=g+a(k,j)*a(i,k)
  43. 14 continue
  44. endif
  45. e(j)=g/h
  46. f=f+e(j)*a(i,j)
  47. 15 continue
  48. hh=f/(h+h)
  49. do 17 j=1,l
  50. f=a(i,j)
  51. g=e(j)-hh*f
  52. e(j)=g
  53. do 16 k=1,j
  54. a(j,k)=a(j,k)-f*e(k)-g*a(i,k)
  55. 16 continue
  56. 17 continue
  57. endif
  58. else
  59. e(i)=a(i,l)
  60. endif
  61. d(i)=h
  62. 18 continue
  63. d(1)=0.
  64. e(1)=0.
  65. do 23 i=1,n
  66. l=i-1
  67. if(d(i).ne.0.)then
  68. do 21 j=1,l
  69. g=0.
  70. do 19 k=1,l
  71. g=g+a(i,k)*a(k,j)
  72. 19 continue
  73. do 20 k=1,l
  74. a(k,j)=a(k,j)-g*a(k,i)
  75. 20 continue
  76. 21 continue
  77. endif
  78. d(i)=a(i,i)
  79. a(i,i)=1.
  80. if(l.ge.1)then
  81. do 22 j=1,l
  82. a(i,j)=0.
  83. a(j,i)=0.
  84. 22 continue
  85. endif
  86. 23 continue
  87. do 911 i=2,n
  88. e(i-1)=e(i)
  89. 911 continue
  90. e(n)=0.
  91. do 915 l=1,n
  92. iter=0
  93. 91 do 912 m=l,n-1
  94. dd=abs(d(m))+abs(d(m+1))
  95. if (abs(e(m))+dd.eq.dd) go to 92
  96. 912 continue
  97. m=n
  98. 92 if(m.ne.l)then
  99. if(iter.eq.30) then
  100. write(6,*) 'too many iterations'
  101. call erreur(5)
  102. return
  103. endif
  104. iter=iter+1
  105. g=(d(l+1)-d(l))/(2.*e(l))
  106. r=sqrt(g**2+1.)
  107. g=d(m)-d(l)+e(l)/(g+sign(r,g))
  108. s=1.
  109. c=1.
  110. p=0.
  111. do 914 i=m-1,l,-1
  112. f=s*e(i)
  113. b=c*e(i)
  114. if(abs(f).ge.abs(g))then
  115. c=g/f
  116. r=sqrt(c**2+1.)
  117. e(i+1)=f*r
  118. s=1./r
  119. c=c*s
  120. else
  121. s=f/g
  122. r=sqrt(s**2+1.)
  123. e(i+1)=g*r
  124. c=1./r
  125. s=s*c
  126. endif
  127. g=d(i+1)-p
  128. r=(d(i)-g)*s+2.*c*b
  129. p=s*r
  130. d(i+1)=g+p
  131. g=c*r-b
  132. do 913 k=1,n
  133. f=a(k,i+1)
  134. a(k,i+1)=s*a(k,i)+c*f
  135. a(k,i)=c*a(k,i)-s*f
  136. 913 continue
  137. 914 continue
  138. d(l)=d(l)-p
  139. e(l)=g
  140. e(m)=0.
  141. go to 91
  142. endif
  143. 915 continue
  144. return
  145. end
  146.  
  147.  
  148.  

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