Télécharger bairst.eso

Retour à la liste

Numérotation des lignes :

  1. C BAIRST SOURCE CHAT 05/01/12 21:31:31 5004
  2. subroutine bairst(mlree1,iretr,ireti,ireel,kerre)
  3. IMPLICIT INTEGER(I-N)
  4. IMPLICIT REAL*8(A-H,O-Z)
  5. C
  6. C recherche des racines reel d'un polynome de degre N par la methode
  7. C de lin-bairstow ( procedure fournie par Mr Valbuena CERN)
  8. C
  9. -INC SMLREEL
  10. pointeur mlree4.mlreel,mlree5.mlreel
  11. segini,mlreel=mlree1
  12. kerre=0
  13. jg = prog(/1)
  14. segini mlree1
  15. segini mlree2
  16. segini mlree3
  17. segini mlree4
  18. if(ireel.eq.0) segini mlree5
  19. ndeg= prog(/1) - 1
  20. k = 0
  21. n=jg - 1
  22. do 10 ilj=1,10000
  23. * write(6,*) ' boucle 10 ilj', ilj
  24. r = prog(2) / prog(3)
  25. s = prog(1)/prog(3)
  26. rmin=abs (r)
  27. do 20 ijl=1,100000
  28. * write(6,*) ' boucle 20 ijl n ', ijl,n
  29. mlree1.prog(n+1)=prog(n+1)
  30. mlree2.prog(n+1)=0.d0
  31. mlree3.prog(n+1)=0.d0
  32. mlree1.prog(n) = prog(n) - r * mlree1.prog(n+1)
  33. mlree2.prog(n) = -mlree1.prog(n+1)
  34. mlree3.prog(n) = 0.d0
  35. do 30 j=2,n-2
  36. mlree1.prog(n-j+1)=prog(n-j+1)-r*mlree1.prog(n-j+2)
  37. $ - s*mlree1.prog(n-j+3)
  38. mlree2.prog(n-j+1)=-mlree1.prog(n-j+2)-r*mlree2.prog(n-j+2)
  39. $ - s*mlree2.prog(n-j+3)
  40. mlree3.prog(n-j+1)= -mlree1.prog(n-j+3)-r*mlree3.prog(n-j+2)
  41. $ -s*mlree3.prog(n-j+3)
  42. 30 continue
  43. r1=prog(2)- r * mlree1.prog(3) -s*mlree1.prog(4)
  44. s1=prog(1) - s*mlree1.prog(3)
  45. t=-mlree1.prog(3) -r*mlree2.prog(3) -s*mlree2.prog(4)
  46. u=-mlree1.prog(4)-s*mlree3.prog(4) -r*mlree3.prog(3)
  47. v=-s*mlree2.prog(3)
  48. w=-mlree1.prog(3) -s*mlree3.prog(3)
  49. r2=(s1*u-r1*w)/(t*w-u*v)
  50. s2=(v*r1-t*s1)/(t*w-u*v)
  51. s=s+s2
  52. r=r+r2
  53. * write(6,*) ' boucle 20 ijl n r r2', ijl,n, r,r2
  54. * write(6,*) ' r2 ' , r2
  55. if( abs(r2).GT.0.0000001) then
  56. * write(6,*) ' boucle 20 ijl n r r2', ijl,n, r,r2
  57. if(abs(r2).lt.rmin) then
  58. nite=ijl
  59. rmin=abs(r2)
  60. rvrai=r-r2
  61. svrai=s-s2
  62. endif
  63. go to 20
  64. else
  65. * write(6,*) ' boucle 20 ijl n r r2', ijl,n, r,r2
  66. go to 21
  67. endif
  68. 20 continue
  69. kerre=944
  70. * write(6,*) ' non convergence rmin degré n ite' ,rmin,n,nite
  71. r=rvrai
  72. s=svrai
  73. 21 continue
  74. do 40 klm=1,100000
  75. * write(6,*) ' klm n k' , klm,n,k
  76. * write(6,*) ' klm n ' , klm,n
  77. g=r*r-4*s
  78. if(g.gt.0.) then
  79. k = k + 1
  80. preel=-r/2.d0
  81. pima=sqrt(g) /2.d0
  82. mlree4.prog(k)= preel + pima
  83. k = k + 1
  84. mlree4.prog(k)=preel - pima
  85. * write(6,*) ' sol reel ' , mlree4.prog(k-1),mlree4.prog(k)
  86. else
  87. if(ireel.eq.0) then
  88. k = k+1
  89. preel = -r/2.d0
  90. pima = sqrt (-1.*g) / 2.
  91. mlree4.prog(k)=preel
  92. mlree5.prog(k)=pima
  93. k = k +1
  94. mlree4.prog(j<{pan style="color: #009900;">)=preel
  95. mlree5.prog(k)=-pima
  96. * write(6,*) ' sol ima ' , mlree4.prog(k-1),mlree5.prog(k-1)
  97. endif
  98. endif
  99. n = n -2
  100. if(n.eq.0) then
  101. go to 11
  102. endif
  103. do 50 i=n,0,-1
  104. * write(6,*) ' boucle 50 i',i
  105. prog(i+1)=mlree1.prog(i+3)
  106. 50 continue
  107. if(n.GT.2) go to 10
  108. if(n.eq.2) then
  109. r = prog(n)/prog(n+1)
  110. s = prog(n-1) / prog(n+1)
  111. go to 40
  112. endif
  113. go to 41
  114. 40 continue
  115. 41 continue
  116. if(n.lt.2) then
  117. if(prog(n+1).ne.0.)then
  118. k = k +1
  119. mlree4.prog(k)= -prog(n)/prog(n+1)
  120. endif
  121. endif
  122. go to 11
  123. 10 continue
  124. 11 continue
  125. jg = k
  126. segadj mlree4
  127. if(ireel.eq.0) segadj mlree5
  128. segsup mlree1,mlree2,mlree3,mlreel
  129. segdes mlree4
  130. iretr=mlree4
  131. if(ireel.eq.0) then
  132. segadj mlree5
  133. segdes mlree5
  134. ireti=mlree5
  135. endif
  136. return
  137. end
  138.  
  139.  
  140.  
  141.  
  142.  

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