Sources Esope | Cast3M

zroots
C ZROOTS    SOURCE    CHAT      05/01/13    04:24:56     5004
      subroutine zroots (a,m,roots,ad,coef,racr,raci,kerre)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
      dimension racr(*),raci(*),coef(*)
      complex a(*),roots(*)
      complex ad(*),x,b,c,xval
 
-INC PPARAM
-INC CCOPTIO
       eps=1.e-7
      maxm=m+1
      do 1 j=1,m+1
       ad(j)=a(j)
   1  continue
*       write(6,*) ' coef',(coef(i),i=1,17)
      j=m
      do 2 k=m,1,-1
      x = cmplx(0.,0.)
*      write(6,*) ' entree dans laguer'
      call laguer(ad,j,x,its,kerre)
 
      if(kerre.ne.0) return
* recherche de la valeur excate par newton et de la derivée
* calcul de la valeur
      xr=real(x)
      xi= aimag(x)
*      write(6,*) ' sortie de laguer x ',x
*      write(6,*) ' appel a zernew'
      call zernew(xr,xi,coef,m,xvr,xvi,kerre)
      if(kerre.ne.0) return
*      write(6,*) ' sortie de zernew xvr xvi',xvr,xvi
      if(abs(xvi).le.2.*eps**2*abs(xvr)) xvi=0.d0
      x = cmplx(xvr,xvi)
      roots(j)=x
      racr(j)=xvr
      raci(j)=xvi
      b=ad(j+1)
      do 3 jj=j,1,-1
      c=ad(jj)
      ad(jj)=b
      b=x*b+c
   3  continue
      if(abs(xvi).gt.2.*eps**2*abs(xvr))then
*     if(abs(xvi).gt.eps*abs(xvr))then
* il existe surement la valeur complexe conjuguée
      j=j-1
      x=cmplx(real(x),-1.*aimag(x))
      roots(j)=x
      racr(j)=xvr
      raci(j)=-xvi
      b=ad(j+1)
      do 30 jj=j,1,-1
      c= ad(jj)
      ad(jj)=b
      b=x*b+c
   30 continue
      endif
      j=j-1
      if(j.eq.0) go to 21
   2  continue
   21 continue
      return
      end