Esope source files

   1 : C LAGUER    SOURCE    PV        07/11/23    21:17:36     5978           
   2 :       subroutine laguer(a,m,x,its,kerre)
   3 :       IMPLICIT INTEGER(I-N)
   4 :       IMPLICIT REAL*8 (A-H,O-Z)
   5 :       dimension frac(8)
   6 :       real*8 err,abp,abm,abx,frac
   7 :       complex dx,x1,b,d,f,g,h,sq,gp,gm,g2
   8 :       complex a(m+1),x
   9 :       data frac/.5,.25,.75,.13,.38,.62,.88,1./
  10 : *      save frac
  11 : *      write(6,*) ' entree dans laguer m ',m
  12 :       epss = 2.e-7
  13 :       mr = 8
  14 :       mt=10
  15 :       maxit = mt*mr*3
  16 :       do 1 iter = 1,maxit
  17 :       its=iter
  18 :       b=a(m+1)
  19 :       err = abs(b)
  20 :       d=cmplx(0.,0.)
  21 :       f=cmplx(0.,0.) 
  22 :       abx = abs(x) 
  23 :       do 2 j=m,1,-1
  24 :           f = x*f+d
  25 :           d=x*d+b
  26 :           b=x*b+a(j)
  27 :           err=abs(b)+abx*err 
  28 :     2  continue
  29 :        err=epss*err
  30 :       if(abs(b).le.err) then
  31 :          return
  32 :       else
  33 :          g=d/b
  34 :          g2=g*g
  35 :          h=g2-2.*f/b
  36 :          sq=sqrt((m-1)*(m*h-g2))
  37 :          gp=g+sq
  38 :          gm=g-sq
  39 :          abp=abs(gp)
  40 :          abm=abs(gm)
  41 :          if(abp.lt.abm) gp=gm
  42 :          if(max(abp,abm).gt.0.) then
  43 :            dx=m/gp
  44 :          else
  45 :            dx=exp(cmplx(log(1.+abx),float(iter)))
  46 :          endif
  47 :       endif
  48 :       x1=x-dx
  49 :       if(x.eq.x1) return
  50 :       if(mod(iter,MT).ne.0) then
  51 :         x=x1
  52 :       else 
  53 :         x=x-dx*frac(iter/mt)
  54 :       endif
  55 :    1  continue
  56 :       kerre =944
  57 :       return
  58 :       end
  59 :   
  60 :  
  61 :  
  62 :  
  63 :  
  64 :