!
!         this routine is a supplement to the journal article,
!      Algorithm 726: ORTHPOL - A Package of Routines for Generating
!        Orthogonal Polynomials and Gauss-Type Quadrature Rules.
!                         Walter Gautschi,
!      ACM Transactions on Mathematical Software, 20, 21-62(1994)
!
!            A footnote in Gautschi's article reveals that
!            this particular routine came from Gene Golub.
!
! Given  n  and a measure  dlambda, this routine generates the n-point
! Gaussian quadrature formula
!
!     integral over supp(dlambda) of f(x)dlambda(x)
!
!        = sum from k=1 to k=n of w(k)f(x(k)) + R(n;f).
!
! The nodes are returned as  zero(k)=x(k) and the weights as
! weight(k)=w(k), k=1,2,...,n. The user has to supply the recursion
! coefficients  alpha(k), beta(k), k=0,1,2,...,n-1, for the measure
! dlambda. The routine computes the nodes as eigenvalues, and the
! weights in term of the first component of the respective normalized
! eigenvectors of the n-th order Jacobi matrix associated with  dlambda.
! It uses a translation and adaptation of the algol procedure  imtql2,
! Numer. Math. 12, 1968, 377-383, by Martin and Wilkinson, as modified
! by Dubrulle, Numer. Math. 15, 1970, 450. See also Handbook for
! Autom. Comput., vol. 2 - Linear Algebra, pp.241-248, and the eispack
! routine  imtql2.
!
!        Input:  n - - the number of points in the Gaussian quadrature
!                      formula; type integer
!                alpha,beta - - arrays of dimension  n  to be filled
!                      with the values of  alpha(k-1), beta(k-1), k=1,2,
!                      ...,n
!                eps - the relative accuracy desired in the nodes
!                      and weights
!
!        Output: zero- array of dimension  n  containing the Gaussian
!                      nodes (in increasing order)  zero(k)=x(k), k=1,2,
!                      ...,n
!                weight - array of dimension  n  containing the
!                      Gaussian weights  weight(k)=w(k), k=1,2,...,n
!                ierr- an error flag equal to  0  on normal return,
!                      equal to  i  if the QR algorithm does not
!                      converge within 30 iterations on evaluating the
!                      i-th eigenvalue, equal to  -1  if  n  is not in
!                      range, and equal to  -2  if one of the beta's is
!                      negative.
!
! The array  wrk  is needed for working space.
!
      subroutine gauss(n,alpha,beta,eps,zero,weight,ierr,wrk)
      implicit none
      real*8 alpha(n),beta(n),zero(n),weight(n),wrk(n)
      real*8 dg, dr, dp, ds, dc, db, df
      integer*8 ierr, i, ii, j, k, l, m, mml, n
      if(n.lt.1) then
        ierr=-1
        return
      end if
!
      ierr=0
      zero(1)=alpha(1)
      if(beta(1).lt.0.0D+00) then
        ierr=-2
        return
      end if
      weight(1)=beta(1)
      if (n.eq.1) return
!
      weight(1)=1.0D+00
      wrk(n)=0.0D+00
      do 100 k=2,n
        zero(k)=alpha(k)
        if(beta(k).lt.0.0D+00) then
          ierr=-2
          return
        end if
        wrk(k-1)=sqrt(beta(k))
        weight(k)=0.0D+00
  100 continue
!
      do 240 l=1,n
        j=0
!
! Look for a small subdiagonal element.
!
  105   do 110 m=l,n
          if(m.eq.n) go to 120
          if(abs(wrk(m)).le.eps*(abs(zero(m))+abs(zero(m+1)))) go to 120
  110   continue
  120   dp=zero(l)
        if(m.eq.l) go to 240
        if(j.eq.30) go to 400
        j=j+1
!
! Form shift.
!
        dg=(zero(l+1)-dp)/(2.0D+00*wrk(l))
        dr=sqrt(dg*dg+1.0D+00)
        dg=zero(m)-dp+wrk(l)/(dg+sign(dr,dg))
        ds=1.0D+00
        dc=1.0D+00
        dp=0.0D+00
        mml=m-l

!
! For i=m-1 step -1 until l do ...
!
        do 200 ii=1,mml
          i=m-ii
          df=ds*wrk(i)
          db=dc*wrk(i)
          if(abs(df).lt.abs(dg)) go to 150
          dc=dg/df
          dr=sqrt(dc*dc+1.0D+00)
          wrk(i+1)=df*dr
          ds=1.0D+00/dr
          dc=dc*ds
          go to 160
  150     ds=df/dg
          dr=sqrt(ds*ds+1.0D+00)
          wrk(i+1)=dg*dr
          dc=1.0D+00/dr
          ds=ds*dc
  160     dg=zero(i+1)-dp
          dr=(zero(i)-dg)*ds+2.0D+00*dc*db
          dp=ds*dr
          zero(i+1)=dg+dp
          dg=dc*dr-db
!
! Form first component of vector.
!
          df=weight(i+1)
          weight(i+1)=ds*weight(i)+dc*df
          weight(i)=dc*weight(i)-ds*df
  200   continue
        zero(l)=zero(l)-dp
        wrk(l)=dg
        wrk(m)=0.0D+00
        go to 105
  240 continue
!
! Order eigenvalues and eigenvectors.
!
      do 300 ii=2,n
        i=ii-1
        k=i
        dp=zero(i)
        do 260 j=ii,n
          if(zero(j).ge.dp) go to 260
          k=j
          dp=zero(j)
  260   continue
        if(k.eq.i) go to 300
        zero(k)=zero(i)
        zero(i)=dp
        dp=weight(i)
        weight(i)=weight(k)
        weight(k)=dp

  300 continue
      do 310 k=1,n
        weight(k)=beta(1)*weight(k)*weight(k)
  310 continue
      return
!
! Set error - no convergence to an eigenvalue after 30 iterations.
!
  400 ierr=l
      return
      end