59 lines
1.3 KiB
Fortran
59 lines
1.3 KiB
Fortran
SUBROUTINE laguer(a,m,x,its)
|
|
C ============================
|
|
C
|
|
C Routine from Numerical Recipees
|
|
C
|
|
INCLUDE 'IMPLIC.FOR'
|
|
COMPLEX*16 a(m+1),x
|
|
PARAMETER (EPSS=2.e-7,MR=8,MT=10,MAXIT=MT*MR)
|
|
REAL frac(MR)
|
|
COMPLEX*16 dx,x1,b,d,f,g,h,sq,gp,gm,g2
|
|
SAVE frac
|
|
DATA frac /.5,.25,.75,.13,.38,.62,.88,1./
|
|
C
|
|
do iter=1,MAXIT
|
|
its=iter
|
|
b=a(m+1)
|
|
err=abs(b)
|
|
d=dcmplx(0.d0,0.d0)
|
|
f=dcmplx(0.d0,0.d0)
|
|
abx=abs(x)
|
|
do j=m,1,-1
|
|
f=x*f+d
|
|
d=x*d+b
|
|
b=x*b+a(j)
|
|
err=abs(b)+abx*err
|
|
end do
|
|
err=EPSS*err
|
|
if(abs(b).le.err) then
|
|
return
|
|
else
|
|
g=d/b
|
|
g2=g*g
|
|
h=g2-2.d0*f/b
|
|
sq=sqrt(dble(m-1)*(dble(m)*h-g2))
|
|
gp=g+sq
|
|
gm=g-sq
|
|
abp=abs(gp)
|
|
abm=abs(gm)
|
|
if(abp.lt.abm) gp=gm
|
|
if (max(abp,abm).gt.0.D0) then
|
|
dx=dble(m)/gp
|
|
else
|
|
dx=exp(dcmplx(log(1.d0+abx),dble(iter)))
|
|
endif
|
|
endif
|
|
x1=x-dx
|
|
if(x.eq.x1)return
|
|
if (mod(iter,MT).ne.0) then
|
|
x=x1
|
|
else
|
|
x=x-dx*frac(iter/MT)
|
|
endif
|
|
end do
|
|
c
|
|
write(6,601) x,x1
|
|
601 format(' too many iterations in laguer, x,x1 ',1p2e9.1)
|
|
return
|
|
END
|