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(k)=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 :
© Cast3M 2003 - All rights reserved.
Disclaimer