Esope source files

   1 : C ZROOTS    SOURCE    CHAT      05/01/13    04:24:56     5004           
   2 :       subroutine zroots (a,m,roots,ad,coef,racr,raci,kerre)
   3 :       IMPLICIT INTEGER(I-N)
   4 :       IMPLICIT REAL*8(A-H,O-Z)
   5 :       dimension racr(*),raci(*),coef(*)
   6 :       complex a(*),roots(*)
   7 :       complex ad(*),x,b,c,xval
   8 : -INC CCOPTIO
   9 :        eps=1.e-7
  10 :       maxm=m+1
  11 :       do 1 j=1,m+1
  12 :        ad(j)=a(j)
  13 :    1  continue
  14 : *       write(6,*) ' coef',(coef(i),i=1,17)
  15 :       j=m
  16 :       do 2 k=m,1,-1
  17 :       x = cmplx(0.,0.)
  18 : *      write(6,*) ' entree dans laguer'
  18 :       call laguer(ad,j,x,its,kerre)  
  19 : 
  20 :       if(kerre.ne.0) return
  21 : * recherche de la valeur excate par newton et de la derivée
  22 : * calcul de la valeur
  23 :       xr=real(x)
  24 :       xi= aimag(x)
  25 : *      write(6,*) ' sortie de laguer x ',x
  26 : *      write(6,*) ' appel a zernew'
  26 :       call zernew(xr,xi,coef,m,xvr,xvi,kerre)
  27 :       if(kerre.ne.0) return
  28 : *      write(6,*) ' sortie de zernew xvr xvi',xvr,xvi
  29 :       if(abs(xvi).le.2.*eps**2*abs(xvr)) xvi=0.d0
  30 :       x = cmplx(xvr,xvi)
  31 :       roots(j)=x
  32 :       racr(j)=xvr
  33 :       raci(j)=xvi
  34 :       b=ad(j+1)
  35 :       do 3 jj=j,1,-1
  36 :       c=ad(jj)
  37 :       ad(jj)=b
  38 :       b=x*b+c
  39 :    3  continue
  40 :       if(abs(xvi).gt.2.*eps**2*abs(xvr))then
  41 : *     if(abs(xvi).gt.eps*abs(xvr))then
  42 : * il existe surement la valeur complexe conjuguée
  43 :       j=j-1
  44 :       x=cmplx(real(x),-1.*aimag(x))
  45 :       roots(j)=x
  46 :       racr(j)=xvr
  47 :       raci(j)=-xvi      
  48 :       b=ad(j+1)
  49 :       do 30 jj=j,1,-1
  50 :       c= ad(jj)
  51 :       ad(jj)=b
  52 :       b=x*b+c
  53 :    30 continue  
  54 :       endif 
  55 :       j=j-1
  56 :       if(j.eq.0) go to 21
  57 :    2  continue
  58 :    21 continue
  59 :       return
  60 :       end
  61 :  
  62 : 
  63 : 
  64 : 
  65 : 
  66 : 
  67 :  
  68 :  
  69 :