      program main

c*********************************************************************72
c
cc toms726_test() tests toms726().
c
c  Modified:
c
c    26 January 2007
c
c  Author:
c
c    John Burkardt
c
      implicit none

      call timestamp ( )
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'toms726_test():'
      write ( *, '(a)' ) '  FORTRAN77 version'
      write ( *, '(a)' ) '  Test toms726().'

      call test01 ( )
      call test02 ( )
      call test03 ( )
c
c  This test is causing floating point exceptions.
c
      if ( .false. ) then
        call test04 ( )
      end if
      call test05 ( )
      call test06 ( )
      call test07 ( )
      call test08 ( )
      call test09 ( )

      call test10 ( )
      call test11 ( )
c
c  Terminate.
c
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'toms726_test():'
      write ( *, '(a)' ) '  Normal end of execution.'
      write ( *, '(a)' ) ' '
      call timestamp ( )

      stop
      end
      subroutine test01 ( )

c*********************************************************************72
c
cc TEST01
c
c This test generates the first n beta-coefficients in the recurrence
c relation for the orthogonal polynomials relative to the weight
c function
c
c         ((1-om2*x**2)*(1-x**2))**(-1/2)  on (-1,1)
c
c for om2=.1(.2).9,.99,.999, both in single and double precision,
c using modified moments if  modmom=.true.  and ordinary moments
c otherwise. In the former case, n=80, in the latter, n=20. Printed
c are the double-precision values of the coefficients along with the
c relative errors of the single-precision values.
c
c  Modified:
c
c    17 March 2008
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 fnu(160),f(80),f0(80),rr(80),a(159),b(159),
     *alpha(80),beta(80),s(80),s0(160),s1(160),s2(160)
      double precision doom2(7),deps,d1mach,dom2,dnu(160),d(80),
     *d0(80),drr(80),da(159),db(159),dalpha(80),dbeta(80),ds(80),
     *ds0(160),ds1(160),ds2(160)
      logical modmom

      data doom2/.1d0,.3d0,.5d0,.7d0,.9d0,.99d0,.999d0/

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TEST01'

      write(*,1)
    1 format(/)
      modmom=.true.
      eps=r1mach(3)
      deps=d1mach(3)
      if(modmom) then
        n=80
      else
        n=20
      end if
      ndm1=2*n-1
      do 30 iom=1,7
        dom2=doom2(iom)
        om2=sngl(dom2)
c
c Compute the modified resp. ordinary moments using Eqs. (3.7) and (3.9)
c of the companion paper. On machines with limited exponent range, some
c of the high-order modified moments may underflow, without this having
c any deteriorating effect on the accuracy.
c
        call mm_01_r4(n,eps,modmom,om2,fnu,ierr,f,f0,rr)
        call mm_01_r8(n,deps,modmom,dom2,dnu,iderr,d,d0,drr)
        if(ierr.ne.0 .or. iderr.ne.0) then
          write(*,2) ierr,iderr,om2
    2     format(/5x,'ierr in mm_01_r4 = ',i1,
     &    '  iderr in mm_01_r8 = ',i1,
     *      '  for om2 = ',f8.4/)
          goto 30
        end if
c
c Generate the recursion coefficients for the polynomials defining the
c modified resp. ordinary moments.
c
        if(modmom) then
          call recur(ndm1,3,0.,0.,a,b,ierr)
          call drecur(ndm1,3,0.d0,0.d0,da,db,iderr)
        else
          do k=1,ndm1
            a(k)=0.
            b(k)=0.
            da(k)=0.d0
            db(k)=0.d0
          end do
        end if
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,ierr,s0,s1,s2)
c
c On machines with limited single-precision exponent range, the routine
c CHEB_R4  may generate an underflow exception, which however is harmless
c and can be ignored.
c
        call cheb_r8(n,da,db,dnu,dalpha,dbeta,ds,iderr,ds0,ds1,ds2)
        write(*,3) ierr,iderr
    3   format(/5x,'ierr in cheb_r4 = ',i3,'  iderr in cheb_r8 = ',i3/)
        write(*,4)
    4   format(/5x,'k',14x,'dbeta(k)'/)
        do 20 k=1,n
          km1=k-1
          if(iderr.eq.0 .or. km1.lt.abs(iderr)) then
            if(ierr.eq.0 .or. km1.lt.abs(ierr)) then
              errb=sngl(dabs(dble(beta(k))-dbeta(k))/dbeta(k))
              if(k.eq.1) then
                write(*,5) km1,dbeta(k),errb,om2
    5           format(1x,i5,d36.28,e12.4,'   om2 =',f6.3)
              else
                write(*,6) km1,dbeta(k),errb
    6           format(1x,i5,d36.28,e12.4)
              end if
            else
              write(*,7) km1,dbeta(k)
    7         format(1x,i5,d36.28)
            end if
          end if
   20   continue
        write(*,1)
   30 continue
      return
      end
      subroutine mm_01_r4(n,eps,modmom,om2,fnu,ierr,f,f0,rr)

c*********************************************************************72
c
cc MM_01_R4 generates the modified (Chebyshev) resp. ordinary
c moments of the weight function
c
c          ((1-om2*x**2)*(1-x**2))**(-1/2)  on (-1,1)
c
c using Eqs. (3.7) resp. (3.9) of the companion paper.
c
      dimension fnu(*),f(n),f0(n),rr(n)
      logical modmom
c
c The array  fnu  is assumed to have dimension  2*n.
c
      ierr=0
      nd=2*n
      ndm1=nd-1
      pi=4.*atan(1.)
c
c Compute the Fourier coefficients of ((1-om2*sin(theta)**2))**(-1/2)
c as minimal solution of a three-term recurrence relation as described
c on pp.310-311 of W. Gautschi,On generating orthogonal polynomials'',
c SIAM J. Sci. Statist. Comput. 3, 1982, 289-317.
c
      q=om2/(2.-om2+2.*sqrt(1.-om2))
      q1=(1.+q*q)/q
      do k=1,n
        f(k)=0.
      end do
      nu=nd
   20 nu=nu+10
      do 30 k=1,n
        f0(k)=f(k)
   30 continue
      if(nu.gt.500) then
        ierr=1
        return
      end if
      r=0.
      s=0.
      do 40 k=1,nu
        n1=nu-k+1
        fn1=real(n1)
        r=-(fn1-.5)/(fn1*q1+(fn1+.5)*r)
        s=r*(2.+s)
        if(n1.le.n) rr(n1)=r
   40 continue
      c0=1./(1.+s)
      f(1)=rr(1)*c0
      if(n.gt.1) then
        do 50 k=2,n
          f(k)=rr(k)*f(k-1)
   50   continue
      end if
      do 60 k=1,n
        if(abs(f(k)-f0(k)).gt.eps*abs(f(k))) goto 20
   60 continue
c
c Compute the desired modified resp. ordinary moments in term of
c the above Fourier coefficients.
c
      fnu(1)=pi*c0
      if(n.eq.1) return
      fnu(2)=0.
      if(n.eq.2) return
      if(modmom) then
        c=2.*pi
        do 70 k=3,ndm1,2
          k1=(k-1)/2
          c=-.25*c
          fnu(k)=c*f(k1)
          fnu(k+1)=0.
   70   continue
      else
        c=.5*pi
        fnu(3)=c*(c0-f(1))
        fnu(4)=0.
        c=-c
        do 90 k=5,ndm1,2
          k1=(k-1)/2
          k1m1=k1-1
          c=-.25*c
          c1=1.
          sum=f(k1)
          do 80 i=1,k1m1
            c1=-c1*real(2*k1-i+1)/real(i)
            sum=sum+c1*f(k1-i)
   80     continue
          c1=-c1*real(k1+1)/real(2*k1)
          sum=sum+c1*c0
          fnu(k)=c*sum
          fnu(k+1)=0.
   90   continue
      end if

      return
      end
      subroutine mm_01_r8(n,deps,modmom,dom2,dnu,ierrd,d,d0,drr)

c*********************************************************************72
c
cc MM_01_R8 is a double-precision version of the routine MM_01_R4.
c
      double precision deps,dom2,dnu(*),d(n),d0(n),drr(n),dpi,dq,
     *dq1,dr,ds,dn1,dc0,dc,dc1,dsum
      logical modmom
c
c The array  dnu  is assumed to have dimension  2*n.
c
      ierrd=0
      nd=2*n
      ndm1=nd-1
      dpi=4.d0*datan(1.d0)
      dq=dom2/(2.d0-dom2+2.d0*dsqrt(1.d0-dom2))
      dq1=(1.d0+dq*dq)/dq
      do k=1,n
        d(k)=0.d0
      end do
      nud=nd
   20 nud=nud+10
      do 30 k=1,n
        d0(k)=d(k)
   30 continue
      if(nud.gt.1000) then
        ierrd=1
        return
      end if
      dr=0.d0
      ds=0.d0
      do 40 k=1,nud
        n1=nud-k+1
        dn1=dble(n1)
        dr=-(dn1-.5d0)/(dn1*dq1+(dn1+.5d0)*dr)
        ds=dr*(2.d0+ds)
        if(n1.le.n) drr(n1)=dr
   40 continue
      dc0=1.d0/(1.d0+ds)
      d(1)=drr(1)*dc0
      if(n.gt.1) then
        do 50 k=2,n
          d(k)=drr(k)*d(k-1)
   50   continue
      end if
      do 60 k=1,n
        if(dabs(d(k)-d0(k)).gt.deps*dabs(d(k))) goto 20
   60 continue
      dnu(1)=dpi*dc0
      if(n.eq.1) return
      dnu(2)=0.d0
      if(n.eq.2) return
      if(modmom) then
        dc=2.d0*dpi
        do 70 k=3,ndm1,2
          k1=(k-1)/2
          dc=-.25d0*dc
          dnu(k)=dc*d(k1)
          dnu(k+1)=0.d0
   70   continue
      else
        dc=.5d0*dpi
        dnu(3)=dc*(dc0-d(1))
        dnu(4)=0.d0
        dc=-dc
        do 90 k=5,ndm1,2
          k1=(k-1)/2
          k1m1=k1-1
          dc=-.25d0*dc
          dc1=1.d0
          dsum=d(k1)
          do 80 i=1,k1m1
            dc1=-dc1*dble(2*k1-i+1)/dble(i)
            dsum=dsum+dc1*d(k1-i)
   80     continue
          dc1=-dc1*dble(k1+1)/dble(2*k1)
          dsum=dsum+dc1*dc0
          dnu(k)=dc*dsum
          dnu(k+1)=0.d0
   90   continue
      end if

      return
      end
      subroutine test02 ( )

c*********************************************************************72
c
c This test generates the first n recursion coefficients for the
c orthogonal polynomials relative to the weight function
c
c        (x**sigma)*ln(1/x)  on (0,1],  sigma = -.5, 0, .5,
c
c where n=100 when using modified (Legendre) moments, and n=12 when
c using ordinary moments. It prints the double-precision values of the
c coefficients as well as the relative errors of the respective single-
c precision values and the maximum relative errors.
c
      dimension a(199),b(199),fnu(200),alpha(100),beta(100),s(100),
     *s0(200),s1(200),s2(200)
      double precision dsigma,da(199),db(199),dnu(200),dalpha(100),
     *dbeta(100),ds(100),ds0(200),ds1(200),ds2(200)
      logical modmom,intexp

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TEST02'

      modmom=.true.
c
c Generate the recursion coefficients for the polynomials defining the
c modified resp. ordinary moments.
c
      if(modmom) then
        n=100
        ndm1=2*n-1
        call recur(ndm1,2,0.,0.,a,b,ierr)
        call drecur(ndm1,2,0.d0,0.d0,da,db,iderr)
      else
        n=12
        ndm1=2*n-1
        do k=1,ndm1
          a(k)=0.
          b(k)=0.
          da(k)=0.d0
          db(k)=0.d0
        end do
      end if
      do 30 is=1,3
        dsigma=-.5d0+.5d0*dble(is-1)
        sigma=sngl(dsigma)
        if(is.eq.2) then
          intexp=.true.
        else
          intexp=.false.
        end if
c
c Compute the modified resp. ordinary moments using Eqs. (3.12) and
c (3.11) of the companion paper. On machines with limited exponent
c range, some of the high-order modified moments may underflow, without
c this having any deteriorating effect on the accuracy.
c
        call mm_02_r4(n,modmom,intexp,sigma,fnu)
        call mm_02_r8(n,modmom,intexp,dsigma,dnu)
c
c Compute the desired recursion coefficients by means of the modified
c Chebyshev algorithm; for the latter, see, e.g., Section 2.4 of
c W. Gautschi, On generating orthogonal polynomials'', SIAM J. Sci.
c Statist. Comput. 3, 1982, 289-317.
c
        call cheb_r4(n,a,b,fnu,alpha,beta,s,ierr,s0,s1,s2)
c
c On machines with limited single-precision exponent range, the routine
c CHEB_R4  may generate an underflow exception, which however is harmless
c and can be ignored.
c
        call cheb_r8(n,da,db,dnu,dalpha,dbeta,ds,iderr,ds0,ds1,ds2)
        write(*,1) ierr,iderr
    1   format(/6x,'ierr in cheb_r4 = ',i4,' iderr in cheb_r8 = ',i4/)
c
c Compute and print the relative errors and their maxima.
c
        eamax=0.
        ebmax=0.
        write(*,2) sigma
    2   format(/32x,'sigma =',f5.1)
        write(*,3)
    3   format(/3x,'k',13x,'dalpha(k)',25x,'dbeta(k)'/)
        do 20 k=1,n
          km1=k-1
          if(iderr.eq.0 .or. km1.lt.abs(iderr)) then
            write(*,4) km1,dalpha(k),dbeta(k)
    4       format(1x,i3,2d33.25)
            if(ierr.eq.0 .or. km1.lt.abs(ierr)) then
              erra=sngl(dabs(dble(alpha(k))-dalpha(k))/dalpha(k))
              errb=sngl(dabs(dble(beta(k))-dbeta(k))/dbeta(k))
              write(*,5) erra,errb
    5         format(4x,e12.4,21x,e12.4)
              if(erra.gt.eamax) then
                eamax=erra
                kamax=km1
              end if
              if(errb.gt.ebmax) then
                ebmax=errb
                kbmax=km1
              end if
            end if
          end if
   20   continue
        write(*,6) eamax,kamax,ebmax,kbmax
    6   format(/6x,'eamax =',e11.4,' at',i3,9x,'ebmax =',e11.4,
     *    ' at',i3//)
   30 continue

      return
      end
      subroutine mm_02_r4(n,modmom,intexp,sigma,fnu)

c*********************************************************************72
c
c This generates the first  2*n  modified moments (if modmom=.true.)
c relative to shifted monic Legendre polynomials, using Eq. (3.12) of
c the companion paper, and the first  2*n  ordinary moments (if modmom
c =.false.) by Eq. (3.11), of the weight function
c
c          (x**sigma)*ln(1/x)  on (0,1],   sigma > -1,
c
c for sigma an integer (if intexp=.true.) or a real number (if intexp
c =.false.). In either case, the input variable  sigma  is of type real.
c
      dimension fnu(*)
      logical modmom,intexp
c
c The array  fnu  is assumed to have dimension  2*n.
c
      nd=2*n
      sigp1=sigma+1.
      if(modmom) then
        isigma=int(sigma)
        isigp1=isigma+1
        isigp2=isigma+2
        isigp3=isigma+3
        if(intexp .and. isigp1.lt.nd) then
          kmax=isigp1
        else
          kmax=nd
        end if
        c=1.
        do 20 k=1,kmax
          km1=k-1
          fk=real(k)
          p=1.
          s=1./sigp1
          if(kmax.gt.1) then
            do i=1,km1
              fi=real(i)
              p=(sigp1-fi)*p/(sigp1+fi)
              s=s+1./(sigp1+fi)-1./(sigp1-fi)
            end do
          end if
          fnu(k)=c*s*p/sigp1
          c=fk*c/(4.*fk-2.)
   20   continue
        if(.not.intexp .or. isigp1.ge.nd) return
        q=-.5
        if(isigma.gt.0) then
          do 30 iq=1,isigma
            fiq=real(iq)
            q=fiq*fiq*q/((2.*fiq+1.)*(2.*fiq+2.))
   30     continue
        end if
        fnu(isigp2)=c*q
        if(isigp2.eq.nd) return
        do 40 k=isigp3,nd
          km1=k-1
          fkm1=real(km1)
          fnu(k)=-fkm1*(fkm1-sigp1)*fnu(km1)/((4.*fkm1-2.)*
     *      (fkm1+sigp1))
   40   continue
        return
      else
        do 50 k=1,nd
          fkm1=real(k-1)
          fnu(k)=(1./(sigp1+fkm1))**2
   50   continue
      end if
      end
      subroutine mm_02_r8(n,modmom,intexp,dsigma,dnu)

c*********************************************************************72
c
cc MM_02_R8 is a double-precision version of the routine MM_02_R4.
c
      double precision dsigma,dnu(*),dsigp1,dc,dk,dp,ds,di,dq,diq,dkm1
      logical modmom,intexp
c
c The array  dnu  is assumed to have dimension  2*n.
c
      nd=2*n
      dsigp1=dsigma+1.d0
      if(modmom) then
        isigma=idint(dsigma)
        isigp1=isigma+1
        isigp2=isigma+2
        isigp3=isigma+3
        if(intexp .and. isigp1.lt.nd) then
          kmax=isigp1
        else
          kmax=nd
        end if
        dc=1.d0
        do 20 k=1,kmax
          km1=k-1
          dk=dble(k)
          dp=1.d0
          ds=1.d0/dsigp1
          if(kmax.gt.1) then
            do i=1,km1
              di=dble(i)
              dp=(dsigp1-di)*dp/(dsigp1+di)
              ds=ds+1.d0/(dsigp1+di)-1.d0/(dsigp1-di)
            end do
          end if
          dnu(k)=dc*ds*dp/dsigp1
          dc=dk*dc/(4.d0*dk-2.d0)
   20   continue
        if(.not.intexp .or. isigp1.ge.nd) return
        dq=-.5d0
        if(isigma.gt.0) then
          do 30 iq=1,isigma
            diq=dble(iq)
            dq=diq*diq*dq/((2.d0*diq+1.d0)*(2.d0*diq+2.d0))
   30     continue
        end if
        dnu(isigp2)=dc*dq
        if(isigp2.eq.nd) return
        do 40 k=isigp3,nd
          km1=k-1
          dkm1=dble(km1)
          dnu(k)=-dkm1*(dkm1-dsigp1)*dnu(km1)/((4.d0*dkm1-2.d0)*
     *      (dkm1+dsigp1))
   40   continue
        return
      else
        do 50 k=1,nd
          dkm1=dble(k-1)
          dnu(k)=(1.d0/(dsigp1+dkm1))**2
   50   continue
      end if
      return
      end
      subroutine test03 ( )

c*********************************************************************72
c
      dimension x(320),w(320),bexact(320),p0(320),p1(320),p2(320),
     *als(320),bes(320),all(320),bel(320)
c
c This test applies both the Stieltjes procedure (cf. Section 2.1 of
c W. Gautschi, On generating orthogonal polynomials'', SIAM J. Sci.
c Statist. Comput. 3, 1982, 289-317) and the Lanczos algorithm (cf.
c W.B. Gragg and W.J. Harrod, The numerically stable reconstruction of
c Jacobi matrices from spectral data'', Numer. Math. 44, 1984, 317-335)
c to generate the  N  recursion coefficients, N = 40, 80, 160, 320,
c for the (monic) orthogonal polynomials relative to the discrete inner
c product supported on  N  equally spaced points on [-1,1] (including
c the end points) and having equal weights 2/N. The routine prints the
c absolute errors in the alpha-coefficients and the relative errors in
c the beta-coefficients, these being computed using known formulae for
c the coefficients. The maxima of these errors are also printed.
c
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TEST03'

      ncap=20
      do 30 incap=1,4
        ncap=2*ncap
        fncap=real(ncap)
        fncm1=real(ncap-1)
c
c Generate the abscissae and weights of the discrete inner product and
c the exact beta-coefficients.
c
        x(1)=-1.
        w(1)=2./fncap
        bexact(1)=2.
        do k=2,ncap
          fkm1=real(k-1)
          x(k)=-1.+2.*fkm1/fncm1
          w(k)=w(1)
          bexact(k)=(1.+1./fncm1)**2*(1.-(fkm1/fncap)**2)/
     *      (4.-(1./fkm1)**2)
        end do
c
c Compute the desired coefficients, first by the Stieltjes procedure,
c and then by the Lanczos algorithm. Indicate via the error flag  ierrs
c whether a critical underflow condition has arisen in Stieltjes's
c procedure. (There may, in addition, occur harmless underflow, which
c the routine  sti  does not test for.)
c
        call sti(ncap,ncap,x,w,als,bes,ierrs,p0,p1,p2)
        call lancz(ncap,ncap,x,w,all,bel,ierrl,p0,p1)
        write(*,1) ierrs,ierrl
    1   format(/5x,'ierr in sti = ',i4,11x,'ierr in lancz = ',i3/)
c
c Compute and print the absolute errors of the alpha-coefficients and
c the relative errors of the beta-coefficients as well as the maximum
c respective errors.
c
        erralm=0.
        errblm=0.
        write(*,2)
    2   format(5x,'k',4x,'erra',8x,'errb',10x,'erra',8x,'errb'/)
        do 20 in=1,ncap
          inm1=in-1
          erras=abs(als(in))
          errbs=abs((bes(in)-bexact(in))/bexact(in))
          erral=abs(all(in))
          errbl=abs((bel(in)-bexact(in))/bexact(in))
          if(erral.gt.erralm) erralm=erral
          if(errbl.gt.errblm) errblm=errbl
          if(ierrs.eq.0 .or. inm1.lt.abs(ierrs)) then
            if(in.eq.1) then
              write(*,3) inm1,erras,errbs,erral,errbl,ncap
    3         format(1x,i5,2e12.4,2x,2e12.4,'   N =',i4)
            else
              write(*,4) inm1,erras,errbs,erral,errbl
    4         format(1x,i5,2e12.4,2x,2e12.4)
            end if
          else
            if(in.eq.1) then
              write(*,5) inm1,erral,errbl,ncap
    5         format(1x,i5,26x,2e12.4,'   N =',i4)
            else
              write(*,6) inm1,erral,errbl
    6         format(1x,i5,26x,2e12.4)
            end if
          end if
   20   continue
        write(*,7) erralm,errblm
    7   format(/32x,2e12.4//)
   30 continue
      return
      end
      subroutine test04 ( )

c*********************************************************************72
c
cc TEST04 is a test of the routine  mcdis, which is applied to generate
c the first 80 recurrence coefficients for the orthogonal polynomials
c belonging to the weight function
c
c       (1-x**2)**(-1/2) + c   on (-1,1),  c = 1, 10, 100.
c
c The corresponding inner product is discretized by applying the
c Gauss-Chebyshev quadrature rule to the first term of the weight
c function, and the Gauss-Legendre rule to the second term. In
c addition to the beta-coefficients (all alpha's are zero), the
c routine prints the variables  ncap  and  kount  to confirm
c convergence after one iteration.
c
      external qchle_r4

      dimension a(81),b(81),e(81),xp(1),yp(1),endl(1),endr(1),
     *xfer(1),wfer(1),alpha(80),beta(80),be(80),x(81),w(81),
     *xm(162),wm(162)
      real p0(162)
      real p1(162)
      real p2(162)
      real betap(80,3)
      logical finl,finr

      common /common04/ c,a,b,e,epsma

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TEST04'
c
c epsma is the machine single precision.
c
      epsma=r1mach(3)
      iq=1
      idelta=2
      irout=1
      n=80
      mc=2
      mp=0
      ncapm=81
      eps=5000.*epsma
      iem=0
      c=.1

      do ic=1,3

        c=10.*c
c
c Compute the desired recursion coefficients. On machines with limited
c exponent range, harmless underflow may occur in the routine  gauss
c used in  QCHLE_R4  to generate the Gauss-Legendre quadrature rule.
c
        call mcdis(n,ncapm,mc,mp,xp,yp,qchle_r4,eps,iq,idelta,irout,finl,
     *    finr,endl,endr,xfer,wfer,alpha,beta,ncap,kount,ierr,ie,be,
     *    x,w,xm,wm,p0,p1,p2)

        write(*,2) ncap,kount,ierr,ie,c
    2   format(3x,'ncap = ',i3,' kount = ',i2,' ierr = ',i3,
     *    ' ie = ',i3,' for c = ',f5.1)

        iem = max ( iem, abs ( ie ) )

        if(ie.ne.0 .and. abs(ie).le.n) then
          call mcdis(abs(ie)-1,ncapm,mc,mp,xp,yp,qchle_r4,eps,iq,
     &    idelta, irout,finl,finr,endl,endr,xfer,wfer,alpha,beta,
     &    ncap,kount,ierr,ie,be,x,w,xm,wm,p0,p1,p2)
          write(*,2) ncap,kount,ierr,ie,c
          write ( * , '(a)' ) ' '
        end if
c
c Assemble the results in an array.
c
        do k=1,n
          if(ie.eq.0 .or. k-1.lt.abs(ie)) then
            betap(k,ic)=beta(k)
          else
            betap(k,ic)=0.
          end if
        end do

      end do
c
c Print the results.
c
      write(*,3)
    3 format(//3x,'k',4x,'beta(k), c=1',6x,'beta(k), c=10',4x,
     *  'beta(k), c=100'/)
      do 30 k=1,n
        km1=k-1
        if(iem.eq.0 .or. km1.lt.iem) then
          write(*,4) km1,betap(k,1),betap(k,2),betap(k,3)
        end if
    4   format(1x,i3,3e18.10)
   30 continue

      return
      end
      subroutine qchle_r4(n,x,w,i,ierr)

c*********************************************************************72
c
cc QCHLE_R4 returns Gauss-Chebyshev or Gauss-Legendre rules.
c
      dimension x(n),w(n),a(81),b(81),e(81)

      common /common04/ c,a,b,e,epsma

      fn=real(n)
      pi=4.*atan(1.)

      if(i.eq.1) then

        do 10 k=1,n
          fk=real(k)
          x(k)=cos((2.*fk-1.)*pi/(2.*fn))
          w(k)=pi/fn
   10   continue

      else

        call recur(n,1,0.,0.,a,b,ierr)
        call gauss(n,a,b,epsma,x,w,ierr,e)
        do 20 k=1,n
          w(k)=c*w(k)
   20   continue

      end if

      return
      end
      subroutine test05 ( )

c*********************************************************************72
c
      external qjac_r4
      dimension a(41),b(41),e(41),xp(1),yp(1),endl(1),endr(1),xfer(1),
     *wfer(1),alpha(40), beta(40),be(40),x(41),w(41),xm(42),wm(42),
     *p0(42),p1(42),p2(42)
      double precision dy,dalj,dbej,dalpbe,da(41),db(41),daex(40),
     *dbex(40),dnum,dd,dkm1,dden
      logical finl,finr

      common /common05/ a,b,e,epsma
c
c This test generates the first 40 recursion coefficients for
c polynomials orthogonal with respect to the Jacobi weight function
c with parameters  alj = -.8(.2)1., bej = -.8(.2)1.  and an added mass
c point of strength  y = .5, 1, 2, 4, 8  at the left end point. It
c also computes the maximum relative errors (absolute errors for alpha-
c coefficients near zero) of the computed coefficients by comparing
c them against the exact coefficients known analytically.
c
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TEST05'

      write(*,1)
    1 format(/)
c
c epsma is the machine single precision
c
      epsma=r1mach(3)
      iq=1
      idelta=2
      irout=1
      n=40
      mc=1
      mp=1
      xp(1)=-1.
      ncapm=41
      eps=5000.*epsma
      y=.25
      do 40 iy=1,5
        y=2.*y
        dy=dble(y)
        yp(1)=y
        write(*,2) y
    2   format(/1x,'y = ',f6.2/)
        write(*,3)
    3   format(2x,'alj',2x,'bej',5x,'erra',7x,'errb',6x,'alpha',
     *    4x,'beta',4x,'ka',2x,'kb',1x,'ierr',1x,'ie',1x,'it'/)
        do 30 ia=1,10
          alj=-1.+.2*real(ia)
          dalj=dble(alj)
          do 20 ib=1,10
            bej=-1.+.2*real(ib)
            alpbe=alj+bej
            dbej=dble(bej)
            dalpbe=dalj+dbej
c
c Generate the Jacobi recurrence coefficients.
c
            call recur(ncapm,6,alj,bej,a,b,ierr)
            call drecur(ncapm,6,dalj,dbej,da,db,iderr)
c
c Compute the desired recursion coefficients.
c
            call mcdis(n,ncapm,mc,mp,xp,yp,qjac_r4,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
c Compute the exact coefficients by Eqs. (4.19)-(4.21) of the companion
c paper along with the relative errors (absolute errors for alpha-
c coefficients close to zero).
c
            daex(1)=(da(1)-dy)/(1.d0+dy)
            dbex(1)=1.d0+dy
            erra=abs(sngl(dble(alpha(1))-daex(1)))
            if(abs(sngl(daex(1))).gt.eps) erra=erra/abs(sngl(daex(1)))
            errb=abs(sngl((dble(beta(1))-dbex(1))/dbex(1)))
            erram=erra
            alpham=alpha(1)
            errbm=errb
            betam=beta(1)
            kam=0
            kbm=0
            dnum=1.d0+dy
            dd=1.d0
            do 10 k=2,n
              km1=k-1
              dkm1=dble(km1)
              dden=dnum
              if(k.gt.2) dd=(dbej+dkm1)*(dalpbe+dkm1)*dd/((dalj+dkm1
     *          -1.d0)*(dkm1-1.d0))
              dnum=(1.d0+(dbej+dkm1+1.d0)*(dalpbe+dkm1+1.d0)*dy*dd/
     *          (dkm1*(dalj+dkm1)))/(1.d0+dy*dd)
              daex(k)=da(k)+2.d0*dkm1*(dkm1+dalj)*(dnum-1.d0)/
     *          ((dalpbe+2.d0*dkm1)*(dalpbe+2.d0*dkm1+1.d0))
     *          +2.d0*(dbej+dkm1+1.d0)*(dalpbe+dkm1+1.d0)*((1.d0/dnum)
     *          -1.d0)/((dalpbe+2.d0*dkm1+1.d0)*(dalpbe+2.d0*dkm1+2.d0))
              dbex(k)=dnum*db(k)/dden
              erra=abs(sngl(dble(alpha(k))-daex(k)))
              if(abs(sngl(daex(k))).gt.eps) erra=erra/abs(sngl(daex(k)))
              errb=abs(sngl((dble(beta(k))-dbex(k))/dbex(k)))
              if(erra.gt.erram) then
                erram=erra
                alpham=alpha(k)
                kam=km1
              end if
              if(errb.gt.errbm) then
                errbm=errb
                betam=beta(k)
                kbm=km1
              end if
   10       continue
c
c Print the results.
c
            write(*,4) alj,bej,erram,errbm,alpham,betam,kam,kbm,ierr,
     *        ie,kount
    4       format(1x,2f5.2,2e11.4,2f9.6,4i4,i2)
   20     continue
          write(*,1)
   30   continue
   40 continue
      return
      end
      subroutine qjac_r4(n,x,w,i,ierr)

c*********************************************************************72
c
cc QJAC_R4 returns a Gauss-Jacobi rule.
c
      dimension x(n),w(n),a(41),b(41),e(41)
      common /common05/ a,b,e,epsma

      call gauss(n,a,b,epsma,x,w,ierr,e)
      do 10 k=1,n
        w(k)=w(k)/b(1)
   10 continue
      return
      end
      subroutine test06 ( )

c*********************************************************************72
c
      external qlag,dqlag
      dimension a(100),b(100),e(100),xp(1),yp(1),endl(1),endr(1),
     *xfer(1),wfer(1),alpha(40),beta(40),be(40),x(100),w(100),
     *xm(200),wm(200),p0(200),p1(200),p2(200)
      double precision d1mach,da(300),db(300),de(300),depsma,dxp(1),
     *dyp(1),dendl(1),dendr(1),dxfer(1),dwfer(1),dalpha(40),dbeta(40),
     *dbe(40),dx(300),dw(300),dxm(600),dwm(600),dp0(600),dp1(600),
     *dp2(600),deps
      logical finl,finr,finld,finrd
      common /common06/ a,b,e,epsma
      common/common06_2/da,db,de,depsma
c
c This test generates in single and double precision the first 40
c recursion coefficients of the orthogonal polynomials belonging to
c the logistic density function
c
c              exp(-x)/((1+exp(-x))**2)  on (-oo, oo).
c
c It prints the double-precision beta-coefficients (the alpha's being
c all zero) along with the absolute and relative errors of the alpha-
c coefficients resp. beta-coefficients.
c
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TEST06'

      write(*,1)
    1 format(/)
c
c epsma and depsma are the machine single and double precision.
c
      iq=1
      idelta=1
      irout=1
      epsma=r1mach(3)
      depsma=d1mach(3)
      n=40
      mc=2
      mp=0
      ncapm=100
      ncapmm=mc*ncapm+mp
      ncpmd=300
      ncpmmd=mc*ncpmd+mp
      eps=5000.*epsma
      deps=1000.*depsma
c
c Compute the desired coefficients. On machines with limited exponent
c range, some of the weights in the Gauss-Laguerre quadrature rule may
c underflow.
c
      call mcdis(n,ncapm,mc,mp,xp,yp,qlag,eps,iq,idelta,irout,finl,
     *  finr,endl,endr,xfer,wfer,alpha,beta,ncap,kount,ierr,ie,be,x,w,
     *  xm,wm,p0,p1,p2)
      call dmcdis(n,ncpmd,mc,mp,dxp,dyp,dqlag,deps,iq,idelta,irout,
     *  finld,finrd,dendl,dendr,dxfer,dwfer,dalpha,dbeta,ncapd,kountd,
     *  ierrd,ied,dbe,dx,dw,dxm,dwm,dp0,dp1,dp2)
      write(*,2) ncap,kount,ierr,ie
    2 format(/1x,'ncap = ',i3,' kount = ',i2,' ierr = ',i3,' ie = ',i5)
      write(*,3) ncapd,kountd,ierrd,ied
    3 format(1x,'ncapd= ',i3,' kountd= ',i2,' ierrd= ',i3,' ied= ',i5/)
c
c Print the results.
c
      write(*,4)
    4 format(/5x,'k',13x,'dbeta(k)',17x,'erra',8x,'errb'/)
      do 10 k=1,n
        km1=k-1
        erra=abs(alpha(k))
        errb=abs(sngl((dble(beta(k))-dbeta(k))/dbeta(k)))
        if(ied.eq.0 .or. km1.lt.ied) then
          if(ie.eq.0 .or. km1.lt.ie) then
            write(*,5) km1,dbeta(k),erra,errb
    5       format(1x,i5,d33.25,2e12.4)
          else
            write(*,6) km1,dbeta(k)
    6       format(1x,i5,d33.25)
          end if
        end if
   10 continue
      return
      end
      subroutine qlag(n,x,w,i,ierr)

c*********************************************************************72
c
      dimension x(n),w(n),a(100),b(100),e(100)
      common/common06/a,b,e,epsma
      call recur(n,7,0.,0.,a,b,ierr)
      call gauss(n,a,b,epsma,x,w,ierr,e)
      do 10 k=1,n
        w(k)=w(k)/((1.+exp(-x(k)))**2)
        if(i.eq.1) x(k)=-x(k)
   10 continue
      return
      end
      subroutine dqlag(n,dx,dw,i,ierr)

c*********************************************************************72
c
      double precision dx(n),dw(n),da(300),db(300),de(300),depsma
      common/common06_2/da,db,de,depsma
      call drecur(n,7,0.d0,0.d0,da,db,ierr)
      call dgauss(n,da,db,depsma,dx,dw,ierr,de)
      do 10 k=1,n
        dw(k)=dw(k)/((1.d0+dexp(-dx(k)))**2)
        if(i.eq.1) dx(k)=-dx(k)
   10 continue
      return
      end
      subroutine test07 ( )

c*********************************************************************72
c
      external quad,dquad
      dimension endl(4),endr(4),xp(1),yp(1),xfer(100),wfer(100),
     *alpha(40),beta(40),be(40),x(100),w(100),xm(400),wm(400),p0(400),
     *p1(400),p2(400)
      double precision d1mach,depsma,dendl(4),dendr(4),dxp(1),dyp(1),
     *dxfer(250),dwfer(250),dalpha(40),dbeta(40),dbe(40),dx(250),
     *dw(250),dxm(1000),dwm(1000),dp0(1000),dp1(1000),dp2(1000),di,deps
      logical finl,finr,finld,finrd
c
c This test generates in single and double precision the first 40
c recurrence coefficients of the orthogonal polynomials for the half-
c range Hermite weight function
c
c                   exp(-x**2)  on  (0,oo).
c
c Printed are the double-precision values of the alpha- and beta-
c coefficients along with the respective relative errors of the single-
c precision values.
c
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TEST07'

      finl=.true.
      finr=.false.
      finld=.true.
      finrd=.false.
      epsma=r1mach(3)
      depsma=d1mach(3)
c
c epsma and depsma are the machine single and double precision.
c
      iq=2
      idelta=1
      irout=1
      n=40
      mc=4
      mcd=4
      mp=0
      ncapm=100
      ncpmd=250
c
c Set up the partition for the discretization of the inner product.
c
      do 10 i=1,4
        fi=real(i)
        di=dble(fi)
        endl(i)=3.*(fi-1.)
        endr(i)=3.*fi
        dendl(i)=3.d0*(di-1.d0)
        dendr(i)=3.d0*di
   10 continue
      eps=50.*epsma
      deps=1000.*depsma
c
c Compute the desired recursion coefficients by the multiple-component
c discretization procedure. On the third and fourth subinterval of the
c partition, the quadrature weights produced by  qgp  resp.  dqgp  may
c underflow, on the fourth subinterval even on machines with large
c exponent range.
c
      call 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)
      call dmcdis(n,ncpmd,mcd,mp,dxp,dyp,dquad,deps,iq,idelta,irout,
     *  finld,finrd,dendl,dendr,dxfer,dwfer,dalpha,dbeta,ncapd,kountd,
     *  ierrd,ied,dbe,dx,dw,dxm,dwm,dp0,dp1,dp2)
      write(*,1) ncap,kount,ierr,ie
    1 format(/1x,'ncap = ',i3,' kount = ',i2,' ierr = ',i3,' ie = ',i3)
      write(*,2) ncapd,kountd,ierrd,ied
    2 format(1x,'ncapd =',i3,' kountd =',i2,' ierrd =',i3,' ied =',i3/)
c
c Print the results.
c
      write(*,3)
    3 format(/5x,'k',13x,'dalpha(k)',25x,'dbeta(k)')
      write(*,4)
    4 format(11x,'erra',29x,'errb')
      do 20 k=1,n
        km1=k-1
        erra=abs(sngl((dble(alpha(k))-dalpha(k))/dalpha(k)))
        errb=abs(sngl((dble(beta(k))-dbeta(k))/dbeta(k)))
        write(*,5) km1,dalpha(k),dbeta(k)
    5   format(1x,i5,2d33.25)
        write(*,6) erra,errb
    6   format(6x,e12.4,21x,e12.4)
   20 continue
      return
      end
      subroutine quad(n,x,w,i,ierr)

c*********************************************************************72
c
      dimension x(n),w(n)
      print *,' User has selected the wrong SP-quadrature routine.'
      stop
      end
      subroutine dquad(n,dx,dw,i,ierr)

c*********************************************************************72
c
      double precision dx(n),dw(n)
      print *,' User has selected the wrong DP-quadrature routine.'
      stop
      end
      function wf(x,i)

c*********************************************************************72
c
      wf=exp(-x*x)
      return
      end
      double precision function dwf(dx,i)

c*********************************************************************72
c
      double precision dx
      dwf=dexp(-dx*dx)
      return
      end
      subroutine test08 ( )

c*********************************************************************72
c
      external qcheb
      dimension oom2(7),xp(1),yp(1),endl(1),endr(1),xfer(1),wfer(1),
     *a(79),b(79),fnu(80),alpha(40),beta(40),be(40),x(500),w(500),
     *xm(500),wm(500),s(40),s0(80),s1(80),s2(80)
      logical finl,finr
      common /common08/ om2
      data oom2/.1,.3,.5,.7,.9,.99,.999/
c
c This test reproduces the results of  test1  in single precision,
c for n=40, using the routine  mccheb  in place of CHEB_R4  and
c Gauss-Chebyshev quadrature to discretize the modified moments.
c
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TEST08'

      write(*,1)
    1 format(/)
      epsma=r1mach(3)
c
c epsma is the machine single precision.
c
      iq=1
      idelta=1
      n=40
      ndm1=2*n-1
      mc=1
      mp=0
      ncapm=500
      eps=100.*epsma
c
c Generate the recurrence coefficients for the (Chebyshev) polynomials
c defining the modified moments.
c
      call recur(ndm1,3,0.,0.,a,b,ierr)
c
c Compute the desired recursion coefficients by the discretized
c Chebyshev algorithm.
c
      do 20 iom=1,7
        om2=oom2(iom)
        call mccheb(n,ncapm,mc,mp,xp,yp,qcheb,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
c On machines with limited single-precision exponent range, the routine
c CHEB_R4 may have generated an underflow exception, which however is
c harmless and can be ignored.
c
        write(*,2) ncap,kount,ierr
    2   format(/'ncap=',i3,'  kount=',i3,'  ierr=',i3/)
c
c Print the results.
c
        write(*,3)
    3   format(5x,'k',6x,'beta(k)'/)
        do 10 k=1,n
          km1=k-1
          if(k.eq.1) then
            write(*,4) km1,beta(k),om2
    4       format(1x,i5,e18.10,'   om2 =',f6.3)
          else
            write(*,5) km1,beta(k)
    5     format(1x,i5,e18.10)
          end if
   10   continue
        write(*,1)
   20 continue
      return
      end
      subroutine qcheb(n,x,w,i,ierr)

c*********************************************************************72
c
      dimension x(n),w(n)
      common/common08/om2
      fn=real(n)
      pi=4.*atan(1.)
      do 10 k=1,n
        fk=real(k)
        x(k)=cos((2.*fk-1.)*pi/(2.*fn))
        w(k)=pi/(fn*sqrt(1.-om2*x(k)**2))
   10 continue
      return
      end
      subroutine test09 ( )

c*********************************************************************72
c
      dimension a(199),b(199),fnu(200),alpha(100),beta(100),s(100),
     *s0(200),s1(200),s2(200),alphc(100),betc(100)
      double precision dsigma,da(199),db(199),dnu(200),dalpha(100),
     *dbeta(100),ds(100),ds0(200),ds1(200),ds2(200),dalphc(100),
     *dbetc(100)
c
c This test recomputes the results of  test2  for  sigma=.5  by applying
c the routine  chri  with  iopt=1, x=0  to the weight function with
c parameter  sigma=-.5. Printed are the relative discrepancies in both
c single and double precision between these results and those obtained
c in  test 2  by the modified Chebyshev algorithm. The test is embedded
c in the routine  test2, from which all print statements have been
c removed.
c
      logical modmom,intexp

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TEST09'

      modmom=.true.
c
c Generate the recursion coefficients for the polynomials defining the
c modified resp. ordinary moments.
c
      if(modmom) then
        n=100
        ndm1=2*n-1
        call recur(ndm1,2,0.,0.,a,b,ierr)
        call drecur(ndm1,2,0.d0,0.d0,da,db,iderr)
      else
        n=12
        ndm1=2*n-1
        do 10 k=1,ndm1
          a(k)=0.
          b(k)=0.
          da(k)=0.d0
          db(k)=0.d0
   10   continue
      end if
      do 30 is=1,3
        dsigma=-.5d0+.5d0*dble(is-1)
        sigma=sngl(dsigma)
        if(is.eq.2) then
          intexp=.true.
        else
          intexp=.false.
        end if
c
c Compute the modified resp. ordinary moments using Eqs. (3.12) and
c (3.11) of the companion paper. On machines with limited exponent
c range, some of the high-order modified moments may underflow, without
c this having any deteriorating effect on the accuracy.
c
        call mm_09_r4(n,modmom,intexp,sigma,fnu)
        call mm_09_r8(n,modmom,intexp,dsigma,dnu)
c
c Compute the desired recursion coefficients by means of the modified
c Chebyshev algorithm; for the latter, see, e.g., Section 2.4 of
c W. Gautschi, On generating orthogonal polynomials'', SIAM J. Sci.
c Statist. Comput. 3, 1982, 289-317.
c
        call cheb_r4(n,a,b,fnu,alpha,beta,s,ierr,s0,s1,s2)
c
c On machines with limited single-precision exponent range, the routine
c CHEB_R4 may generate an underflow exception, which however is harmless
c and can be ignored.
c
        call cheb_r8(n,da,db,dnu,dalpha,dbeta,ds,iderr,ds0,ds1,ds2)
c
c Up to this point the code is identical with the one of  test2.
c
        if(is.eq.1) then
          write(*,1) ierr,iderr
    1     format(/1x,'ierr in cheb_r4 = ',i4,' iderr in cheb_r8 = ',i4/)
          if(ierr.ne.0) then
            nc=abs(ierr)
          else
            nc=n
          end if
          if(iderr.ne.0) then
            ncd=abs(iderr)
          else
            ncd=n
          end if
c
c Compute the desired recursion coefficients by a modification
c algorithm.
c
          nm1=nc-1
          call chri(nm1,1,alpha,beta,0.,0.,0.,0.,alphc,betc,ierr)
          nm1=ncd-1
          call chri_r8(nm1,1,dalpha,dbeta,0.d0,0.d0,0.d0,0.d0,dalphc,
     *      dbetc,iderr)
        end if
        if(is.eq.3) then
          write(*,2)
    2     format(/1x,'test of the results for sigma=1/2'/)
          np=nc
          if(ncd.lt.nc) np=ncd
          nm1=np-1
c
c Compute and print the relative discrepancies between the results of
c the modified Chebyshev algorithm and the modification algorithm.
c
          write(*,3)
    3     format(3x,'k',2x,'err alpha',4x,'err beta',12x,'err dalpha',
     *      2x,'err dbeta'/)
          do 15 k=1,nm1
            km1=k-1
            errac=abs(alpha(k)-alphc(k))/alpha(k)
            errbc=abs(beta(k)-betc(k))/beta(k)
            errdac=sngl(dabs(dalpha(k)-dalphc(k))/dalpha(k))
            errdbc=sngl(dabs(dbeta(k)-dbetc(k))/dbetc(k))
            write(*,4) km1,errac,errbc,errdac,errdbc
    4       format(1x,i3,2e12.4,9x,2e12.4)
   15     continue
          write(*,5)
    5     format(/1x,'end of test'/)
        end if
c
c The rest of the code is essentially the same as the corresponding
c piece of code in  test2  with all print statements removed.
c
        eamax=0.
        ebmax=0.
        do 20 k=1,n
          km1=k-1
          erra=sngl(dabs(dble(alpha(k))-dalpha(k))/dalpha(k))
          errb=sngl(dabs(dble(beta(k))-dbeta(k))/dbeta(k))
          if(erra.gt.eamax) then
            eamax=erra
            kamax=km1
          end if
          if(errb.gt.ebmax) then
            ebmax=errb
            kbmax=km1
          end if
   20   continue
   30 continue
      return
      end
      subroutine mm_09_r4(n,modmom,intexp,sigma,fnu)

c*********************************************************************72
c
c
cc MM_09_R4 generates the first  2*n  modified moments (if modmom=.true.)
c relative to shifted monic Legendre polynomials, using Eq. (3.12) of
c the companion paper, and the first  2*n  ordinary moments (if modmom
c =.false.) by Eq. (3.11), of the weight function
c
c          (x**sigma)*ln(1/x)  on (0,1],   sigma > -1,
c
c for sigma an integer (if intexp=.true.) or a real number (if intexp
c =.false.). In either case, the input variable sigma is of type real.
c
      dimension fnu(*)
      logical modmom,intexp
c
c The array  fnu  is assumed to have dimension  2*n.
c
      nd=2*n
      sigp1=sigma+1.
      if(modmom) then
        isigma=int(sigma)
        isigp1=isigma+1
        isigp2=isigma+2
        isigp3=isigma+3
        if(intexp .and. isigp1.lt.nd) then
          kmax=isigp1
        else
          kmax=nd
        end if
        c=1.
        do 20 k=1,kmax
          km1=k-1
          fk=real(k)
          p=1.
          s=1./sigp1
          if(kmax.gt.1) then
            do 10 i=1,km1
              fi=real(i)
              p=(sigp1-fi)*p/(sigp1+fi)
              s=s+1./(sigp1+fi)-1./(sigp1-fi)
   10       continue
          end if
          fnu(k)=c*s*p/sigp1
          c=fk*c/(4.*fk-2.)
   20   continue
        if(.not.intexp .or. isigp1.ge.nd) return
        q=-.5
        if(isigma.gt.0) then
          do 30 iq=1,isigma
            fiq=real(iq)
            q=fiq*fiq*q/((2.*fiq+1.)*(2.*fiq+2.))
   30     continue
        end if
        fnu(isigp2)=c*q
        if(isigp2.eq.nd) return
        do 40 k=isigp3,nd
          km1=k-1
          fkm1=real(km1)
          fnu(k)=-fkm1*(fkm1-sigp1)*fnu(km1)/((4.*fkm1-2.)*
     *      (fkm1+sigp1))
   40   continue
        return
      else
        do 50 k=1,nd
          fkm1=real(k-1)
          fnu(k)=(1./(sigp1+fkm1))**2
   50   continue
      end if
      end
      subroutine mm_09_r8(n,modmom,intexp,dsigma,dnu)

c*********************************************************************72
c
c
cc MM_09_R8 is a double-precision version of the routine MM_09_R4.
c
      double precision dsigma,dnu(*),dsigp1,dc,dk,dp,ds,di,dq,diq,dkm1
      logical modmom,intexp
c
c The array  dnu  is assumed to have dimension  2*n.
c
      nd=2*n
      dsigp1=dsigma+1.d0
      if(modmom) then
        isigma=idint(dsigma)
        isigp1=isigma+1
        isigp2=isigma+2
        isigp3=isigma+3
        if(intexp .and. isigp1.lt.nd) then
          kmax=isigp1
        else
          kmax=nd
        end if
        dc=1.d0
        do 20 k=1,kmax
          km1=k-1
          dk=dble(k)
          dp=1.d0
          ds=1.d0/dsigp1
          if(kmax.gt.1) then
            do 10 i=1,km1
              di=dble(i)
              dp=(dsigp1-di)*dp/(dsigp1+di)
              ds=ds+1.d0/(dsigp1+di)-1.d0/(dsigp1-di)
   10       continue
          end if
          dnu(k)=dc*ds*dp/dsigp1
          dc=dk*dc/(4.d0*dk-2.d0)
   20   continue
        if(.not.intexp .or. isigp1.ge.nd) return
        dq=-.5d0
        if(isigma.gt.0) then
          do 30 iq=1,isigma
            diq=dble(iq)
            dq=diq*diq*dq/((2.d0*diq+1.d0)*(2.d0*diq+2.d0))
   30     continue
        end if
        dnu(isigp2)=dc*dq
        if(isigp2.eq.nd) return
        do 40 k=isigp3,nd
          km1=k-1
          dkm1=dble(km1)
          dnu(k)=-dkm1*(dkm1-dsigp1)*dnu(km1)/((4.d0*dkm1-2.d0)*
     *      (dkm1+dsigp1))
   40   continue
        return
      else
        do 50 k=1,nd
          dkm1=dble(k-1)
          dnu(k)=(1.d0/(dsigp1+dkm1))**2
   50   continue
      end if
      end
      subroutine test10 ( )

c*********************************************************************72
c
      dimension a(31),b(31),alpha(31),beta(31),z(11),w(11),e(11),
     *a1(31),b1(31),betap(20,12),erram(12),errbm(12)
      double precision depsma,d1mach,da(31),db(31),dalpha(31),dbeta(31),
     *dz(11),dw(11),de(11),da1(31),db1(31)
c
c epsma and depsma are the machine single and double precision.
c

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TEST10'

      epsma=r1mach(3)
      depsma=d1mach(3)
c
c This test applies the routines  indp  and  dindp  to generate the
c first 20 recursion coefficients of the induced Legendre polynomials
c pind(k,m)(.), m=0,1,2,...,11, that is, of the polynomials orthogonal
c relative to the weight function
c
c                 [p(m)(x)]**2    on [-1,1],
c
c where  p(m)(.)  is the (monic) Legendre polynomial of degree m.
c (When m=0, then  pind(k,0)(.)=p(k)(.).) The routine also prints the
c absolute and relative errors, respectively, of the alpha- and beta-
c coefficients.
c
      n=20
      do 20 im=1,12
        m=im-1
        npm=n+m
c
c Generate the Legendre recurrence coefficients required in the
c routines  indp  and dindp.
c
        call recur(npm,1,0.,0.,a,b,ierr)
        call drecur(npm,1,0.d0,0.d0,da,db,ierr)
c
c Compute the desired recursion coefficients.
c
        call indp(n,m,a,b,epsma,alpha,beta,ierr,z,w,e,a1,b1)
        call dindp(n,m,da,db,depsma,dalpha,dbeta,ierr,dz,dw,
     *    de,da1,db1)
c
c Compute and print the respective errors.
c
        erram(im)=0.
        errbm(im)=0.
        do 10 k=1,n
          erra=sngl(dabs(dble(alpha(k))-dalpha(k)))
          errb=sngl(dabs((dble(beta(k))-dbeta(k))/dbeta(k)))
          if(erra.gt.erram(im)) erram(im)=erra
          if(errb.gt.errbm(im)) errbm(im)=errb
          betap(k,im)=sngl(dbeta(k))
   10   continue
   20 continue
      do 40 ip=1,3
        ip4=1+4*(ip-1)
        write(*,1) ip4-1,ip4,ip4+1,ip4+2
    1   format(5x,'k',2x,'m=',i1,'  beta(k)',2x,'m=',i1,'  beta(k)',2x,
     *    'm=',i2,' beta(k)',2x,'m=',i2,' beta(k)'/)
        do 30 k=1,n
          km1=k-1
          write(*,2) km1,betap(k,ip4),betap(k,ip4+1),betap(k,ip4+2),
     *      betap(k,ip4+3)
    2     format(1x,i5,4f14.10)
   30   continue
        write(*,3) erram(ip4),erram(ip4+1),erram(ip4+2),erram(ip4+3)
    3   format(/4x,'erra',e12.4,3e14.4)
        write(*,4) errbm(ip4),errbm(ip4+1),errbm(ip4+2),errbm(ip4+3)
    4   format(4x,'errb',e12.4,3e14.4//)
   40 continue
      return
      end
      subroutine indp(n,m,a,b,eps,alpha,beta,ierr,z,w,e,a1,b1)

c*********************************************************************72
c
c
c If  p(m)(.)  denotes the (monic) orthogonal polynomial of degree  m
c relative to the weight function  w(x), then the corresponding m-th
c induced orthogonal polynomials  pind(k,m)(.), k=0,1,2,..., are those
c orthogonal with respect to the weight function
c
c                   (p(m)(x)**2)*w(x).
c
c (For background on induced orthogonal polynomials, including an
c algorithm for generating their recursion coefficients, see W. Gautschi
c and S. Li,A set of orthogonal polynomials induced by a given
c orthogonal polynomial'', Aequationes Math., to appear.) This routine
c obtains the first n recurrence coefficients of the m-th induced
c orthogonal polynomials by an m-fold application of the routine  chri
c with  iopt=7, the shifts taken being, in succession, the zeros of
c p(m)(.).
c
      dimension a(*),b(*),alpha(*),beta(*),z(m),w(m),e(m),
     *a1(*),b1(*)
c
c The arrays  a,b,alpha,beta,a1,b1  are assumed to have dimension  n+m.
c
      npm=n+m
      do 10 k=1,npm
        alpha(k)=a(k)
        beta(k)=b(k)
   10 continue
      if(m.eq.0) return
      call gauss(m,a,b,eps,z,w,ierr,e)
      do 30 imu=1,m
        mi=npm-imu
        do 20 k=1,mi+1
          a1(k)=alpha(k)
          b1(k)=beta(k)
   20   continue
        x=z(imu)
        call chri(mi,7,a1,b1,x,0.,0.,0.,alpha,beta,ierrc)
   30 continue
      return
      end
      subroutine dindp(n,m,da,db,deps,dalpha,dbeta,ierr,dz,dw,de,
     *da1,db1)

c*********************************************************************72
c
c
c This is a double-precision version of the routine  indp.
c
      double precision da(*),db(*),deps,dalpha(*),dbeta(*),
     *dz(m),dw(m),de(m),da1(*),db1(*),dx
c
c The arrays  da,db,dalpha,dbeta,da1,db1  are assumed to have
c dimension  n+m.
c
      npm=n+m
      do 10 k=1,npm
        dalpha(k)=da(k)
        dbeta(k)=db(k)
   10 continue
      if(m.eq.0) return
      call dgauss(m,da,db,deps,dz,dw,ierr,de)
      do 30 imu=1,m
        mi=npm-imu
        do 20 k=1,mi+1
          da1(k)=dalpha(k)
          db1(k)=dbeta(k)
   20   continue
        dx=dz(imu)
        call chri_r8(mi,7,da1,db1,dx,0.d0,0.d0,0.d0,dalpha,
     *    dbeta,ierrc)
   30 continue
      return
      end
      subroutine test11 ( )

c*********************************************************************72
c
      complex rho,rold,z,e
      dimension xx(5),rr(5),a(500),b(500),alpha(40),beta(40),alphc(40),
     *betc(40),fnu(80),rho(80),rold(80),s(40),s0(80),s1(80),s2(80),
     *alphr(40),betr(40),alphcr(40),betcr(40)
      double precision d1mach,depsma,deps,dal,dbe,dhi,da(800),db(800),
     *dx,dy,dalpha(40),dbeta(40),dnu(80),drhor(80),drhoi(80),
     *droldr(80),droldi(80),ds(40),ds0(80),ds1(80),ds2(80),dhr,
     *dalphc(40),dbetc(40),dalphr(40),dbetr(40),dalcr(40),dbecr(40)
      data xx/1.001,1.01,1.04,1.07,1.1/
      data rr/1.05,1.1625,1.275,1.3875,1.5/
c
c This test is to illustrate the dissimilar performance of the routines
c chri  and  gchri  in the case of division of the Jacobi weight
c function  w(t;alj,bej)  with parameters  alj,bej  by either a linear
c divisor  t-x  or a quadratic divisor  (t-x)**2 + y**2 . In either
c case, the parameters selected are  alj=-.8(.4).8, bej=alj(.4).8.  In
c the former case, x = -1.001, -1.01, -1.04, -1.07 and -1.1, whereas in
c the latter case,  x and  y  are chosen to lie, regularly spaced, on
c the upper half of an ellipse with foci at +1 and -1 and sum of the
c semiaxes equal to  rho = 1.05, 1.1625, 1.275, 1.3875 and 1.5. The
c routines are run in both single and double precision with n=40, the
c results of the latter being used to calculate, and print, the maximum
c absolute and relative error of the single-precision alpha- and beta-
c coefficients, respectively. Also printed are the starting recurrence
c indexes required in the backward recurrence schemes of  gchri,dgchri
c to achieve single- resp. double-precision accuracy. This information
c is contained in the first line of each 3-line block of the output,
c where in the case of quadratic divisors only average values (averaged
c over the upper half of the respective ellipse) are shown. The second
c and third line of each 3-line block display the maximum
c reconstruction error'', that is, the maximum errors in the alpha's
c and beta's if the coefficients produced by  gchri,chri  and  dgchri,
c chri_r8  are fed back to the routines  chri  and  chri_r8  with  iopt=1
c to recover the original recursion coefficients in single and double
c precision.
c
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TEST11'

      write(*,1)
    1 format(/)
      epsma=r1mach(3)
      depsma=d1mach(3)
c
c epsma and depsma are the machine single and double precision.
c
      n=40
      np1=n+1
      nm1=n-1
      nd=2*n
      ndm1=nd-1
      numax=500
      numaxd=800
      eps=10.*epsma
      deps=100.d0*depsma
      epsd=sngl(deps)
      ipoly=6
      do 70 ial=1,5
        al=-.8+.4*real(ial-1)
        dal=dble(al)
        ibemax=6-ial
        do 60 ibe=1,ibemax
          be=al+.4*real(ibe-1)
          dbe=dble(be)
          write(*,2) al,be
    2     format(///1x,'al = ',f6.2,'  be = ',f6.2//)
          hi=0.
          dhi=0.d0
c
c Generate the Jacobi recurrence coefficients to be used in the
c backward recurrence algorithm of the routines  gchri  and  dgchri.
c
          call recur(numax,ipoly,al,be,a,b,ierr)
          call drecur(numaxd,ipoly,dal,dbe,da,db,ierrd)
          write(*,3)
    3     format(30x,'gchri',20x,'chri')
          write(*,4)
    4     format(5x,'x',5x,'nu0',2x,'nud0',4x,'erra',8x,'errb',10x,
     *      'erra',7x,'errb'/)
          do 20 ix=1,5
            x=-xx(ix)
            dx=dble(x)
            y=0.
            dy=0.d0
            z=cmplx(x,y)
c
c Compute the starting index for backward recurrence.
c
            nu0=nu0jac(ndm1,z,eps)
            nu0d=nu0jac(ndm1,z,epsd)
c
c Generate the recurrence coefficients for the Jacobi weight function
c divided by a linear divisor, using the routines  gchri,dgchri.
c
            call gchri(n,1,nu0,numax,eps,a,b,x,y,alpha,beta,nu,ierrg,
     *        ierrc,fnu,rho,rold,s,s0,s1,s2)
c
c On machines with limited single-precision exponent range, the routine
c CHEB_R4 used in  gchri  may have generated an underflow exception,
c which however is harmless and can be ignored.
c
            call dgchri(n,1,nu0d,numaxd,deps,da,db,dx,dy,dalpha,dbeta,
     *        nud,ierrgd,ierrcd,dnu,drhor,drhoi,droldr,droldi,ds,ds0,
     *        ds1,ds2)
            if(ierrg.ne.0 .or. ierrc.ne.0 .or. ierrgd.ne.0 .or.ierrcd
     *        .ne.0) then
              write(*,5) ierrg,ierrgd,al,be,x
    5         format(/1x,'ierrg in gchri = ',i4,' ierrg in dgchri = ',
     *          i4,' for al = ',f6.2,' be = ',f6.2,' x = ',f7.4)
              write(*,6) ierrc,ierrcd,al,be,x
    6         format(1x,'ierrc in gchri = ',i4,' ierrc in dgchri = ',
     *          i4,' for al = ',f6.2,' be = ',f6.2,' x = ',f7.4/)
              goto 20
            end if
c
c Generate the recurrence coefficients for the Jacobi weight function
c divided by a linear divisor, using the routines  chri, chri_r8.
c
            hr=real(rho(1))
            dhr=drhor(1)
            call chri(n,4,a,b,x,y,hr,hi,alphc,betc,ierr)
            call chri_r8(n,4,da,db,dx,dy,dhr,dhi,dalphc,dbetc,ierr)
c
c Do the reconstruction.
c
            call chri(nm1,1,alpha,beta,x,y,0.,0.,alphr,betr,ierr)
            call chri_r8(nm1,1,dalpha,dbeta,dx,dy,0.d0,0.d0,dalphr,
     &        dbetr, ierr)
            call chri(nm1,1,alphc,betc,x,y,0.,0.,alphcr,betcr,ierr)
            call chri_r8(nm1,1,dalphc,dbetc,dx,dy,0.d0,0.d0,dalcr,dbecr,
     *        ierr)
c
c Compute and print the maximum errors.
c
            erragm=0.
            errbgm=0.
            erracm=0.
            errbcm=0.
            erram=0.
            errbm=0.
            errdam=0.
            errdbm=0.
            eracrm=0.
            erbcrm=0.
            edacrm=0.
            edbcrm=0.
            do 10 k=1,n
              km1=k-1
              errag=abs(sngl(dble(alpha(k))-dalpha(k)))
              errbg=abs(sngl((dble(beta(k))-dbeta(k))/dbeta(k)))
              errac=abs(sngl(dble(alphc(k))-dalphc(k)))
              errbc=abs(sngl((dble(betc(k))-dbetc(k))/dbetc(k)))
              if(k.lt.n) then
                erra=abs(alphr(k)-a(k))
                errb=abs((betr(k)-b(k))/b(k))
                errda=abs(sngl(dalphr(k)-da(k)))
                errdb=abs(sngl((dbetr(k)-db(k))/db(k)))
                eracr=abs(alphcr(k)-a(k))
                erbcr=abs((betcr(k)-b(k))/b(k))
                edacr=abs(sngl(dalcr(k)-da(k)))
                edbcr=abs(sngl((dbecr(k)-db(k))/db(k)))
                if(erra.gt.erram) erram=erra
                if(errb.gt.errbm) errbm=errb
                if(errda.gt.errdam) errdam=errda
                if(errdb.gt.errdbm) errdbm=errdb
                if(eracr.gt.eracrm) eracrm=eracr
                if(erbcr.gt.erbcrm) erbcrm=erbcr
                if(edacr.gt.edacrm) edacrm=edacr
                if(edbcr.gt.edbcrm) edbcrm=edbcr
              end if
              if(errag.gt.erragm) erragm=errag
              if(errbg.gt.errbgm) errbgm=errbg
              if(errac.gt.erracm) erracm=errac
              if(errbc.gt.errbcm) errbcm=errbc
   10       continue
            write(*,7) x,nu0,nu0d,erragm,errbgm,erracm,errbcm
    7       format(/1x,f7.4,2i6,2e12.4,2x,2e12.4)
            if(ix.eq.1) then
              write(*,8) erram,errbm,errdam,errdbm
    8         format(11x,'reconstr.',2e12.4,2x,2e12.4)
              write(*,9) eracrm,erbcrm,edacrm,edbcrm
    9         format(12x,'errors',2x,2e12.4,2x,2e12.4)
            else
              write(*,11) erram,errbm,errdam,errdbm
   11         format(20x,2e12.4,2x,2e12.4)
              write(*,11) eracrm,erbcrm,edacrm,edbcrm
            end if
   20     continue
          write(*,1)
          ndiv=20
          ndivm1=ndiv-1
          fndiv=real(ndiv)
          fndm1=real(ndivm1)
          pi=4.*atan(1.)
          write(*,3)
          write(*,12)
   12     format(4x,'rho',4x,'nu0',2x,'nud0',4x,'erra',8x,'errb',10x,
     *      'erra',7x,'errb'/)
          do 50 ir=1,5
            r=rr(ir)
            agmv=0.
            bgmv=0.
            acmv=0.
            bcmv=0.
            amv=0.
            bmv=0.
            damv=0.
            dbmv=0.
            acrmv=0.
            bcrmv=0.
            dacrmv=0.
            dbcrmv=0.
            nu0v=0
            nu0dv=0
            do 40 ith=1,ndivm1
c
c Generate the points on the ellipse.
c
              theta=pi*real(ith)/fndiv
              e=cmplx(cos(theta),sin(theta))
              z=.5*(r*e+1./(r*e))
              x=real(z)
              y=aimag(z)
              dx=dble(x)
              dy=dble(y)
c
c Compute the starting index for backward recurrence.
c
              nu0=nu0jac(ndm1,z,eps)
              nu0d=nu0jac(ndm1,z,epsd)
c
c Generate the recurrence coefficients for the Jacobi weight function
c divided by a quadratic divisor, using the routines  gchri,dgchri.
c
              call gchri(n,2,nu0,numax,eps,a,b,x,y,alpha,beta,nu,ierrg,
     *          ierrc,fnu,rho,rold,s,s0,s1,s2)
              call dgchri(n,2,nu0d,numaxd,deps,da,db,dx,dy,dalpha,dbeta,
     *          nud,ierrgd,ierrcd,dnu,drhor,drhoi,droldr,droldi,ds,ds0,
     *          ds1,ds2)
              if(ierrg.ne.0 .or.ierrc.ne.0 .or. ierrgd.ne.0 .or. ierrcd
     *          .ne.0) then
                write(*,5) ierrg,ierrgd,al,be,x
                write(*,6) ierrc,ierrcd,al,be,x
                goto 40
              end if
              nu0v=nu0v+nu0
              nu0dv=nu0dv+nu0d
c
c Generate the recurrence coefficients for the Jacobi weight function
c divided by a quadratic divisor, using the routines  chri, chri_r8.
c
              hr=real(rho(1))
              hi=aimag(rho(1))
              dhr=drhor(1)
              dhi=drhoi(1)
              call chri(n,5,a,b,x,y,hr,hi,alphc,betc,ierr)
              call chri_r8(n,5,da,db,dx,dy,dhr,dhi,dalphc,dbetc,ierr)
c
c Do the reconstruction.
c
              call chri(nm1,2,alpha,beta,x,y,0.,0.,alphr,betr,ierr)
              call chri_r8(nm1,2,dalpha,dbeta,dx,dy,0.d0,0.d0,dalphr,
     *          dbetr,ierr)
              call chri(nm1,2,alphc,betc,x,y,0.,0.,alphcr,betcr,ierr)
              call chri_r8(nm1,2,dalphc,dbetc,dx,dy,0.d0,0.d0,dalcr,
     *          dbecr,ierr)
c
c Compute and print the maximum average errors.
c
              erragm=0.
              errbgm=0.
              erracm=0.
              errbcm=0.
              erram=0.
              errbm=0.
              errdam=0.
              errdbm=0.
              eracrm=0.
              erbcrm=0.
              edacrm=0.
              edbcrm=0.
              do 30 k=1,n
                km1=k-1
                errag=abs(sngl(dble(alpha(k))-dalpha(k)))
                errbg=abs(sngl((dble(beta(k))-dbeta(k))/dbeta(k)))
                errac=abs(sngl(dble(alphc(k))-dalphc(k)))
                errbc=abs(sngl((dble(betc(k))-dbetc(k))/dbetc(k)))
                if(k.lt.n) then
                  erra=abs(alphr(k)-a(k))
                  errb=abs((betr(k)-b(k))/b(k))
                  errda=abs(sngl(dalphr(k)-da(k)))
                  errdb=abs(sngl((dbetr(k)-db(k))/db(k)))
                  eracr=abs(alphcr(k)-a(k))
                  erbcr=abs((betcr(k)-b(k))/b(k))
                  edacr=abs(sngl(dalcr(k)-da(k)))
                  edbcr=abs(sngl((dbecr(k)-db(k))/db(k)))
                  if(erra.gt.erram) erram=erra
                  if(errb.gt.errbm) errbm=errb
                  if(errda.gt.errdam) errdam=errda
                  if(errdb.gt.errdbm) errdbm=errdb
                  if(eracr.gt.eracrm) eracrm=eracr
                  if(erbcr.gt.erbcrm) erbcrm=erbcr
                  if(edacr.gt.edacrm) edacrm=edacr
                  if(edbcr.gt.edbcrm) edbcrm=edbcr
                end if
                if(errag.gt.erragm) erragm=errag
                if(errbg.gt.errbgm) errbgm=errbg
                if(errac.gt.erracm) erracm=errac
                if(errbc.gt.errbcm) errbcm=errbc
   30         continue
              agmv=agmv+erragm
              bgmv=bgmv+errbgm
              acmv=acmv+erracm
              bcmv=bcmv+errbcm
              amv=amv+erram
              bmv=bmv+errbm
              damv=damv+errdam
              dbmv=dbmv+errdbm
              acrmv=acrmv+eracrm
              bcrmv=bcrmv+erbcrm
              dacrmv=dacrmv+edacrm
              dbcrmv=dbcrmv+edbcrm
   40       continue
            nu0=real(nu0v)/fndm1
            nu0d=real(nu0dv)/fndm1
            erragm=agmv/fndm1
            errbgm=bgmv/fndm1
            erracm=acmv/fndm1
            errbcm=bcmv/fndm1
            erram=amv/fndm1
            errbm=bmv/fndm1
            errdam=damv/fndm1
            errdbm=dbmv/fndm1
            eracrm=acrmv/fndm1
            erbcrm=bcrmv/fndm1
            edacrm=dacrmv/fndm1
            edbcrm=dbcrmv/fndm1
            if(ir.eq.1) then
              write(*,7) r,nu0,nu0d,erragm,errbgm,erracm,errbcm
              write(*,8) erram,errbm,errdam,errdbm
              write(*,9) eracrm,erbcrm,edacrm,edbcrm
            else
              write(*,7) r,nu0,nu0d,erragm,errbgm,erracm,errbcm
              write(*,11) erram,errbm,errdam,errdbm
              write(*,11) eracrm,erbcrm,edacrm,edbcrm
            end if
   50     continue
   60   continue
   70 continue

      return
      end
