subroutine cheb_r4 ( n, a, b, fnu, alpha, beta, s, ierr, s0, s1, & s2 ) c*********************************************************************72 c cc CHEB_R4 generates recursion coefficients ALPHA and BETA. c c Given a set of polynomials p(0),p(1),...,p(2*n-1) satisfying c c p(k+1)(x)=(x-a(k))*p(k)(x)-b(k)*p(k-1)(x), c k=0,1,...,2*n-2, c c p(-1)(x)=0, p(0)(x)=1, c c and associated modified moments c c fnu(k)=integral of p(k)(x)*dlambda(x), c k=0,1,...,2*n-1, c c this subroutine uses the modified Chebyshev algorithm to generate the c recursion coefficients alpha(k),beta(k), k=0,1,...,n-1, for the c polynomials pi(k) orthogonal with respect to the integration c measure dlambda(x), i.e., c c pi(k+1)(x)=(x-alpha(k))*pi(k)(x)-beta(k)*pi(k-1)(x), c k=0,1,...,n-1, c c pi(-1)(x)=0, pi(0)(x)=1. c c Reference: c c Walter Gautschi, c On Generating Orthogonal Polynomials, c SIAM Journal on Scientific and Statistical Computing, c Volume 3, Number 3, 1982, pages 289-317. c c Input: n - - the number of recursion coefficients desired c a,b-- arrays of dimension 2*n-1 to be filled with the c values of a(k-1),b(k-1), k=1,2,...,2*n-1 c fnu-- array of dimension 2*n to be filled with the c values of the modified moments fnu(k-1), k=1,2, c ...,2*n c Output: alpha,beta-- arrays containing, respectively, the c recursion coefficients alpha(k-1),beta(k-1), c k=1,2,...,n, where beta(0) is the total mass. c s - - array containing the normalization factors c s(k)=integral [pi(k)(x)]**2 dlambda(x), k=0,1, c 2,...,n-1. c ierr- an error flag, equal to 0 on normal return, c equal to 1 if abs(fnu(0)) is less than the c machine zero, equal to 2 if n is out of range, c equal to -k if s(k), k=0,1,2,...,n-1, is about c to underflow, and equal to +k if it is about to c overflow. c c The arrays s0,s1,s2 are needed for working space. c c On machines with limited exponent range, the occurrence of underflow c [overflow] in the computation of the alpha's and beta's can often c be avoided by multiplying all modified moments by a sufficiently large c [small] scaling factor and dividing the new beta(0) by the same c scaling factor. c c The routine uses the function subroutine r1mach. c dimension a(*),b(*),fnu(*),alpha(n),beta(n),s(n),s0(*),s1(*), *s2(*) c c The arrays a,b are assumed to have dimension 2*n-1, the arrays c fnu,s0,s1,s2 dimension 2*n. c nd=2*n tiny=10.0*r1mach(1) huge=0.1*r1mach(2) ierr=0 if(abs(fnu(1)).lt.tiny) then ierr=1 return end if if(n.lt.1) then ierr=2 return end if c c Initialization c alpha(1)=a(1)+fnu(2)/fnu(1) beta(1)=fnu(1) if(n.eq.1) return s(1)=fnu(1) do l=1,nd s0(l)=0.0 s1(l)=fnu(l) end do c c Continuation c do k=2,n lk=nd-k+1 do l=k,lk c c The quantities s2(l) for l > k are auxiliary quantities which may c be zero or may become so small as to underflow, without however c causing any harm. c s2(l)=s1(l+1)-(alpha(k-1)-a(l))*s1(l)-beta(k-1)*s0(l) * +b(l)*s1(l-1) if(l.eq.k) s(k)=s2(k) c c Check impending underflow or overflow c if(abs(s(k)).lt.tiny) then ierr=-(k-1) return else if(abs(s(k)).gt.huge) then ierr=k-1 return end if end do c c Compute the alpha- and beta-coefficient c alpha(k)=a(k)+(s2(k+1)/s2(k))-(s1(k)/s1(k-1)) beta(k)=s2(k)/s1(k-1) do l=k,lk s0(l)=s1(l) s1(l)=s2(l) end do end do return end subroutine chri ( n, iopt, a, b, x, y, hr, hi, alpha, beta, ierr ) c*********************************************************************72 c cc CHRI_R4 implements the Christoffel or generalized Christoffel theorem. c c c This subroutine implements the Christoffel or generalized Christoffel c theorem. In all cases except iopt=7, it uses nonlinear recurrence c algorithms described in W. Gautschi,An algorithmic implementation c of the generalized Christoffel theorem'', Numerical Integration c (G. Haemmerlin, ed.), Birkhaeuser, Basel, 1982, pp. 89-106. The case c iopt=7 incorporates a QR step with shift x in the manner of c J. Kautsky and G.H. Golub, On the calculation of Jacobi matrices'', c Linear Algebra Appl. 52/53, 1983, 439-455, using the algorithm of c Eq. (67.11) on p. 567 in J.H. Wilkinson,The Algebraic Eigenvalue c Problem'', Clarendon Press, Oxford, 1965. Given the recursion c coefficients a(k),b(k), k=0,1,...,n, for the (monic) orthogonal c polynomials with respect to some measure dlambda(t), it generates c the recursion coefficients alpha(k),beta(k), k=0,1,...,n-1, for the c measure c c (t-x)dlambda(t) if iopt=1 c [(t-x)**2+y**2]dlambda(t) if iopt=2 c (t**2+y**2)dlambda(t) with if iopt=3 c dlambda(t) and supp(dlambda) c symmetric with respect to c the origin c dlambda(t)/(t-x) if iopt=4 c dlambda(t)/[(t-x)**2+y**2] if iopt=5 c dlambda(t)/(t**2+y**2) with if iopt=6 c dlambda(t) and supp(dlambda) c symmetric with respect to c the origin c [(t-x)**2]dlambda(t) if iopt=7 c c c Input: n - - - the number of recurrence coefficients c desired; type integer c iopt - - an integer selecting the desired weight c distribution c a,b - - arrays of dimension n+1 containing the c recursion coefficients a(k-1),b(k-1),k=1,2, c ...,n+1, of the polynomials orthogonal with c respect to the given measure dlambda(t) c x,y - - real parameters defining the linear and c quadratic factors, or divisors, of dlambda(t) c hr,hi - the real and imaginary part, respectively, of c the integral of dlambda(t)/(z-t), where z=x+iy; c the parameter hr is used only if iopt=4 or c 5, the parameter hi only if iopt=5 or 6 c c Output: alpha,beta - - arrays of dimension n containing the c desired recursion coefficients alpha(k-1), c beta(k-1), k=1,2,...,n c c It is assumed that n is larger than or equal to 2. Otherwise, the c routine exits immediately with the error flag ierr set equal to 1. c If iopt is not between 1 and 7, the routine exits with ierr=2. c c The routine uses the function subroutine r1mach to evaluate the c constant eps, which is used only if iopt=7. c dimension a(*),b(*),alpha(n),beta(n) c c The arrays a,b are assumed to have dimension n+1. c eps=5.*r1mach(3) c c The quantity eps is a constant slightly larger than the machine c precision. c ierr=0 if(n.lt.2) then ierr=1 return end if c c What follows implements Eq. (3.7) of W. Gautschi, op. cit. c if (iopt.eq.1) then e=0. do 10 k=1,n q=a(k)-e-x beta(k)=q*e e=b(k+1)/q alpha(k)=x+q+e 10 continue c c Set the first beta-coefficient as discussed in Section 5.1 of the c companion paper. c beta(1)=b(1)*(a(1)-x) return c c What follows implements Eq. (4.7) of W. Gautschi, op. cit. c else if(iopt.eq.2) then s=x-a(1) t=y eio=0. do 20 k=1,n d=s*s+t*t er=-b(k+1)*s/d ei=b(k+1)*t/d s=x+er-a(k+1) t=y+ei alpha(k)=x+t*er/ei-s*ei/t beta(k)=t*eio*(1.+(er/ei)**2) eio=ei 20 continue c c Set the first beta-coefficient. c beta(1)=b(1)*(b(2)+(a(1)-x)**2+y*y) return c c What follows implements Eq. (4.8) of W. Gautschi, op. cit. c else if(iopt.eq.3) then t=y eio=0. do 30 k=1,n ei=b(k+1)/t t=y+ei alpha(k)=0. beta(k)=t*eio eio=ei 30 continue c c Set the first beta-coefficient. c beta(1)=b(1)*(b(2)+y*y) return c c What follows implements Eqs. (5.1),(5.2) of W. Gautschi, op. cit. c else if(iopt.eq.4) then alpha(1)=x-b(1)/hr beta(1)=-hr q=-b(1)/hr do 40 k=2,n e=a(k-1)-x-q beta(k)=q*e q=b(k)/e alpha(k)=q+e+x 40 continue return c c What follows implements Eq. (5.8) of W. Gautschi, op. cit. c else if(iopt.eq.5) then nm1=n-1 d=hr*hr+hi*hi eroo=a(1)-x+b(1)*hr/d eioo=-b(1)*hi/d-y alpha(1)=x+hr*y/hi beta(1)=-hi/y alpha(2)=x-b(1)*hi*eroo/(d*eioo)+hr*eioo/hi beta(2)=y*eioo*(1.+(hr/hi)**2) if(n.eq.2) return so=b(2)/(eroo**2+eioo**2) ero=a(2)-x-so*eroo eio=so*eioo-y alpha(3)=x+eroo*eio/eioo+so*eioo*ero/eio beta(3)=-b(1)*hi*eio*(1.+(eroo/eioo)**2)/d if(n.eq.3) return do 50 k=3,nm1 s=b(k)/(ero**2+eio**2) er=a(k)-x-s*ero ei=s*eio-y alpha(k+1)=x+ero*ei/eio+s*eio*er/ei beta(k+1)=so*eioo*ei*(1.+(ero/eio)**2) eroo=ero eioo=eio ero=er eio=ei so=s 50 continue return c c What follows implements Eq. (5.9) of W. Gautschi, op. cit. c else if(iopt.eq.6) then nm1=n-1 eoo=-b(1)/hi-y eo=b(2)/eoo-y alpha(1)=0. beta(1)=-hi/y alpha(2)=0. beta(2)=y*eoo if(n.eq.2) return alpha(3)=0. beta(3)=-b(1)*eo/hi if(n.eq.3) return do 60 k=3,nm1 e=b(k)/eo-y beta(k+1)=b(k-1)*e/eoo alpha(k+1)=0. eoo=eo eo=e 60 continue return c c What follows implements a QR step with shift x. c else if(iopt.eq.7) then u=0. c=1. c0=0. do 70 k=1,n gamma=a(k)-x-u cm1=c0 c0=c if(abs(c0).gt.eps) then p2=(gamma**2)/c0 else p2=cm1*b(k) end if if(k.gt.1) beta(k)=s*(p2+b(k+1)) s=b(k+1)/(p2+b(k+1)) c=p2/(p2+b(k+1)) u=s*(gamma+a(k+1)-x) alpha(k)=gamma+u+x 70 continue beta(1)=b(1)*(b(2)+(x-a(1))**2) return else ierr=2 return end if end subroutine cheb_r8 ( n, da, db, dnu, dalpha, dbeta, ds, iderr, & ds0, ds1, ds2 ) c*********************************************************************72 c cc CHEB_R8 generates recursion coefficients ALPHA and BETA. c c Reference: c c Walter Gautschi, c On Generating Orthogonal Polynomials, c SIAM Journal on Scientific and Statistical Computing, c Volume 3, Number 3, 1982, pages 289-317. c double precision da,db,dnu,dalpha,dbeta,ds,ds0,ds1,ds2,dtiny, *d1mach,dhuge dimension da(*),db(*),dnu(*),dalpha(n),dbeta(n),ds(n), *ds0(*),ds1(*),ds2(*) c c The arrays da,db are assumed to have dimension 2*n-1, the arrays c dnu,ds0,ds1,ds2 dimension 2*n. c nd=2*n dtiny=10.d0*d1mach(1) dhuge=.1d0*d1mach(2) iderr=0 if(dabs(dnu(1)).lt.dtiny) then iderr=1 return end if if(n.lt.1) then iderr=2 return end if dalpha(1)=da(1)+dnu(2)/dnu(1) dbeta(1)=dnu(1) if(n.eq.1) return ds(1)=dnu(1) do 10 l=1,nd ds0(l)=0.d0 ds1(l)=dnu(l) 10 continue do 40 k=2,n lk=nd-k+1 do 20 l=k,lk ds2(l)=ds1(l+1)-(dalpha(k-1)-da(l))*ds1(l)-dbeta(k-1)*ds0(l) * +db(l)*ds1(l-1) if(l.eq.k) ds(k)=ds2(k) 20 continue if(dabs(ds(k)).lt.dtiny) then iderr=-(k-1) return else if(dabs(ds(k)).gt.dhuge) then iderr=k-1 return end if dalpha(k)=da(k)+(ds2(k+1)/ds2(k))-(ds1(k)/ds1(k-1)) dbeta(k)=ds2(k)/ds1(k-1) do 30 l=k,lk ds0(l)=ds1(l) ds1(l)=ds2(l) 30 continue 40 continue return end subroutine chri_r8(n,iopt,da,db,dx,dy,dhr,dhi,dalpha,dbeta,ierr) c*********************************************************************72 c cc CHRI_R8 implements the Christoffel or generalized Christoffel theorem. c double precision da,db,dx,dy,dhr,dhi,dalpha,dbeta,deps,d1mach, *de,dq,ds,dt,deio,dd,der,dei,deroo,deioo,dso,dero,deoo,deo,du,dc, *dc0,dgam,dcm1,dp2 dimension da(*),db(*),dalpha(n),dbeta(n) c c The arrays da,db are assumed to have dimension n+1. c deps=5.d0*d1mach(3) ierr=0 if(n.lt.2) then ierr=1 return end if if(iopt.eq.1) then de=0.d0 do 10 k=1,n dq=da(k)-de-dx dbeta(k)=dq*de de=db(k+1)/dq dalpha(k)=dx+dq+de 10 continue dbeta(1)=db(1)*(da(1)-dx) return else if(iopt.eq.2) then ds=dx-da(1) dt=dy deio=0.d0 do 20 k=1,n dd=ds*ds+dt*dt der=-db(k+1)*ds/dd dei=db(k+1)*dt/dd ds=dx+der-da(k+1) dt=dy+dei dalpha(k)=dx+dt*der/dei-ds*dei/dt dbeta(k)=dt*deio*(1.d0+(der/dei)**2) deio=dei 20 continue dbeta(1)=db(1)*(db(2)+(da(1)-dx)**2+dy*dy) return else if(iopt.eq.3) then dt=dy deio=0.d0 do 30 k=1,n dei=db(k+1)/dt dt=dy+dei dalpha(k)=0.d0 dbeta(k)=dt*deio deio=dei 30 continue dbeta(1)=db(1)*(db(2)+dy*dy) return else if(iopt.eq.4) then dalpha(1)=dx-db(1)/dhr dbeta(1)=-dhr dq=-db(1)/dhr do 40 k=2,n de=da(k-1)-dx-dq dbeta(k)=dq*de dq=db(k)/de dalpha(k)=dq+de+dx 40 continue return else if(iopt.eq.5) then nm1=n-1 dd=dhr*dhr+dhi*dhi deroo=da(1)-dx+db(1)*dhr/dd deioo=-db(1)*dhi/dd-dy dalpha(1)=dx+dhr*dy/dhi dbeta(1)=-dhi/dy dalpha(2)=dx-db(1)*dhi*deroo/(dd*deioo)+dhr*deioo/dhi dbeta(2)=dy*deioo*(1.d0+(dhr/dhi)**2) if(n.eq.2) return dso=db(2)/(deroo**2+deioo**2) dero=da(2)-dx-dso*deroo deio=dso*deioo-dy dalpha(3)=dx+deroo*deio/deioo+dso*deioo*dero/deio dbeta(3)=-db(1)*dhi*deio*(1.d0+(deroo/deioo)**2)/dd if(n.eq.3) return do 50 k=3,nm1 ds=db(k)/(dero**2+deio**2) der=da(k)-dx-ds*dero dei=ds*deio-dy dalpha(k+1)=dx+dero*dei/deio+ds*deio*der/dei dbeta(k+1)=dso*deioo*dei*(1.d0+(dero/deio)**2) deroo=dero deioo=deio dero=der deio=dei dso=ds 50 continue return else if(iopt.eq.6) then nm1=n-1 deoo=-db(1)/dhi-dy deo=db(2)/deoo-dy dalpha(1)=0.d0 dbeta(1)=-dhi/dy dalpha(2)=0.d0 dbeta(2)=dy*deoo if(n.eq.2) return dalpha(3)=0.d0 dbeta(3)=-db(1)*deo/dhi if(n.eq.3) return do 60 k=3,nm1 de=db(k)/deo-dy dbeta(k+1)=db(k-1)*de/deoo dalpha(k+1)=0.d0 deoo=deo deo=de 60 continue return else if(iopt.eq.7) then du=0.d0 dc=1.d0 dc0=0.d0 do 70 k=1,n dgam=da(k)-dx-du dcm1=dc0 dc0=dc if(dabs(dc0).gt.deps) then dp2=(dgam**2)/dc0 else dp2=dcm1*db(k) end if if(k.gt.1) dbeta(k)=ds*(dp2+db(k+1)) ds=db(k+1)/(dp2+db(k+1)) dc=dp2/(dp2+db(k+1)) du=ds*(dgam+da(k+1)-dx) dalpha(k)=dgam+du+dx 70 continue dbeta(1)=db(1)*(db(2)+(dx-da(1))**2) return else ierr=2 return end if end subroutine dgauss(n,dalpha,dbeta,deps,dzero,dweigh,ierr,de) c*********************************************************************72 c cc DGAUSS is a double-precision version of the routine gauss. c double precision dalpha,dbeta,deps,dzero,dweigh,de,dp,dg,dr, *ds,dc,df,db dimension dalpha(n),dbeta(n),dzero(n),dweigh(n),de(n) if(n.lt.1) then ierr=-1 return end if ierr=0 dzero(1)=dalpha(1) if(dbeta(1).lt.0.d0) then ierr=-2 return end if dweigh(1)=dbeta(1) if (n.eq.1) return dweigh(1)=1.d0 de(n)=0.d0 do 100 k=2,n dzero(k)=dalpha(k) if(dbeta(k).lt.0.d0) then ierr=-2 return end if de(k-1)=dsqrt(dbeta(k)) dweigh(k)=0.d0 100 continue do 240 l=1,n j=0 105 do 110 m=l,n if(m.eq.n) goto 120 if(dabs(de(m)).le.deps*(dabs(dzero(m))+dabs(dzero(m+1)))) * goto 120 110 continue 120 dp=dzero(l) if(m.eq.l) goto 240 if(j.eq.30) goto 400 j=j+1 dg=(dzero(l+1)-dp)/(2.d0*de(l)) dr=dsqrt(dg*dg+1.d0) dg=dzero(m)-dp+de(l)/(dg+dsign(dr,dg)) ds=1.d0 dc=1.d0 dp=0.d0 mml=m-l do 200 ii=1,mml i=m-ii df=ds*de(i) db=dc*de(i) if(dabs(df).lt.dabs(dg)) goto 150 dc=dg/df dr=dsqrt(dc*dc+1.d0) de(i+1)=df*dr ds=1.d0/dr dc=dc*ds goto 160 150 ds=df/dg dr=dsqrt(ds*ds+1.d0) de(i+1)=dg*dr dc=1.d0/dr ds=ds*dc 160 dg=dzero(i+1)-dp dr=(dzero(i)-dg)*ds+2.d0*dc*db dp=ds*dr dzero(i+1)=dg+dp dg=dc*dr-db df=dweigh(i+1) dweigh(i+1)=ds*dweigh(i)+dc*df dweigh(i)=dc*dweigh(i)-ds*df 200 continue dzero(l)=dzero(l)-dp de(l)=dg de(m)=0.d0 goto 105 240 continue do 300 ii=2,n i=ii-1 k=i dp=dzero(i) do 260 j=ii,n if(dzero(j).ge.dp) goto 260 k=j dp=dzero(j) 260 continue if(k.eq.i) goto 300 dzero(k)=dzero(i) dzero(i)=dp dp=dweigh(i) dweigh(i)=dweigh(k) dweigh(k)=dp 300 continue do 310 k=1,n dweigh(k)=dbeta(1)*dweigh(k)*dweigh(k) 310 continue return 400 ierr=l return end subroutine dgchri(n,iopt,nu0,numax,deps,da,db,dx,dy,dalpha,dbeta, *nu,ierr,ierrc,dnu,drhor,drhoi,droldr,droldi,ds,ds0,ds1,ds2) c*********************************************************************72 c cc DGCHRI is a double-precision version of the routine gchri. c double precision deps,da(numax),db(numax),dx,dy,dalpha(n), *dbeta(n),dnu(*),drhor(*),drhoi(*),droldr(*),droldi(*), *ds(n),ds0(*),ds1(*),ds2(*) c c The arrays dnu,drhor,drhoi,droldr,droldi,ds0,ds1,ds2 are assumed c to have dimension 2*n. c if(n.lt.1) then ierr=-1 return end if ierr=0 nd=2*n ndm1=nd-1 if(iopt.eq.1) then call dknum(ndm1,nu0,numax,dx,dy,deps,da,db,drhor,drhoi,nu, * ierr,droldr,droldi) do 10 k=1,nd dnu(k)=-drhor(k) 10 continue call cheb_r8 ( n,da,db,dnu,dalpha,dbeta,ds,ierrc,ds0,ds1,ds2) return else if(iopt.eq.2) then dy=dabs(dy) call dknum(ndm1,nu0,numax,dx,dy,deps,da,db,drhor,drhoi,nu, * ierr,droldr,droldi) do 20 k=1,nd dnu(k)=-drhoi(k)/dy 20 continue call cheb_r8 (n,da,db,dnu,dalpha,dbeta,ds,ierrc,ds0,ds1,ds2) return else ierr=1 return end if end subroutine dkern(n,nu0,numax,dx,dy,deps,da,db,dkerr,dkeri, * nu,ierr,droldr,droldi) c*********************************************************************72 c cc DKERN is a double-precision version of the routine kern. c double precision dx,dy,deps,da(numax),db(numax),dkerr(*), * dkeri(*),droldr(*),droldi(*),dp0r,dp0i,dpr,dpi,dpm1r, * dpm1i,dden,dt c c The arrays dkerr,dkeri,droldr,droldi are assumed to have c dimension n+1. c call dknum(n,nu0,numax,dx,dy,deps,da,db,dkerr,dkeri,nu,ierr, * droldr,droldi) if(ierr.ne.0) return dp0r=0.d0 dp0i=0.d0 dpr=1.d0 dpi=0.d0 do 10 k=1,n dpm1r=dp0r dpm1i=dp0i dp0r=dpr dp0i=dpi dpr=(dx-da(k))*dp0r-dy*dp0i-db(k)*dpm1r dpi=(dx-da(k))*dp0i+dy*dp0r-db(k)*dpm1i dden=dpr**2+dpi**2 dt=(dkerr(k+1)*dpr+dkeri(k+1)*dpi)/dden dkeri(k+1)=(dkeri(k+1)*dpr-dkerr(k+1)*dpi)/dden dkerr(k+1)=dt 10 continue return end subroutine dknum(n,nu0,numax,dx,dy,deps,da,db,drhor,drhoi,nu, *ierr,droldr,droldi) c*********************************************************************72 c cc DKNUM is a double-precision version of the routine knum. c double precision dx,dy,deps,da(numax),db(numax),drhor(*), *drhoi(*),droldr(*),droldi(*),drr,dri,dden,dt c c The arrays drhor,drhoi,droldr,droldi are assumed to have c dimension n+1. c ierr=0 np1=n+1 if(nu0.gt.numax) then ierr=nu0 return end if if(nu0.lt.np1) nu0=np1 nu=nu0-5 do 10 k=1,np1 drhor(k)=0.d0 drhoi(k)=0.d0 10 continue 20 nu=nu+5 if(nu.gt.numax) then ierr=numax goto 60 end if do 30 k=1,np1 droldr(k)=drhor(k) droldi(k)=drhoi(k) 30 continue drr=0.d0 dri=0.d0 do 40 j=1,nu j1=nu-j+1 dden=(dx-da(j1)-drr)**2+(dy-dri)**2 drr=db(j1)*(dx-da(j1)-drr)/dden dri=-db(j1)*(dy-dri)/dden if(j1.le.np1) then drhor(j1)=drr drhoi(j1)=dri end if 40 continue do 50 k=1,np1 C Following statement replaced -- authors remark c c if((drhor(k)-droldr(k))**2+(drhoi(k)-droldi(k))**2.gt. c * deps*(drhor(k)**2+drhoi(k)**2)) goto 20 if((drhor(k)-droldr(k))**2+(drhoi(k)-droldi(k))**2.gt. * (deps**2)*(drhor(k)**2+drhoi(k)**2)) goto 20 50 continue 60 if(n.eq.0) return do 70 k=2,np1 dt=drhor(k)*drhor(k-1)-drhoi(k)*drhoi(k-1) drhoi(k)=drhor(k)*drhoi(k-1)+drhoi(k)*drhor(k-1) drhor(k)=dt 70 continue return end subroutine dlancz(n,ncap,dx,dw,dalpha,dbeta,ierr,dp0,dp1) c*********************************************************************72 c cc DLANCZ is a double-precision version of the routine lancz. c double precision dx(ncap),dw(ncap),dalpha(n),dbeta(n), *dp0(ncap),dp1(ncap),dpi,dgam,dsig,dt,dxlam,drho,dtmp, *dtsig,dtk if(n.le.0 .or. n.gt.ncap) then ierr=1 return else ierr=0 end if do 10 i=1,ncap dp0(i)=dx(i) dp1(i)=0.d0 10 continue dp1(1)=dw(1) do 30 i=1,ncap-1 dpi=dw(i+1) dgam=1.d0 dsig=0.d0 dt=0.d0 dxlam=dx(i+1) do 20 k=1,i+1 drho=dp1(k)+dpi dtmp=dgam*drho dtsig=dsig if(drho.le.0.d0) then dgam=1.d0 dsig=0.d0 else dgam=dp1(k)/drho dsig=dpi/drho end if dtk=dsig*(dp0(k)-dxlam)-dgam*dt dp0(k)=dp0(k)-(dtk-dt) dt=dtk if(dsig.le.0.d0) then dpi=dtsig*dp1(k) else dpi=(dt**2)/dsig end if dtsig=dsig dp1(k)=dtmp 20 continue 30 continue do 40 k=1,n dalpha(k)=dp0(k) dbeta(k)=dp1(k) 40 continue return end subroutine dlob(n,dalpha,dbeta,dleft,dright,dzero,dweigh, *ierr,de,da,db) c*********************************************************************72 c cc DLOB is a double-precision version of the routine lob. c double precision dleft,dright,depsma,dp0l,dp0r,dp1l,dp1r,dpm1l, *dpm1r,ddet,dalpha(*),dbeta(*),dzero(*),dweigh(*),de(*),da(*), *db(*),d1mach c c The arrays dalpha,dbeta,dzero,dweigh,de,da,db are assumed to have c dimension n+2. c depsma=d1mach(3) c c depsma is the machine double precision. c np1=n+1 np2=n+2 do 10 k=1,np2 da(k)=dalpha(k) db(k)=dbeta(k) 10 continue dp0l=0.d0 dp0r=0.d0 dp1l=1.d0 dp1r=1.d0 do 20 k=1,np1 dpm1l=dp0l dp0l=dp1l dpm1r=dp0r dp0r=dp1r dp1l=(dleft-da(k))*dp0l-db(k)*dpm1l dp1r=(dright-da(k))*dp0r-db(k)*dpm1r 20 continue ddet=dp1l*dp0r-dp1r*dp0l da(np2)=(dleft*dp1l*dp0r-dright*dp1r*dp0l)/ddet db(np2)=(dright-dleft)*dp1l*dp1r/ddet call dgauss(np2,da,db,depsma,dzero,dweigh,ierr,de) return end subroutine dmcdis(n,ncapm,mc,mp,dxp,dyp,dquad,deps,iq,idelta, *irout,finld,finrd,dendl,dendr,dxfer,dwfer,dalpha,dbeta,ncap, *kount,ierrd,ied,dbe,dx,dw,dxm,dwm,dp0,dp1,dp2) c*********************************************************************72 c cc DMCDIS is a double-precision version of the routine mcdis. c double precision dxp,dyp,deps,dendl,dendr,dxfer,dwfer,dalpha, *dbeta,dbe,dx,dw,dxm,dwm,dp0,dp1,dp2 dimension dxp(*),dyp(*),dendl(mc),dendr(mc),dxfer(ncapm), *dwfer(ncapm),dalpha(n),dbeta(n),dbe(n),dx(ncapm),dw(ncapm), *dxm(*),dwm(*),dp0(*),dp1(*),dp2(*) logical finld,finrd c c The arrays dxp,dyp are assumed to have dimension mp if mp > 0, c the arrays dxm,dwm,dp0,dp1,dp2 dimension mc*ncapm+mp. c if(idelta.le.0) idelta=1 if(n.lt.1) then ierrd=-1 return end if incap=1 kount=-1 ierr=0 do 10 k=1,n dbeta(k)=0.d0 10 continue ncap=(2*n-1)/idelta 20 do 30 k=1,n dbe(k)=dbeta(k) 30 continue kount=kount+1 if(kount.gt.1) incap=2**(kount/5)*n ncap=ncap+incap if(ncap.gt.ncapm) then ierrd=ncapm return end if mtncap=mc*ncap do 50 i=1,mc im1tn=(i-1)*ncap if(iq.eq.1) then call dquad(ncap,dx,dw,i,ierr) else call dqgp(ncap,dx,dw,i,ierr,mc,finld,finrd,dendl,dendr,dxfer, * dwfer) end if if(ierr.ne.0) then ierrd=i return end if do 40 k=1,ncap dxm(im1tn+k)=dx(k) dwm(im1tn+k)=dw(k) 40 continue 50 continue if(mp.ne.0) then do 60 k=1,mp dxm(mtncap+k)=dxp(k) dwm(mtncap+k)=dyp(k) 60 continue end if if(irout.eq.1) then call dsti(n,mtncap+mp,dxm,dwm,dalpha,dbeta,ied,dp0,dp1,dp2) else call dlancz(n,mtncap+mp,dxm,dwm,dalpha,dbeta,ied,dp0,dp1) end if do 70 k=1,n if(dabs(dbeta(k)-dbe(k)).gt.deps*dabs(dbeta(k))) goto 20 70 continue return end subroutine dmcheb(n,ncapm,mcd,mp,dxp,dyp,dquad,deps,iq, *idelta,finld,finrd,dendl,dendr,dxfer,dwfer,da,db,dnu,dalpha, *dbeta,ncap,kount,ierrd,dbe,dx,dw,dxm,dwm,ds,ds0,ds1,ds2) c*********************************************************************72 c cc DMCHEB is a double-precision version of the routine mccheb. c double precision dxp,dyp,deps,dendl,dendr,dxfer,dwfer,da,db, *dnu,dalpha,dbeta,dbe,dx,dw,dxm,dwm,ds,ds0,ds1,ds2,dsum,dp1, *dp,dpm1 dimension dxp(*),dyp(*),dendl(mcd),dendr(mcd),dxfer(ncapm), *dwfer(ncapm),da(*),db(*),dnu(*),dalpha(n),dbeta(n),dbe(n), *dx(ncapm),dw(ncapm),dxm(*),dwm(*),ds(n),ds0(*),ds1(*),ds2(*) logical finld,finrd c c The arrays dxp,dyp are assumed to have dimension mp if mp > 0, c the arrays da,db dimension 2*n-1, the arrays dnu,ds0,ds1,ds2 c dimension 2*n, and the arrays dxm,dwm dimension mc*ncapm+mp. c nd=2*n if(idelta.le.0) idelta=1 if(n.lt.1) then ierrd=-1 return end if incap=1 kount=-1 ierrd=0 do 10 k=1,n dbeta(k)=0.d0 10 continue ncap=(nd-1)/idelta 20 do 30 k=1,n dbe(k)=dbeta(k) 30 continue kount=kount+1 if(kount.gt.1) incap=2**(kount/5)*n ncap=ncap+incap if(ncap.gt.ncapm) then ierrd=ncapm return end if mtncap=mcd*ncap do 50 i=1,mcd im1tn=(i-1)*ncap if(iq.eq.1) then call dquad(ncap,dx,dw,i,ierr) else call dqgp(ncap,dx,dw,i,ierr,mcd,finld,finrd,dendl,dendr, * dxfer,dwfer) end if if(ierr.ne.0) then ierrd=i return end if do 40 k=1,ncap dxm(im1tn+k)=dx(k) dwm(im1tn+k)=dw(k) 40 continue 50 continue if(mp.ne.0) then do 60 k=1,mp dxm(mtncap+k)=dxp(k) dwm(mtncap+k)=dyp(k) 60 continue end if mtnpmp=mtncap+mp do 90 k=1,nd km1=k-1 dsum=0.d0 do 80 i=1,mtnpmp dp1=0.d0 dp=1.d0 if(k.gt.1) then do 70 l=1,km1 dpm1=dp1 dp1=dp dp=(dxm(i)-da(l))*dp1-db(l)*dpm1 70 continue end if dsum=dsum+dwm(i)*dp 80 continue dnu(k)=dsum 90 continue call cheb_r8 (n,da,db,dnu,dalpha,dbeta,ds,ierr,ds0,ds1,ds2) do 100 k=1,n if(dabs(dbeta(k)-dbe(k)).gt.deps*dabs(dbeta(k))) goto 20 100 continue return end subroutine dqgp(n,dx,dw,i,ierr,mcd,finld,finrd,dendl,dendr, * dxfer,dwfer) c*********************************************************************72 c cc DQGP is a double-precision version of the routine qgp. c c The user has to supply the routine c c double precision function dwf(dx,i), c c which evaluates the weight function in double precision at the c point dx on the i-th component interval. c double precision dx,dw,dendl,dendr,dxfer,dwfer,dphi,dphi1,dwf dimension dx(n),dw(n),dendl(mcd),dendr(mcd),dxfer(*),dwfer(*) logical finld,finrd c c The arrays dxfer,dwfer are dimensioned in the routine dmcdis. c ierr=0 if(i.eq.1) call dfejer(n,dxfer,dwfer) if(i.gt.1 .and. i.lt.mcd) goto 60 if(mcd.eq.1) then if(finld.and.finrd) goto 60 if(finld) goto 20 if(finrd) goto 40 do 10 k=1,n call dsymtr(dxfer(k),dphi,dphi1) dx(k)=dphi dw(k)=dwfer(k)*dwf(dphi,i)*dphi1 10 continue return else if((i.eq.1.and.finld).or.(i.eq.mcd.and.finrd)) goto 60 if(i.eq.1) goto 40 end if 20 do 30 k=1,n call dtr(dxfer(k),dphi,dphi1) dx(k)=dendl(mcd)+dphi dw(k)=dwfer(k)*dwf(dx(k),mcd)*dphi1 30 continue return 40 do 50 k=1,n call dtr(-dxfer(k),dphi,dphi1) dx(k)=dendr(1)-dphi dw(k)=dwfer(k)*dwf(dx(k),1)*dphi1 50 continue return 60 do 70 k=1,n dx(k)=.5d0*((dendr(i)-dendl(i))*dxfer(k)+dendr(i)+dendl(i)) dw(k)=.5d0*(dendr(i)-dendl(i))*dwfer(k)*dwf(dx(k),i) 70 continue return end subroutine dsymtr(dt,dphi,dphi1) c*********************************************************************72 c cc DSYMTR is a double-precision version of symtr. c double precision dt,dphi,dphi1,dt2 dt2=dt*dt dphi=dt/(1.-dt2) dphi1=(dt2+1.d0)/(dt2-1.d0)**2 return end subroutine dtr(dt,dphi,dphi1) c*********************************************************************72 c cc DTR is a double-precision version of tr. c double precision dt,dphi,dphi1 dphi=(1.d0+dt)/(1.d0-dt) dphi1=2.d0/(dt-1.d0)**2 return end subroutine dfejer(n,dx,dw) c*********************************************************************72 c cc DFEJER is a double-precision version of fejer. c double precision dx,dw,dpi,dn,dc1,dc0,dt,dsum,dc2 dimension dx(n),dw(n) dpi=4.d0*datan(1.d0) nh=n/2 np1h=(n+1)/2 dn=dble(n) do 10 k=1,nh dx(n+1-k)=dcos(.5d0*dble(2*k-1)*dpi/dn) dx(k)=-dx(n+1-k) 10 continue if(2*nh.ne.n) dx(np1h)=0.d0 do 30 k=1,np1h dc1=1.d0 dc0=2.d0*dx(k)*dx(k)-1.d0 dt=2.d0*dc0 dsum=dc0/3.d0 do 20 m=2,nh dc2=dc1 dc1=dc0 dc0=dt*dc1-dc2 dsum=dsum+dc0/dble(4*m*m-1) 20 continue dw(k)=2.d0*(1.d0-2.d0*dsum)/dn dw(n+1-k)=dw(k) 30 continue return end subroutine dradau(n,dalpha,dbeta,dend,dzero,dweigh,ierr,de, *da,db) c*********************************************************************72 c c This is a double-precision version of the routine radau. c double precision dend,depsma,dp0,dp1,dpm1,dalpha(*),dbeta(*), *dzero(*),dweigh(*),de(*),da(*),db(*),d1mach c c The arrays dalpha,dbeta,dzero,dweigh,de,da,db are assumed to have c dimension n+1. c depsma=d1mach(3) c c depsma is the machine double precision. c np1=n+1 do 10 k=1,np1 da(k)=dalpha(k) db(k)=dbeta(k) 10 continue dp0=0.d0 dp1=1.d0 do 20 k=1,n dpm1=dp0 dp0=dp1 dp1=(dend-da(k))*dp0-db(k)*dpm1 20 continue da(np1)=dend-db(np1)*dp0/dp1 call dgauss(np1,da,db,depsma,dzero,dweigh,ierr,de) return end subroutine drecur(n,ipoly,dal,dbe,da,db,iderr) c*********************************************************************72 c c This is a double-precision version of the routine recur. c external dgamma double precision dal,dbe,da,db,dlmach,d1mach,dkm1,dalpbe,dt, *dlga,dal2,dbe2,dgamma dimension da(n),db(n) if(n.lt.1) then iderr=3 return end if dlmach=dlog(d1mach(2)) iderr=0 do 10 k=1,n da(k)=0.d0 10 continue if(ipoly.eq.1) then db(1)=2.d0 if (n.eq.1) return do 20 k=2,n dkm1=dble(k-1) db(k)=1.d0/(4.d0-1.d0/(dkm1*dkm1)) 20 continue return else if(ipoly.eq.2) then da(1)=.5d0 db(1)=1.d0 if(n.eq.1) return do 30 k=2,n da(k)=.5d0 dkm1=dble(k-1) db(k)=.25d0/(4.d0-1.d0/(dkm1*dkm1)) 30 continue return else if(ipoly.eq.3) then db(1)=4.d0*datan(1.d0) if(n.eq.1) return db(2)=.5d0 if(n.eq.2) return do 40 k=3,n db(k)=.25d0 40 continue return else if(ipoly.eq.4) then db(1)=2.d0*datan(1.d0) if(n.eq.1) return do 50 k=2,n db(k)=.25d0 50 continue return else if(ipoly.eq.5) then db(1)=4.d0*datan(1.d0) da(1)=.5d0 if(n.eq.1) return do 60 k=2,n db(k)=.25d0 60 continue return else if(ipoly.eq.6) then if(dal.le.-1.d0 .or. dbe.le.-1.d0) then iderr=1 return else dalpbe=dal+dbe da(1)=(dbe-dal)/(dalpbe+2.d0) dt=(dalpbe+1.d0)*dlog(2.d0)+dlga(dal+1.d0)+dlga(dbe+1.d0)- * dlga(dalpbe+2.d0) if(dt.gt.dlmach) then iderr=2 db(1)=d1mach(2) else db(1)=dexp(dt) end if if(n.eq.1) return dal2=dal*dal dbe2=dbe*dbe da(2)=(dbe2-dal2)/((dalpbe+2.d0)*(dalpbe+4.d0)) db(2)=4.d0*(dal+1.d0)*(dbe+1.d0)/((dalpbe+3.d0)*(dalpbe+ * 2.d0)**2) if(n.eq.2) return do 70 k=3,n dkm1=dble(k-1) da(k)=.25d0*(dbe2-dal2)/(dkm1*dkm1*(1.d0+.5d0*dalpbe/dkm1) * *(1.d0+.5d0*(dalpbe+2.d0)/dkm1)) db(k)=.25d0*(1.d0+dal/dkm1)*(1.d0+dbe/dkm1)*(1.d0+dalpbe/ * dkm1)/((1.d0+.5d0*(dalpbe+1.d0)/dkm1)*(1.d0+.5d0*(dalpbe * -1.d0)/dkm1)*(1.d0+.5d0*dalpbe/dkm1)**2) 70 continue return end if else if(ipoly.eq.7) then if(dal.le.-1.d0) then iderr=1 return else da(1)=dal+1.d0 db(1)=dgamma(dal+1.d0,iderr) if(iderr.eq.2) db(1)=d1mach(2) if(n.eq.1) return do 80 k=2,n dkm1=dble(k-1) da(k)=2.d0*dkm1+dal+1.d0 db(k)=dkm1*(dkm1+dal) 80 continue return end if else if(ipoly.eq.8) then db(1)=dsqrt(4.d0*datan(1.d0)) if(n.eq.1) return do 90 k=2,n db(k)=.5d0*dble(k-1) 90 continue return else iderr=4 end if return end double precision function dlga(dx) c*********************************************************************72 c c This routine evaluates the logarithm of the gamma function by a c combination of recurrence and asymptotic approximation. c c Reference: c c William Cody, Kenneth Hillstrom, c Chebyshev Approximations for the Natural Logarithm of the c Gamma Function, c Mathematics of Computation, c Volume 21, Number 98, April 1967, pages 198-203. c c c The entries in the next data statement are the numerators and c denominators, respectively, of the quantities B[16]/(16*15), c B[14]/(14*13),..., B[2]/(2*1), where B[2n] are the Bernoulli c numbers. c double precision dbnum,dbden,dx,d1mach,dc,dp,dy,dt,ds dimension dbnum(8),dbden(8) data dbnum/-3.617d3,1.d0,-6.91d2,1.d0,-1.d0,1.d0,-1.d0,1.d0/, * dbden/1.224d5,1.56d2,3.6036d5,1.188d3,1.68d3,1.26d3,3.6d2, *1.2d1/ c c The quantity dprec in the next statement is the number of decimal c digits carried in double-precision floating-point arithmetic. c dprec=-alog10(sngl(d1mach(3))) dc=.5d0*dlog(8.d0*datan(1.d0)) dp=1.d0 dy=dx y=sngl(dy) c c The quantity y0 below is the threshold value beyond which asymptotic c evaluation gives sufficient accuracy; see Eq. 6.1.42 in M. Abramowitz c and I.A. Stegun,Handbook of Mathematical Functions''. The constants c are .12118868... = ln(10)/19 and .05390522... = ln(|B[20]|/190)/19. c y0=exp(.121189*dprec+.053905) 10 if(y.gt.y0) goto 20 dp=dy*dp dy=dy+1.d0 y=sngl(dy) goto 10 20 dt=1.d0/(dy*dy) c c The right-hand side of the next assignment statement is B[18]/(18*17). c ds=4.3867d4/2.44188d5 do 30 i=1,8 ds=dt*ds+dbnum(i)/dbden(i) 30 continue dlga=(dy-.5d0)*dlog(dy)-dy+dc+ds/dy-dlog(dp) return end double precision function dgamma(dx,iderr) c*********************************************************************72 c c This evaluates the gamma function for real positive dx, using the c function subroutine dlga. c double precision dx,dlmach,d1mach,dt,dlga dlmach=dlog(d1mach(2)) iderr=0 dt=dlga(dx) if(dt.ge.dlmach) then iderr=2 dgamma=d1mach(2) return else dgamma=dexp(dt) return end if end subroutine dsti(n,ncap,dx,dw,dalpha,dbeta,ierr,dp0,dp1,dp2) c*********************************************************************72 c c This is a double-precision version of the routine sti. c c Reference: c c Walter Gautschi, c On Generating Orthogonal Polynomials, c SIAM Journal on Scientific and Statistical Computing, c Volume 3, Number 3, 1982, pages 289-317. c double precision dx,dw,dalpha,dbeta,dp0,dp1,dp2,dtiny,d1mach, *dhuge,dsum0,dsum1,dsum2,dt dimension dx(ncap),dw(ncap),dalpha(n),dbeta(n),dp0(ncap), *dp1(ncap),dp2(ncap) dtiny=10.d0*d1mach(1) dhuge=.1d0*d1mach(2) ierr=0 if(n.le.0 .or. n.gt.ncap) then ierr=1 return end if nm1=n-1 dsum0=0.d0 dsum1=0.d0 do 10 m=1,ncap dsum0=dsum0+dw(m) dsum1=dsum1+dw(m)*dx(m) 10 continue dalpha(1)=dsum1/dsum0 dbeta(1)=dsum0 if(n.eq.1) return do 20 m=1,ncap dp1(m)=0.d0 dp2(m)=1.d0 20 continue do 40 k=1,nm1 dsum1=0.d0 dsum2=0.d0 do 30 m=1,ncap if(dw(m).eq.0.d0) goto 30 dp0(m)=dp1(m) dp1(m)=dp2(m) dp2(m)=(dx(m)-dalpha(k))*dp1(m)-dbeta(k)*dp0(m) if(dabs(dp2(m)).gt.dhuge .or. dabs(dsum2).gt.dhuge) then ierr=k return end if dt=dw(m)*dp2(m)*dp2(m) dsum1=dsum1+dt dsum2=dsum2+dt*dx(m) 30 continue if(dabs(dsum1).lt.dtiny) then ierr=-k return end if dalpha(k+1)=dsum2/dsum1 dbeta(k+1)=dsum1/dsum0 dsum0=dsum1 40 continue return end subroutine gauss(n,alpha,beta,eps,zero,weight,ierr,e) c*********************************************************************72 c c Given n and a measure dlambda, this routine generates the n-point c Gaussian quadrature formula c c integral over supp(dlambda) of f(x)dlambda(x) c c = sum from k=1 to k=n of w(k)f(x(k)) + R(n;f). c c The nodes are returned as zero(k)=x(k) and the weights as c weight(k)=w(k), k=1,2,...,n. The user has to supply the recursion c coefficients alpha(k), beta(k), k=0,1,2,...,n-1, for the measure c dlambda. The routine computes the nodes as eigenvalues, and the c weights in term of the first component of the respective normalized c eigenvectors of the n-th order Jacobi matrix associated with dlambda. c It uses a translation and adaptation of the algol procedure imtql2, c Numer. Math. 12, 1968, 377-383, by Martin and Wilkinson, as modified c by Dubrulle, Numer. Math. 15, 1970, 450. See also Handbook for c Autom. Comput., vol. 2 - Linear Algebra, pp.241-248, and the eispack c routine imtql2. c c Input: n - - the number of points in the Gaussian quadrature c formula; type integer c alpha,beta - - arrays of dimension n to be filled c with the values of alpha(k-1), beta(k-1), k=1,2, c ...,n c eps - the relative accuracy desired in the nodes c and weights c c Output: zero- array of dimension n containing the Gaussian c nodes (in increasing order) zero(k)=x(k), k=1,2, c ...,n c weight - array of dimension n containing the c Gaussian weights weight(k)=w(k), k=1,2,...,n c ierr- an error flag equal to 0 on normal return, c equal to i if the QR algorithm does not c converge within 30 iterations on evaluating the c i-th eigenvalue, equal to -1 if n is not in c range, and equal to -2 if one of the beta's is c negative. c c The array e is needed for working space. c dimension alpha(n),beta(n),zero(n),weight(n),e(n) if(n.lt.1) then ierr=-1 return end if ierr=0 zero(1)=alpha(1) if(beta(1).lt.0.) then ierr=-2 return end if weight(1)=beta(1) if (n.eq.1) return weight(1)=1. e(n)=0. do 100 k=2,n zero(k)=alpha(k) if(beta(k).lt.0.) then ierr=-2 return end if e(k-1)=sqrt(beta(k)) weight(k)=0. 100 continue do 240 l=1,n j=0 c c Look for a small subdiagonal element. c 105 do 110 m=l,n if(m.eq.n) goto 120 if(abs(e(m)).le.eps*(abs(zero(m))+abs(zero(m+1)))) goto 120 110 continue 120 p=zero(l) if(m.eq.l) goto 240 if(j.eq.30) goto 400 j=j+1 c c Form shift. c g=(zero(l+1)-p)/(2.*e(l)) r=sqrt(g*g+1.) g=zero(m)-p+e(l)/(g+sign(r,g)) s=1. c=1. p=0. mml=m-l c c For i=m-1 step -1 until l do ... c do 200 ii=1,mml i=m-ii f=s*e(i) b=c*e(i) if(abs(f).lt.abs(g)) goto 150 c=g/f r=sqrt(c*c+1.) e(i+1)=f*r s=1./r c=c*s goto 160 150 s=f/g r=sqrt(s*s+1.) e(i+1)=g*r c=1./r s=s*c 160 g=zero(i+1)-p r=(zero(i)-g)*s +2.*c*b p=s*r zero(i+1)=g+p g=c*r-b c c Form first component of vector. c f=weight(i+1) weight(i+1)=s*weight(i)+c*f weight(i)=c*weight(i)-s*f 200 continue zero(l)=zero(l)-p e(l)=g e(m)=0. goto 105 240 continue c c Order eigenvalues and eigenvectors. c do 300 ii=2,n i=ii-1 k=i p=zero(i) do 260 j=ii,n if(zero(j).ge.p) goto 260 k=j p=zero(j) 260 continue if(k.eq.i) goto 300 zero(k)=zero(i) zero(i)=p p=weight(i) weight(i)=weight(k) weight(k)=p 300 continue do 310 k=1,n weight(k)=beta(1)*weight(k)*weight(k) 310 continue return c c Set error - no convergence to an eigenvalue after 30 iterations. c 400 ierr=l return end subroutine gchri(n,iopt,nu0,numax,eps,a,b,x,y,alpha,beta,nu, * ierr,ierrc,fnu,rho,rold,s,s0,s1,s2) c*********************************************************************72 c c This routine implements the generalized Christoffel theorem, using c the method of modified moments (cf. Section 4 of W. Gautschi, c Minimal solutions of three-term recurrence relations and orthogonal c polynomials'', Math. Comp. 36, 1981, 547-554). Given the recursion c coefficients a(k), b(k), k=0,1,...n, for the (monic) orthogonal c polynomials with respect to some measure dlambda(t), it generates c the recursion coefficients alpha(k), beta(k), k=0,1,2,...,n-1 for c the measure c c dlambda(t)/(t-x) if iopt=1 c dlambda(t)/{(t-x)**2+y**2} if iopt=2 c c Input: n - - the number of recurrence coefficients desired; c type integer c iopt - an integer selecting the desired weight distribution c nu0 - an integer estimating the starting backward c recurrence index; in the absence of any better c choice, take nu0 = 3*n c numax - an integer controlling termination of backward c recursion in case of nonconvergence; a conservative c choice is numax = 500 c eps - - a relative error tolerance; type real c a,b - - arrays of dimension numax to be supplied with the c recursion coefficients a(k)=alpha(k-1),b(k)=beta(k), c k=1,2,...,numax, for the measure dlambda c x,y - - real parameters defining the linear and quadratic c divisors of dlambda c c Output: alpha,beta - arrays of dimension n containing the desired c recursion coefficients alpha(k-1), beta(k-1), k=1, c 2,...,n c nu - - the backward recurrence index yielding convergence; c in case of nonconvergence, nu will have the value c numax c ierr - an error flag, where c ierr=0 on normal return c ierr=1 if iopt is neither 1 nor 2 c ierr=nu0 if nu0 > numax c ierr=numax if the backward recurrence algorithm does c not converge c ierr=-1 if n is not in range c ierrc - an error flag inherited from the routine cheb c c The arrays fnu,s,s0,s1,s2 are working space. The routine calls c upon the routines knum and cheb. c complex rho,rold,z dimension a(numax),b(numax),alpha(n),beta(n),fnu(*),rho(*), *rold(*),s(n),s0(*),s1(*),s2(*) c c The arrays fnu,rho,rold,s0,s1,s2 are assumed to have dimension 2*n. c if(n.lt.1) then ierr=-1 return end if ierr=0 nd=2*n ndm1=nd-1 c c Linear divisor c if(iopt.eq.1) then c c Generate the modified moments of dlambda. c z=cmplx(x,0.) call knum(ndm1,nu0,numax,z,eps,a,b,rho,nu,ierr,rold) do 10 k=1,nd fnu(k)=-real(rho(k)) 10 continue c c Compute the desired recursion coefficients by means of the modified c Chebyshev algorithm. c call cheb_r4 (n,a,b,fnu,alpha,beta,s,ierrc,s0,s1,s2) return c c Quadratic divisor c else if(iopt.eq.2) then c c Generate the modified moments of dlambda. c y=abs(y) z=cmplx(x,y) call knum(ndm1,nu0,numax,z,eps,a,b,rho,nu,ierr,rold) do 20 k=1,nd fnu(k)=-aimag(rho(k))/y 20 continue c c Compute the desired recursion coefficients by means of the modified c Chebyshev algorithm. c call cheb_r4 (n,a,b,fnu,alpha,beta,s,ierrc,s0,s1,s2) return else ierr=1 return end if end subroutine kern(n,nu0,numax,z,eps,a,b,ker,nu,ierr,rold) c*********************************************************************72 c c This routine generates the kernels in the Gauss quadrature remainder c term, namely c c K(k)(z)=rho(k)(z)/pi(k)(z), k=0,1,2,...,n, c c where rho(k) are the output quantities of the routine knum, and c pi(k) the (monic) orthogonal polynomials. The results are returned c in the array ker as ker(k)=K(k-1)(z), k=1,2,...,n+1. All the other c input and output parameters have the same meaning as in the routine c knum. c complex z,ker,rold,p0,p,pm1 dimension a(numax),b(numax),ker(*),rold(*) c c The arrays ker,rold are assumed to have dimension n+1. c call knum(n,nu0,numax,z,eps,a,b,ker,nu,ierr,rold) p0=(0.,0.) p=(1.,0.) do 10 k=1,n pm1=p0 p0=p p=(z-a(k))*p0-b(k)*pm1 ker(k+1)=ker(k+1)/p 10 continue return end subroutine knum(n,nu0,numax,z,eps,a,b,rho,nu,ierr,rold) c*********************************************************************72 c c This routine generates c c rho(k)(z)=integral pi(k)(t)dlambda(t)/(z-t), k=0,1,2,...,n, c c where pi(k)(t) is the (monic) k-th degree orthogonal polynomial c with respect to the measure dlambda(t), and the integral is extended c over the support of dlambda. It is assumed that z is a complex c number outside the smallest interval containing the support of c dlambda. The quantities rho(k)(z) are computed as the first n+1 c members of the minimal solution of the basic three-term recurrence c relation c c y(k+1)(z)=(z-a(k))y(k)(z)-b(k)y(k-1)(z), k=0,1,2,..., c c satisfied by the orthogonal polynomials pi(k)(z). c c Input: n - - the largest integer k for which rho(k) is c desired c nu0 - an estimate of the starting backward recurrence c index; if no better estimate is known, set c nu0 = 3*n/2; for Jacobi, Laguerre and Hermite c weight functions, estimates of nu0 are generated c respectively by the routines nu0jac,nu0lag and c nu0her c numax - an integer larger than n cutting off backward c recursion in case of nonconvergence; if nu0 c exceeds numax, then the routine aborts with the c error flag ierr set equal to nu0 c z - - - the variable in rho(k)(z); type complex c eps - - the relative accuracy to which the rho(k) are c desired c a,b - - arrays of dimension numax to be supplied with the c recurrence coefficients a(k-1), b(k-1), k=1,2,..., c numax. c c Output: rho - - an array of dimension n+1 containing the results c rho(k)=rho(k-1)(z), k=1,2,...,n+1; type complex c nu - - the starting backward recurrence index that yields c convergence c ierr - an error flag equal to zero on normal return, equal c to nu0 if nu0 > numax, and equal to numax in c case of nonconvergence. c c The complex array rold of dimension n+1 is used for working space. c complex z,rho,rold,r dimension a(numax),b(numax),rho(*),rold(*) c c The arrays rho,rold are assumed to have dimension n+1. c ierr=0 np1=n+1 if(nu0.gt.numax) then ierr=nu0 return end if if(nu0.lt.np1) nu0=np1 nu=nu0-5 do 10 k=1,np1 rho(k)=(0.,0.) 10 continue 20 nu=nu+5 if(nu.gt.numax) then ierr=numax goto 60 end if do 30 k=1,np1 rold(k)=rho(k) 30 continue r=(0.,0.) do 40 j=1,nu j1=nu-j+1 r=cmplx(b(j1),0.)/(z-cmplx(a(j1),0.)-r) if(j1.le.np1) rho(j1)=r 40 continue do 50 k=1,np1 if(cabs(rho(k)-rold(k)).gt.eps*cabs(rho(k))) goto 20 50 continue 60 if(n.eq.0) return do 70 k=2,np1 rho(k)=rho(k)*rho(k-1) 70 continue return end subroutine lancz(n,ncap,x,w,alpha,beta,ierr,p0,p1) c*********************************************************************72 c c This routine carries out the same task as the routine sti, but c uses the more stable Lanczos method. The meaning of the input c and output parameters is the same as in the routine sti. (This c routine is adapted from the routine RKPW in W.B. Gragg and c W.J. Harrod,The numerically stable reconstruction of Jacobi c matrices from spectral data'', Numer. Math. 44, 1984, 317-335.) c dimension x(ncap),w(ncap),alpha(n),beta(n),p0(ncap),p1(ncap) if(n.le.0 .or. n.gt.ncap) then ierr=1 return else ierr=0 end if do 10 i=1,ncap p0(i)=x(i) p1(i)=0. 10 continue p1(1)=w(1) do 30 i=1,ncap-1 pi=w(i+1) gam=1. sig=0. t=0. xlam=x(i+1) do 20 k=1,i+1 rho=p1(k)+pi tmp=gam*rho tsig=sig if(rho.le.0.) then gam=1. sig=0. else gam=p1(k)/rho sig=pi/rho end if tk=sig*(p0(k)-xlam)-gam*t p0(k)=p0(k)-(tk-t) t=tk if(sig.le.0.) then pi=tsig*p1(k) else pi=(t**2)/sig end if tsig=sig p1(k)=tmp 20 continue 30 continue do 40 k=1,n alpha(k)=p0(k) beta(k)=p1(k) 40 continue return end subroutine lob(n,alpha,beta,aleft,right,zero,weight,ierr,e,a,b) c*********************************************************************72 c c Given n and a measure dlambda, this routine generates the c (n+2)-point Gauss-Lobatto quadrature formula c c integral over supp(dlambda) of f(x)dlambda(x) c c = w(0)f(x(0)) + sum from k=1 to k=n of w(k)f(x(k)) c c + w(n+1)f(x(n+1)) + R(n;f). c c The nodes are returned as zero(k)=x(k), the weights as weight(k) c =w(k), k=0,1,...,n,n+1. The user has to supply the recursion c coefficients alpha(k), beta(k), k=0,1,...,n,n+1, for the measure c dlambda. The nodes and weights are computed in terms of the c eigenvalues and first component of the normalized eigenvectors of c a slightly modified Jacobi matrix of order n+2. The routine calls c upon the subroutine gauss and the function subroutine r1mach. c c Input: n - - the number of interior points in the Gauss-Lobatto c formula; type integer c alpha,beta - arrays of dimension n+2 to be supplied with c the recursion coefficients alpha(k-1), beta(k-1), c k=1,2,...,n+2, of the underlying measure; the c routine does not use alpha(n+2), beta(n+2) c aleft,right - the prescribed left and right endpoints c x(0) and x(n+1) of the Gauss-Lobatto formula c c Output: zero - an array of dimension n+2 containing the nodes (in c increasing order) zero(k)=x(k), k=0,1,...,n,n+1 c weight-an array of dimension n+2 containing the weights c weight(k)=w(k), k=0,1,...,n,n+1 c ierr - an error flag inherited from the routine gauss c c The arrays e,a,b are needed for working space. c dimension alpha(*),beta(*),zero(*),weight(*),e(*),a(*),b(*) c c The arrays alpha,beta,zero,weight,e,a,b are assumed to have c dimension n+2. c epsma=r1mach(3) c c epsma is the machine single precision. c np1=n+1 np2=n+2 do k=1,np2 a(k)=alpha(k) b(k)=beta(k) end do p0l=0.0 p0r=0.0 p1l=1.0 p1r=1.0 do k=1,np1 pm1l=p0l p0l=p1l pm1r=p0r p0r=p1r p1l=(aleft-a(k))*p0l-b(k)*pm1l p1r=(right-a(k))*p0r-b(k)*pm1r end do det=p1l*p0r-p1r*p0l a(np2)=(aleft*p1l*p0r-right*p1r*p0l)/det b(np2)=(right-aleft)*p1l*p1r/det call gauss(np2,a,b,epsma,zero,weight,ierr,e) return end subroutine mccheb(n,ncapm,mc,mp,xp,yp,quad,eps,iq,idelta, *finl,finr,endl,endr,xfer,wfer,a,b,fnu,alpha,beta,ncap,kount, *ierr,be,x,w,xm,wm,s,s0,s1,s2) c*********************************************************************72 c c This is a multiple-component discretized modified Chebyshev c algorithm, basically a modified Chebyshev algorithm in which the c modified moments are discretized in the same manner as the inner c product in the discretization procedure mcdis. The input and c output parameters are as in mcdis. In addition, the arrays a,b c must be filled with the recursion coefficients a(k-1),b(k-1), c k=1,2,...,2*n-1, defining the modified moments. The arrays c be,x,w,xm,wm,s,s0,s1,s2 are used for working space. The routine c calls upon the subroutine cheb. The routine exits immediately with c ierr=-1 if n is not in range. c dimension xp(*),yp(*),endl(mc),endr(mc),xfer(ncapm), *wfer(ncapm),a(*),b(*),fnu(*),alpha(n),beta(n),be(n),x(ncapm), *w(ncapm),xm(*),wm(*),s(n),s0(*),s1(*),s2(*) logical finl,finr c c The arrays xp,yp are assumed to have dimension mp if mp > 0, c the arrays a,b dimension 2*n-1, the arrays fnu,s0,s1,s2 dimension c 2*n, and the arrays xm,wm dimension mc*ncapm+mp. c nd=2*n if(idelta.le.0) idelta=1 if(n.lt.1) then ierr=-1 return end if c c Initialization c incap=1 kount=-1 ierr=0 do 10 k=1,n beta(k)=0. 10 continue ncap=(nd-1)/idelta 20 do 30 k=1,n be(k)=beta(k) 30 continue kount=kount+1 if(kount.gt.1) incap=2**(kount/5)*n ncap=ncap+incap if(ncap.gt.ncapm) then ierr=ncapm return end if c c Discretization of the modified moments c mtncap=mc*ncap do 50 i=1,mc im1tn=(i-1)*ncap if(iq.eq.1) then call quad(ncap,x,w,i,ierr) else call qgp(ncap,x,w,i,ierr,mc,finl,finr,endl,endr,xfer,wfer) end if if(ierr.ne.0) then ierr=i return end if do 40 k=1,ncap xm(im1tn+k)=x(k) wm(im1tn+k)=w(k) 40 continue 50 continue if(mp.ne.0) then do 60 k=1,mp xm(mtncap+k)=xp(k) wm(mtncap+k)=yp(k) 60 continue end if mtnpmp=mtncap+mp do 90 k=1,nd km1=k-1 sum=0. do 80 i=1,mtnpmp p1=0. p=1. if(k.gt.1) then do 70 l=1,km1 pm1=p1 p1=p p=(xm(i)-a(l))*p1-b(l)*pm1 70 continue end if sum=sum+wm(i)*p 80 continue fnu(k)=sum 90 continue c c Computation of the desired recursion coefficients c call cheb_r4 (n,a,b,fnu,alpha,beta,s,ierr,s0,s1,s2) c c In the following statement, the absolute value of the beta's is c used to guard against failure in cases where the routine is applied c to variable-sign weight functions and hence the positivity of the c beta's is not guaranteed. c do 100 k=1,n if(abs(beta(k)-be(k)).gt.eps*abs(beta(k))) goto 20 100 continue return end subroutine mcdis(n,ncapm,mc,mp,xp,yp,quad,eps,iq,idelta,irout, *finl,finr,endl,endr,xfer,wfer,alpha,beta,ncap,kount,ierr,ie, *be,x,w,xm,wm,p0,p1,p2) c*********************************************************************72 c c This is a multiple-component discretization procedure as described in c Section 4.3 of the companion paper. It generates to a relative c accuracy of eps the recursion coefficients alpha(k), beta(k), c k=0,1,...,n-1, for the polynomials orthogonal with respect to a c weight distribution consisting of the sum of mc continuous c components and a discrete component with mp points. The continuous c part of the spectrum is made up of mc weight functions, each c supported on its own interval. These intervals may or may not be c disjoint. The discretization of the inner product on the i-th c interval is furnished either by a user-supplied subroutine quad, c or by the general-purpose subroutine qgp provided in this package, c depending on whether iq is equal, or not equal, to 1, respectively. c The user-supplied routine must have the form quad(n,x,w,i,ierr) and c is assumed to supply the abscissas x(k) and weights w(k), k=1,2, c ...,n, to be used in approximating the i-th inner product c c integral of p(x)*q(x)*wf(x,i)dx c c by the c c sum over k from 1 to n of w(k)*p(x(k))*q(x(k)), c c i=1,2,...,mc. c c The desired recurrence coefficients are then approximated by the c recursion coefficients of the discrete orthogonal polynomials c belonging to the discretized inner product, which in turn are c computed by either the Stieltjes procedure or the Lanczos algorithm c according as irout is equal to, or not equal to 1, respectively. c Two error flags ierr,ie are provided which signal the occurrence c of an error condition in the quadrature process, or in the routine c sti or lancz (whichever is used), respectively. The point spectrum c is given through its abscissas xp and jumps yp. c c If the quadrature routine quad has polynomial degree of exactness c at least id(n) for each i, and if id(n)/n = idelta + O(1/n) as c n goes to infinity, then the procedure is designed to converge after c one iteration, provided idelta is set with the appropriate c integer. Normally, idelta=1 (for interpolatory rules) or idelta=2 c (for Gaussian rules). The default value is idelta=1. c c Input: n - - the number of recursion coefficients desired; c type integer c ncapm - a discretization parameter indicating an upper c limit of the fineness of the discretization; c ncapm=500 will usually be satisfactory; type c integer c mc - - the number of disjoint intervals in the c continuous part of the spectrum; type integer c mp - - the number of points in the discrete part of c the spectrum; type integer. If there is no c point spectrum, set mp=0. c xp - - an array of dimension mp containing the c abscissas of the point spectrum c yp - - an array of dimension mp containing the jumps c of the point spectrum c quad - a subroutine determining the discretization of c the inner product on each component interval, c or a dummy routine if iq is not equal to 1 c (see below) c eps - - the desired relative accuracy of the nonzero c recursion coefficients; type real c iq - - an integer selecting a user-supplied quadrature c routine quad if iq=1 or the ORTHPOL routine c qgp otherwise c idelta - a nonzero integer, typically 1 or 2, inducing c fast convergence in the case of special quadrature c routines c irout - an integer selecting the routine for generating c the recursion coefficients from the discrete c inner product. Specifically, irout=1 selects the c routine sti, whereas any other value selects the c routine lancz c c The logical variables finl,finr, the arrays endl,endr of c dimension mc, and the arrays xfer,wfer of dimension ncapm are c input variables to the subroutine qgp and are used (and hence need c to be properly dimensioned) only if iq is not equal to 1. c c Output: alpha,beta - arrays of dimension n, holding as k-th c element alpha(k-1), beta(k-1), k=1,2,...,n, c respectively c ncap - an integer indicating the fineness of the c discretization that yields convergence within c the eps-tolerance c kount - the number of iterations used c ierr - an error flag, equal to 0 on normal return, c equal to -1 if n is not in the proper range, c equal to i if there is an error condition in c the discretization of the i-th interval, c and equal to ncapm if the discretized c Stieltjes procedure does not converge within the c discretization resolution specified by ncapm c ie - - an error flag inherited from the routine sti c or lancz (whichever is used) c c The array be of dimension n, the arrays x,w of dimension ncapm, c and the arrays xm,wm,p0,p1,p2 of dimension mc*ncapm+mp are used c for working space. The routine calls upon the subroutine sti or c lancz, depending on the choice of irout. c dimension xp(*),yp(*),endl(mc),endr(mc),xfer(ncapm),wfer(ncapm), *alpha(n),beta(n),be(n),x(ncapm),w(ncapm),xm(*),wm(*),p0(*),p1(*), *p2(*) logical finl,finr c c The arrays xp,yp are assumed to have dimension mp if mp > 0, and c the arrays xm,wm,p0,p1,p2 dimension mc*ncapm+mp. c if(idelta.le.0) idelta=1 if(n.lt.1) then ierr=-1 return end if c c Initialization c incap=1 kount=-1 ierr=0 do 10 k=1,n beta(k)=0. 10 continue ncap=(2*n-1)/idelta 20 do 30 k=1,n be(k)=beta(k) 30 continue kount=kount+1 if(kount.gt.1) incap=2**(kount/5)*n ncap=ncap+incap if(ncap.gt.ncapm) then ierr=ncapm return end if c c Discretization of the inner product c mtncap=mc*ncap do 50 i=1,mc im1tn=(i-1)*ncap if(iq.eq.1) then call quad(ncap,x,w,i,ierr) else call qgp(ncap,x,w,i,ierr,mc,finl,finr,endl,endr,xfer, * wfer) end if if(ierr.ne.0) then ierr=i return end if do 40 k=1,ncap xm(im1tn+k)=x(k) wm(im1tn+k)=w(k) 40 continue 50 continue if(mp.ne.0) then do 60 k=1,mp xm(mtncap+k)=xp(k) wm(mtncap+k)=yp(k) 60 continue end if c c Computation of the desired recursion coefficients c if(irout.eq.1) then call sti(n,mtncap+mp,xm,wm,alpha,beta,ie,p0,p1,p2) else call lancz(n,mtncap+mp,xm,wm,alpha,beta,ie,p0,p1) end if c c In the following statement, the absolute value of the beta's is c used to guard against failure in cases where the routine is applied c to variable-sign weight functions and hence the positivity of the c beta's is not guaranteed. c do 70 k=1,n if(abs(beta(k)-be(k)).gt.eps*abs(beta(k))) goto 20 70 continue return end function nu0her(n,z,eps) c*********************************************************************72 c c This is an auxiliary function routine providing a starting backward c recurrence index for the Hermite measure that can be used in place c of nu0 in the routines knum and dknum. c complex z nu0her = int ( 2.*(sqrt(.5*real(n+1))+.25*alog(1./eps)/ * abs(aimag(z)))**2 ) return end function nu0jac(n,z,eps) c*********************************************************************72 c c This is an auxiliary function routine providing a starting backward c recurrence index for the Jacobi measure that can be used in place c of nu0 in the routines knum and dknum. c complex z pi=4.*atan(1.) x=real(z) y=abs(aimag(z)) if(x.lt.1.) then if(x.lt.-1.) angle=.5*(2.*pi+atan(y/(x-1.))+atan(y/(x+1.))) if(x.eq.-1.) angle=.5*(1.5*pi-atan(.5*y)) if(x.gt.-1.) angle=.5*(pi+atan(y/(x-1.))+atan(y/(x+1.))) else if(x.eq.1.) angle=.5*(.5*pi+atan(.5*y)) if(x.gt.1.) angle=.5*(atan(y/(x-1.))+atan(y/(x+1.))) end if x2=x*x y2=y*y r=((x2-y2-1.)**2+4.*x2*y2)**.25 r=sqrt((x+r*cos(angle))**2+(y+r*sin(angle))**2) nu0jac = int ( real(n+1)+.5*alog(1./eps)/alog(r) ) return end function nu0lag(n,z,al,eps) c*********************************************************************72 c c This is an auxiliary function routine providing a starting backward c recurrence index for the Laguerre measure that can be used in place c of nu0 in the routines knum and dknum. c complex z pi=4.*atan(1.) x=real(z) y=aimag(z) phi=.5*pi if(y.lt.0.) phi=1.5*pi if(x.eq.0.) goto 10 phi=atan(y/x) if(y.gt.0. .and. x.gt.0.) goto 10 phi=phi+pi if(x.lt.0.) goto 10 phi=phi+pi 10 nu0lag = int ( (sqrt(real(n+1)+.5*(al+1.))+alog(1./eps)/(4.*(x*x+ * y*y)**.25*cos(.5*(phi-pi))))**2-.5*(al+1.) ) return end subroutine qgp(n,x,w,i,ierr,mc,finl,finr,endl,endr,xfer,wfer) c*********************************************************************72 c c This is a general-purpose discretization routine that can be used c as an alternative to the routine quad in the multiple-component c discretization procedure mcdis. It takes no account of the special c nature of the weight function involved and hence may result in slow c convergence of the discretization procedure. This routine, therefore, c should be used only as a last resort, when no better, more natural c discretization can be found. c c It is assumed that there are mc.ge.1 disjoint component intervals. c The discretization is effected by the Fejer quadrature rule, c suitably transformed to the respective interval. An interval that c extends to minus infinity has to be indexed by 1; one that extends c to plus infinity has to be indexed by mc. c c The output variable ierr is given the value 0. Additional input c parameters and working space used by this routine are as follows: c c mc - the number of component intervals; type integer c finl - a logical variable to be set .true. if the c extreme left interval is finite, and .false. c otherwise c finr - a logical variable to be set .true. if the c extreme right interval is finite, and .false. c otherwise c endl - an array of dimension mc containing the left c endpoints of the component intervals; if the c first of these extends to minus infinity, endl(1) c can be set to an arbitrary value c endr - an array of dimension mc containing the right c endpoints of the component intervals; if the c last of these extends to plus infinity, endr(mc) c can be set to an arbitrary value c xfer,wfer-working arrays holding the Fejer nodes and c weights, respectively, for the interval [-1,1]. c c The user has to supply the routine c c function wf(x,i), c c which evaluates the weight function at the point x on the i-th c component interval. The routine also uses the subroutines fejer, c symtr and tr, which are appended. c dimension x(n),w(n),endl(mc),endr(mc),xfer(*),wfer(*) logical finl,finr c c The arrays xfer,wfer are dimensioned in the routine mcdis. c ierr=0 if(i.eq.1) call fejer(n,xfer,wfer) if(i.gt.1 .and. i.lt.mc) goto 60 if(mc.eq.1) then if(finl.and.finr) goto 60 if(finl) goto 20 if(finr) goto 40 do 10 k=1,n call symtr(xfer(k),phi,phi1) x(k)=phi w(k)=wfer(k)*wf(phi,i)*phi1 10 continue return else if((i.eq.1.and.finl).or.(i.eq.mc.and.finr)) goto 60 if(i.eq.1) goto 40 end if 20 do 30 k=1,n call tr(xfer(k),phi,phi1) x(k)=endl(mc)+phi w(k)=wfer(k)*wf(x(k),mc)*phi1 30 continue return 40 do 50 k=1,n call tr(-xfer(k),phi,phi1) x(k)=endr(1)-phi w(k)=wfer(k)*wf(x(k),1)*phi1 50 continue return 60 do 70 k=1,n x(k)=.5*((endr(i)-endl(i))*xfer(k)+endr(i)+endl(i)) w(k)=.5*(endr(i)-endl(i))*wfer(k)*wf(x(k),i) 70 continue return end subroutine symtr(t,phi,phi1) c*********************************************************************72 c c This implements a particular transformation x=phi(t) mapping c the t-interval [-1,1] to the x-interval [-oo,oo]. c c input: t c output: phi=phi(t) c phi1=derivative of phi(t) c t2=t*t phi=t/(1.-t2) phi1=(t2+1.)/(t2-1.)**2 return end subroutine tr(t,phi,phi1) c*********************************************************************72 c c This implements a particular transformation x=phi(t) mapping c the t-interval [-1,1] to the x-interval [0,oo]. c c input: t c output: phi=phi(t) c phi1=derivative of phi(t) c phi=(1.+t)/(1.-t) phi1=2./(t-1.)**2 return end subroutine fejer(n,x,w) c*********************************************************************72 c c This routine generates the n-point Fejer quadrature rule. c c input: n - the number of quadrature nodes c output: x,w - arrays of dimension n holding the quadrature c nodes and weights, respectively; the nodes c are ordered increasingly c dimension x(n),w(n) pi=4.*atan(1.) nh=n/2 np1h=(n+1)/2 fn=real(n) do 10 k=1,nh x(n+1-k)=cos(.5*real(2*k-1)*pi/fn) x(k)=-x(n+1-k) 10 continue if(2*nh.ne.n) x(np1h)=0. do 30 k=1,np1h c1=1. c0=2.*x(k)*x(k)-1. t=2.*c0 sum=c0/3. do 20 m=2,nh c2=c1 c1=c0 c0=t*c1-c2 sum=sum+c0/real(4*m*m-1) 20 continue w(k)=2.*(1.-2.*sum)/fn w(n+1-k)=w(k) 30 continue return end subroutine radau(n,alpha,beta,end,zero,weight,ierr,e,a,b) c*********************************************************************72 c c Given n and a measure dlambda, this routine generates the c (n+1)-point Gauss-Radau quadrature formula c c integral over supp(dlambda) of f(t)dlambda(t) c c = w(0)f(x(0)) + sum from k=1 to k=n of w(k)f(x(k)) + R(n;f). c c The nodes are returned as zero(k)=x(k), the weights as weight(k) c =w(k), k=0,1,2,...,n. The user has to supply the recursion c coefficients alpha(k), beta(k), k=0,1,2,...,n, for the measure c dlambda. The nodes and weights are computed as eigenvalues and c in terms of the first component of the respective normalized c eigenvectors of a slightly modified Jacobi matrix of order n+1. c To do this, the routine calls upon the subroutine gauss. It also c uses the function subroutine r1mach. c c Input: n - - the number of interior points in the Gauss-Radau c formula; type integer c alpha,beta - arrays of dimension n+1 to be supplied with c the recursion coefficients alpha(k-1), beta(k-1), c k=1,2,...,n+1; the coefficient alpha(n+1) is not c used by the routine c end - the prescribed endpoint x(0) of the Gauss-Radau c formula; type real c c Output: zero - array of dimension n+1 containing the nodes (in c increasing order) zero(k)=x(k), k=0,1,2,...,n c weight-array of dimension n+1 containing the weights c weight(k)=w(k), k=0,1,2,...,n c ierr - an error flag inherited from the routine gauss c c The arrays e,a,b are needed for working space. c dimension alpha(*),beta(*),zero(*),weight(*),e(*),a(*),b(*) c c The arrays alpha,beta,zero,weight,e,a,b are assumed to have c dimension n+1. c epsma=r1mach(3) c c epsma is the machine single precision. c np1=n+1 do 10 k=1,np1 a(k)=alpha(k) b(k)=beta(k) 10 continue p0=0. p1=1. do 20 k=1,n pm1=p0 p0=p1 p1=(end-a(k))*p0-b(k)*pm1 20 continue a(np1)=end-b(np1)*p0/p1 call gauss(np1,a,b,epsma,zero,weight,ierr,e) return end subroutine recur(n,ipoly,al,be,a,b,ierr) c*********************************************************************72 c c This subroutine generates the coefficients a(k),b(k), k=0,1,...,n-1, c in the recurrence relation c c p(k+1)(x)=(x-a(k))*p(k)(x)-b(k)*p(k-1)(x), c k=0,1,...,n-1, c c p(-1)(x)=0, p(0)(x)=1, c c for some classical (monic) orthogonal polynomials, and sets b(0) c equal to the total mass of the weight distribution. The results are c stored in the arrays a,b, which hold, respectively, the coefficients c a(k-1),b(k-1), k=1,2,...,n. c c Input: n - - the number of recursion coefficients desired c ipoly-integer identifying the polynomial as follows: c 1=Legendre polynomial on (-1,1) c 2=Legendre polynomial on (0,1) c 3=Chebyshev polynomial of the first kind c 4=Chebyshev polynomial of the second kind c 5=Jacobi polynomial with parameters al=-.5,be=.5 c 6=Jacobi polynomial with parameters al,be c 7=generalized Laguerre polynomial with c parameter al c 8=Hermite polynomial c al,be-input parameters for Jacobi and generalized c Laguerre polynomials c c Output: a,b - arrays containing, respectively, the recursion c coefficients a(k-1),b(k-1), k=1,2,...,n. c ierr -an error flag, equal to 0 on normal return, c equal to 1 if al or be are out of range c when ipoly=6 or ipoly=7, equal to 2 if b(0) c overflows when ipoly=6 or ipoly=7, equal to 3 c if n is out of range, and equal to 4 if ipoly c is not an admissible integer. In the case ierr=2, c the coefficient b(0) is set equal to the largest c machine-representable number. c c The subroutine calls for the function subroutines r1mach,gamma and c alga. The routines gamma and alga, which are included in this file, c evaluate respectively the gamma function and its logarithm for c positive arguments. They are used only in the cases ipoly=6 and c ipoly=7. c external gamma dimension a(n),b(n) if(n.lt.1) then ierr=3 return end if almach=alog(r1mach(2)) ierr=0 do 10 k=1,n a(k)=0. 10 continue if(ipoly.eq.1) then b(1)=2. if (n.eq.1) return do 20 k=2,n fkm1=real(k-1) b(k)=1./(4.-1./(fkm1*fkm1)) 20 continue return else if (ipoly.eq.2) then a(1)=.5 b(1)=1. if(n.eq.1) return do 30 k=2,n a(k)=.5 fkm1=real(k-1) b(k)=.25/(4.-1./(fkm1*fkm1)) 30 continue return else if(ipoly.eq.3) then b(1)=4.*atan(1.) if(n.eq.1) return b(2)=.5 if(n.eq.2) return do 40 k=3,n b(k)=.25 40 continue return else if(ipoly.eq.4) then b(1)=2.*atan(1.) if(n.eq.1) return do 50 k=2,n b(k)=.25 50 continue return else if(ipoly.eq.5) then b(1)=4.*atan(1.) a(1)=.5 if(n.eq.1) return do 60 k=2,n b(k)=.25 60 continue return else if(ipoly.eq.6) then if(al.le.-1. .or. be.le.-1.) then ierr=1 return else alpbe=al+be a(1)=(be-al)/(alpbe+2.) t=(alpbe+1.)*alog(2.)+alga(al+1.)+alga(be+1.)- * alga(alpbe+2.) if(t.gt.almach) then ierr=2 b(1)=r1mach(2) else b(1)=exp(t) end if if(n.eq.1) return al2=al*al be2=be*be a(2)=(be2-al2)/((alpbe+2.)*(alpbe+4.)) b(2)=4.*(al+1.)*(be+1.)/((alpbe+3.)*(alpbe+2.)**2) if(n.eq.2) return do 70 k=3,n fkm1=real(k-1) a(k)=.25*(be2-al2)/(fkm1*fkm1*(1.+.5*alpbe/fkm1)* * (1.+.5*(alpbe+2.)/fkm1)) b(k)=.25*(1.+al/fkm1)*(1.+be/fkm1)*(1.+alpbe/fkm1)/ * ((1.+.5*(alpbe+1.)/fkm1)*(1.+.5*(alpbe-1.)/fkm1) * *(1.+.5*alpbe/fkm1)**2) 70 continue return end if else if(ipoly.eq.7) then if(al.le.-1.) then ierr=1 return else a(1)=al+1. b(1)=gamma(al+1.,ierr) if(ierr.eq.2) b(1)=r1mach(2) if(n.eq.1) return do 80 k=2,n fkm1=real(k-1) a(k)=2.*fkm1+al+1. b(k)=fkm1*(fkm1+al) 80 continue return end if else if(ipoly.eq.8) then b(1)=sqrt(4.*atan(1.)) if(n.eq.1) return do 90 k=2,n b(k)=.5*real(k-1) 90 continue return else ierr=4 end if end function alga(x) c*********************************************************************72 c c This is an auxiliary function subroutine (not optimized in any c sense) evaluating the logarithm of the gamma function for positive c arguments x. It is called by the subroutine gamma. The integer m0 c in the first executable statement is the smallest integer m such c that 1*3*5* ... *(2*m+1)/(2**m) is greater than or equal to the c largest machine-representable number. The routine is based on a c rational approximation valid on [.5,1.5] due to W.J. Cody and c K.E. Hillstrom; see Math. Comp. 21, 1967, 198-203, in particular the c case n=7 in Table II. For the computation of m0 it calls upon the c function subroutines t and r1mach. The former, appended below, c evaluates the inverse function t = t(y) of y = t ln t. c c Reference: c c William Cody, Kenneth Hillstrom, c Chebyshev Approximations for the Natural Logarithm of the c Gamma Function, c Mathematics of Computation, c Volume 21, Number 98, April 1967, pages 198-203. c dimension cnum(8),cden(8) data cnum/4.120843185847770,85.68982062831317,243.175243524421, *-261.7218583856145,-922.2613728801522,-517.6383498023218, *-77.41064071332953,-2.208843997216182/, *cden/1.,45.64677187585908,377.8372484823942,951.323597679706, *846.0755362020782,262.3083470269460,24.43519662506312, *.4097792921092615/ c c The constants in the statement below are exp(1.) and .5*alog(8.). c m0 = int ( 2.71828*t((alog(r1mach(2))-1.03972)/2.71828) ) xi=aint(x) if(x-xi.gt..5) xi=xi+1. m=ifix(xi)-1 c c Computation of log gamma on the standard interval (1/2,3/2] c xe=x-real(m) snum=cnum(1) sden=cden(1) do 10 k=2,8 snum=xe*snum+cnum(k) sden=xe*sden+cden(k) 10 continue alga=(xe-1.)*snum/sden c c Computation of log gamma on (0,1/2] c if(m.eq.-1) then alga=alga-alog(x) return else if(m.eq.0) then return else c c Computation of log gamma on (3/2,5/2] c p=xe if(m.eq.1) then alga=alga+alog(p) return else c c Computation of log gamma for arguments larger than 5/2 c mm1=m-1 c c The else-clause in the next statement is designed to avoid possible c overflow in the computation of p in the if-clause, at the expense c of computing many logarithms. c if(m.lt.m0) then do 20 k=1,mm1 p=(xe+real(k))*p 20 continue alga=alga+alog(p) return else alga=alga+alog(xe) do 30 k=1,mm1 alga=alga+alog(xe+real(k)) 30 continue return end if end if end if end function gamma(x,ierr) c*********************************************************************72 c c This evaluates the gamma function for real positive x, using the c function subroutines alga and r1mach. In case of overflow, the c routine returns the largest machine-representable number and the c error flag ierr=2. c almach=alog(r1mach(2)) ierr=0 t=alga(x) if(t.ge.almach) then ierr=2 gamma=r1mach(2) return else gamma=exp(t) return end if end function t(y) c*********************************************************************72 c c This evaluates the inverse function t = t(y) of y = t ln t for c nonnegative y to an accuracy of about one percent. For the c approximation used, see pp. 51-52 in W. Gautschi,Computational c aspects of three-term recurrence relations'', SIAM Rev. 9, 1967, c 24-82. c if(y.le.10.) then p=.000057941*y-.00176148 p=y*p+.0208645 p=y*p-.129013 p=y*p+.85777 t=y*p+1.0125 else z=alog(y)-.775 p=(.775-alog(z))/(1.+z) p=1./(1.+p) t=y*p/z end if return end subroutine sti(n,ncap,x,w,alpha,beta,ierr,p0,p1,p2) c*********************************************************************72 c c This routine applies Stieltjes's procedure'' to generate the recursion c coefficients alpha(k), beta(k) , k=0,1,...,n-1, for the discrete c (monic) orthogonal polynomials associated with the inner product c c (f,g)=sum over k from 1 to ncap of w(k)*f(x(k))*g(x(k)). c c The integer n must be between 1 and ncap, inclusive; otherwise, c there is an error exit with ierr=1. The results are stored in the c arrays alpha, beta; the arrays p0, p1, p2 are working arrays. c c If there is a threat of underflow or overflow in the calculation c of the coefficients alpha(k) and beta(k), the routine exits with c the error flag ierr set equal to -k (in the case of underflow) c or +k (in the case of overflow), where k is the recursion index c for which the problem occurs. The former [latter] can often be avoided c by multiplying all weights w(k) by a sufficiently large [small] c scaling factor prior to entering the routine, and, upon exit, divide c the coefficient beta(0) by the same factor. c c This routine should be used with caution if n is relatively close c to ncap, since there is a distinct possibility of numerical c instability developing. (See W. Gautschi,Is the recurrence relation c for orthogonal polynomials always stable?'', BIT, 1993, to appear.) c In that case, the routine lancz should be used. c c Reference: c c Walter Gautschi, c On Generating Orthogonal Polynomials, c SIAM Journal on Scientific and Statistical Computing, c Volume 3, Number 3, 1982, pages 289-317. c dimension x(ncap),w(ncap),alpha(n),beta(n),p0(ncap),p1(ncap), *p2(ncap) tiny=10.*r1mach(1) huge=.1*r1mach(2) ierr=0 if(n.le.0 .or. n.gt.ncap) then ierr=1 return end if nm1=n-1 c c Compute the first alpha- and beta-coefficient. c sum0=0. sum1=0. do 10 m=1,ncap sum0=sum0+w(m) sum1=sum1+w(m)*x(m) 10 continue alpha(1)=sum1/sum0 beta(1)=sum0 if(n.eq.1) return c c Compute the remaining alpha- and beta-coefficients. c do 20 m=1,ncap p1(m)=0. p2(m)=1. 20 continue do 40 k=1,nm1 sum1=0. sum2=0. do 30 m=1,ncap c c The following statement is designed to avoid an overflow condition c in the computation of p2(m) when the weights w(m) go to zero c faster (and underflow) than the p2(m) grow. c if(w(m).eq.0.) goto 30 p0(m)=p1(m) p1(m)=p2(m) p2(m)=(x(m)-alpha(k))*p1(m)-beta(k)*p0(m) c c Check for impending overflow. c if(abs(p2(m)).gt.huge .or. abs(sum2).gt.huge) then ierr=k return end if t=w(m)*p2(m)*p2(m) sum1=sum1+t sum2=sum2+t*x(m) 30 continue c c Check for impending underflow. c if(abs(sum1).lt.tiny) then ierr=-k return end if alpha(k+1)=sum2/sum1 beta(k+1)=sum1/sum0 sum0=sum1 40 continue return end function d1mach ( i ) c*********************************************************************72 c cc D1MACH returns double precision real machine-dependent constants. c c Discussion: c c D1MACH can be used to obtain machine-dependent parameters c for the local machine environment. It is a function c with one input argument, and can be called as follows: c c D = D1MACH ( I ) c c where I=1,...,5. The output value of D above is c determined by the input value of I:. c c D1MACH ( 1) = B**(EMIN-1), the smallest positive magnitude. c D1MACH ( 2) = B**EMAX*(1 - B**(-T)), the largest magnitude. c D1MACH ( 3) = B**(-T), the smallest relative spacing. c D1MACH ( 4) = B**(1-T), the largest relative spacing. c D1MACH ( 5) = LOG10(B) c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 25 April 2007 c c Author: c c Original FORTRAN77 version by Phyllis Fox, Andrew Hall, Norman Schryer. c This version by John Burkardt. c c Reference: c c Phyllis Fox, Andrew Hall, Norman Schryer, c Algorithm 528: c Framework for a Portable Library, c ACM Transactions on Mathematical Software, c Volume 4, Number 2, June 1978, page 176-188. c c Parameters: c c Input, integer I, the index of the desired constant. c c Output, double precision D1MACH, the value of the constant. c implicit none double precision d1mach integer i if ( i < 1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'D1MACH - Fatal error!' write ( *, '(a)' ) ' The input argument I is out of bounds.' write ( *, '(a)' ) ' Legal values satisfy 1 <= I <= 5.' write ( *, '(a,i12)' ) ' I = ', i d1mach = 0.0D+00 stop else if ( i == 1 ) then d1mach = 4.450147717014403D-308 else if ( i == 2 ) then d1mach = 8.988465674311579D+307 else if ( i == 3 ) then d1mach = 1.110223024625157D-016 else if ( i == 4 ) then d1mach = 2.220446049250313D-016 else if ( i == 5 ) then d1mach = 0.301029995663981D+000 else if ( 5 < i ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'D1MACH - Fatal error!' write ( *, '(a)' ) ' The input argument I is out of bounds.' write ( *, '(a)' ) ' Legal values satisfy 1 <= I <= 5.' write ( *, '(a,i12)' ) ' I = ', i d1mach = 0.0D+00 stop end if return end function i1mach ( i ) c*********************************************************************72 c cc I1MACH returns integer machine dependent constants. c c Discussion: c c Input/output unit numbers. c c I1MACH(1) = the standard input unit. c I1MACH(2) = the standard output unit. c I1MACH(3) = the standard punch unit. c I1MACH(4) = the standard error message unit. c c Words. c c I1MACH(5) = the number of bits per integer storage unit. c I1MACH(6) = the number of characters per integer storage unit. c c Integers. c c Assume integers are represented in the S digit base A form: c c Sign * (X(S-1)*A**(S-1) + ... + X(1)*A + X(0)) c c where 0 <= X(1:S-1) < A. c c I1MACH(7) = A, the base. c I1MACH(8) = S, the number of base A digits. c I1MACH(9) = A**S-1, the largest integer. c c Floating point numbers c c Assume floating point numbers are represented in the T digit c base B form: c c Sign * (B**E) * ((X(1)/B) + ... + (X(T)/B**T) ) c c where 0 <= X(I) < B for I=1 to T, 0 < X(1) and EMIN <= E <= EMAX. c c I1MACH(10) = B, the base. c c Single precision c c I1MACH(11) = T, the number of base B digits. c I1MACH(12) = EMIN, the smallest exponent E. c I1MACH(13) = EMAX, the largest exponent E. c c Double precision c c I1MACH(14) = T, the number of base B digits. c I1MACH(15) = EMIN, the smallest exponent E. c I1MACH(16) = EMAX, the largest exponent E. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 25 April 2007 c c Author: c c Original FORTRAN77 version by Phyllis Fox, Andrew Hall, Norman Schryer. c This version by John Burkardt. c c Reference: c c Phyllis Fox, Andrew Hall, Norman Schryer, c Algorithm 528, c Framework for a Portable Library, c ACM Transactions on Mathematical Software, c Volume 4, Number 2, June 1978, page 176-188. c c Parameters: c c Input, integer I, chooses the parameter to be returned. c 1 <= I <= 16. c c Output, integer I1MACH, the value of the chosen parameter. c implicit none integer i integer i1mach if ( i < 1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I1MACH - Fatal error!' write ( *, '(a)' ) ' The input argument I is out of bounds.' write ( *, '(a)' ) ' Legal values satisfy 1 <= I <= 16.' write ( *, '(a,i12)' ) ' I = ', i i1mach = 0 stop else if ( i == 1 ) then i1mach = 5 else if ( i == 2 ) then i1mach = 6 else if ( i == 3 ) then i1mach = 7 else if ( i == 4 ) then i1mach = 6 else if ( i == 5 ) then i1mach = 32 else if ( i == 6 ) then i1mach = 4 else if ( i == 7 ) then i1mach = 2 else if ( i == 8 ) then i1mach = 31 else if ( i == 9 ) then i1mach = 2147483647 else if ( i == 10 ) then i1mach = 2 else if ( i == 11 ) then i1mach = 24 else if ( i == 12 ) then i1mach = -125 else if ( i == 13 ) then i1mach = 128 else if ( i == 14 ) then i1mach = 53 else if ( i == 15 ) then i1mach = -1021 else if ( i == 16 ) then i1mach = 1024 else if ( 16 < i ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I1MACH - Fatal error!' write ( *, '(a)' ) ' The input argument I is out of bounds.' write ( *, '(a)' ) ' Legal values satisfy 1 <= I <= 16.' write ( *, '(a,i12)' ) ' I = ', i i1mach = 0 stop end if return end function r1mach ( i ) c*********************************************************************72 c cc R1MACH returns single precision real machine constants. c c Discussion: c c Assume that single precision real numbers are stored with a mantissa c of T digits in base B, with an exponent whose value must lie c between EMIN and EMAX. Then for values of I between 1 and 5, c R1MACH will return the following values: c c R1MACH(1) = B**(EMIN-1), the smallest positive magnitude. c R1MACH(2) = B**EMAX*(1-B**(-T)), the largest magnitude. c R1MACH(3) = B**(-T), the smallest relative spacing. c R1MACH(4) = B**(1-T), the largest relative spacing. c R1MACH(5) = log10(B) c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 25 April 2007 c c Author: c c Original FORTRAN77 version by Phyllis Fox, Andrew Hall, Norman Schryer. c This version by John Burkardt. c c Reference: c c Phyllis Fox, Andrew Hall, Norman Schryer, c Algorithm 528, c Framework for a Portable Library, c ACM Transactions on Mathematical Software, c Volume 4, Number 2, June 1978, page 176-188. c c Parameters: c c Input, integer I, chooses the parameter to be returned. c 1 <= I <= 5. c c Output, real R1MACH, the value of the chosen parameter. c implicit none integer i real r1mach if ( i < 1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R1MACH - Fatal error!' write ( *, '(a)' ) ' The input argument I is out of bounds.' write ( *, '(a)' ) ' Legal values satisfy 1 <= I <= 5.' write ( *, '(a,i12)' ) ' I = ', i r1mach = 0.0E+00 stop else if ( i == 1 ) then r1mach = 1.1754944E-38 else if ( i == 2 ) then r1mach = 3.4028235E+38 else if ( i == 3 ) then r1mach = 5.9604645E-08 else if ( i == 4 ) then r1mach = 1.1920929E-07 else if ( i == 5 ) then r1mach = 0.3010300E+00 else if ( 5 < i ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R1MACH - Fatal error!' write ( *, '(a)' ) ' The input argument I is out of bounds.' write ( *, '(a)' ) ' Legal values satisfy 1 <= I <= 5.' write ( *, '(a,i12)' ) ' I = ', i r1mach = 0.0E+00 stop end if return end subroutine timestamp ( ) c*********************************************************************72 c cc TIMESTAMP prints out the current YMDHMS date as a timestamp. c c Discussion: c c This FORTRAN77 version is made available for cases where the c FORTRAN90 version cannot be used. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 12 January 2007 c c Author: c c John Burkardt c c Parameters: c c None c implicit none character * ( 8 ) ampm integer d character * ( 8 ) date integer h integer m integer mm character * ( 9 ) month(12) integer n integer s character * ( 10 ) time integer y save month data month / & 'January ', 'February ', 'March ', 'April ', & 'May ', 'June ', 'July ', 'August ', & 'September', 'October ', 'November ', 'December ' / call date_and_time ( date, time ) read ( date, '(i4,i2,i2)' ) y, m, d read ( time, '(i2,i2,i2,1x,i3)' ) h, n, s, mm if ( h .lt. 12 ) then ampm = 'AM' else if ( h .eq. 12 ) then if ( n .eq. 0 .and. s .eq. 0 ) then ampm = 'Noon' else ampm = 'PM' end if else h = h - 12 if ( h .lt. 12 ) then ampm = 'PM' else if ( h .eq. 12 ) then if ( n .eq. 0 .and. s .eq. 0 ) then ampm = 'Midnight' else ampm = 'AM' end if end if end if write ( *, & '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) & d, month(m), y, h, ':', n, ':', s, '.', mm, ampm return end