Sources Esope | Cast3M

bairst
C BAIRST    SOURCE    CHAT      05/01/12    21:31:31     5004
      subroutine bairst(mlree1,iretr,ireti,ireel,kerre)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
C
C recherche des racines reel d'un polynome de degre N par la methode
C  de lin-bairstow ( procedure fournie par Mr Valbuena CERN)
C
-INC SMLREEL
      pointeur mlree4.mlreel,mlree5.mlreel
      segini,mlreel=mlree1
      kerre=0
      jg = prog(/1)
      segini mlree1
      segini mlree2
      segini mlree3
      segini mlree4
      if(ireel.eq.0) segini mlree5
      ndeg= prog(/1) - 1
      k = 0
      n=jg - 1
      do 10 ilj=1,10000
*      write(6,*) ' boucle 10 ilj', ilj
      r = prog(2) / prog(3)
      s = prog(1)/prog(3)
      rmin=abs (r)
       do 20 ijl=1,100000
*       write(6,*) ' boucle 20 ijl  n ', ijl,n
       mlree1.prog(n+1)=prog(n+1)
       mlree2.prog(n+1)=0.d0
       mlree3.prog(n+1)=0.d0
       mlree1.prog(n) = prog(n) - r * mlree1.prog(n+1)
      mlree2.prog(n) = -mlree1.prog(n+1)
      mlree3.prog(n) = 0.d0
      do 30 j=2,n-2
       mlree1.prog(n-j+1)=prog(n-j+1)-r*mlree1.prog(n-j+2)
     $ - s*mlree1.prog(n-j+3)
       mlree2.prog(n-j+1)=-mlree1.prog(n-j+2)-r*mlree2.prog(n-j+2)
     $ - s*mlree2.prog(n-j+3)
      mlree3.prog(n-j+1)= -mlree1.prog(n-j+3)-r*mlree3.prog(n-j+2)
     $ -s*mlree3.prog(n-j+3)
   30 continue
       r1=prog(2)- r * mlree1.prog(3) -s*mlree1.prog(4)
       s1=prog(1) - s*mlree1.prog(3)
       t=-mlree1.prog(3) -r*mlree2.prog(3) -s*mlree2.prog(4)
       u=-mlree1.prog(4)-s*mlree3.prog(4) -r*mlree3.prog(3)
       v=-s*mlree2.prog(3)
       w=-mlree1.prog(3) -s*mlree3.prog(3)
       r2=(s1*u-r1*w)/(t*w-u*v)
       s2=(v*r1-t*s1)/(t*w-u*v)
       s=s+s2
       r=r+r2
*       write(6,*) ' boucle 20 ijl  n r r2', ijl,n, r,r2
*       write(6,*) ' r2 ' , r2
       if( abs(r2).GT.0.0000001) then
*       write(6,*) ' boucle 20 ijl  n r r2', ijl,n, r,r2
          if(abs(r2).lt.rmin) then
            nite=ijl
            rmin=abs(r2)
            rvrai=r-r2
            svrai=s-s2
         endif
         go to 20
       else
*        write(6,*) ' boucle 20 ijl  n r r2', ijl,n, r,r2
         go to 21
       endif
   20  continue
        kerre=944
*       write(6,*) ' non convergence rmin  degré n ite' ,rmin,n,nite
       r=rvrai
       s=svrai
   21  continue
      do 40 klm=1,100000
*      write(6,*) ' klm  n k' , klm,n,k
*       write(6,*) ' klm  n ' , klm,n
       g=r*r-4*s
       if(g.gt.0.) then
         k = k + 1
         preel=-r/2.d0
         pima=sqrt(g) /2.d0
         mlree4.prog(k)= preel + pima
         k = k + 1
         mlree4.prog(k)=preel - pima
*        write(6,*) ' sol reel ' , mlree4.prog(k-1),mlree4.prog(k)
       else
         if(ireel.eq.0) then
          k = k+1
          preel = -r/2.d0
          pima = sqrt (-1.*g) / 2.
          mlree4.prog(k)=preel
          mlree5.prog(k)=pima
          k = k +1
          mlree4.prog(k)=preel
          mlree5.prog(k)=-pima
*        write(6,*) ' sol ima ' , mlree4.prog(k-1),mlree5.prog(k-1)
        endif
      endif
      n = n -2
      if(n.eq.0) then
       go to 11
      endif
       do 50 i=n,0,-1
*       write(6,*) ' boucle 50  i',i
       prog(i+1)=mlree1.prog(i+3)
   50  continue
       if(n.GT.2) go to 10
       if(n.eq.2) then
          r = prog(n)/prog(n+1)
          s = prog(n-1) / prog(n+1)
          go to 40
       endif
       go to 41
   40  continue
   41  continue
       if(n.lt.2) then
        if(prog(n+1).ne.0.)then
         k = k +1
         mlree4.prog(k)= -prog(n)/prog(n+1)
        endif
       endif
       go to 11
   10 continue
   11 continue
      jg = k
      segadj mlree4
      if(ireel.eq.0) segadj mlree5
      segsup mlree1,mlree2,mlree3,mlreel
      segdes mlree4
      iretr=mlree4
      if(ireel.eq.0) then
        segadj mlree5
        segdes mlree5
       ireti=mlree5
      endif
      return
      end