c     QGAUSH
c
c     Reference: Press, et. al.
c
c     Ellen McGrattan, 9-3-94
c
      subroutine qgaush(n,x,w)
      implicit double precision (a-h,o-z)
      parameter(m=50)
      double precision x(n),w(n),p(m),q(m),r(m),d(m),a(m,m),work(10*m)
c
c     Issue warning if m<n+1.
c
      if (m.lt.n+1) 
     c  write(*,*) 'WARNING: need higher value for m in QGAUSH'
c
c     Construct Hermite polynomials recursively: H(i+1)=2*x*H(i)-2*i*H(i-1).
c
      pi   = 3.14159265358979
      p(1) = 1.d0  
      q(1) = 2.d0  
      q(2) = 0.d0  
      do 10 j=1,n-1
        r(1) = 2.d0*q(1)
        r(2) = 2.d0*q(2)
        fj   = float(j)
        do 20 i=3,j+1
          r(i) = 2.d0*(q(i)-fj*p(i-2))
 20     continue
        r(j+2) = -2.d0*float(j)*p(j)
        do 30 i=1,j+2
          p(i) = q(i)
          q(i) = r(i)
 30     continue
 10   continue
c
c     Calculate the roots of H(n)=r(1)*x^n+r(2)*x^(n-1)+ ... + r(n+1).
c
      do 40 j=2,n
        do 50 i=1,n-1
          a(i,j) = 0.d0 
 50     continue
        a(j-1,j) = 1.d0
 40   continue
      a(n,1) = -r(n+1)/r(1)
      do 60 i=2,n  
        a(n-i+1,1) = 0.d0
        a(n,i)     = -r(n-i+2)/r(1)
 60   continue
      call dgeev('N','N',n,a,m,x,p,dum,1,dum,1,work,10*m,info)
c
c     Calculate the derivative of H(n).
c 
      do 70 i=1,n
        d(i) = float(n-i+1)*r(i)
        p(i) = 0.d0
        q(i) = 0.d0
        r(i) = 0.d0
 70   continue
c
c     Calculate the numerator for the weights recursively.
c 
      p(1) = 2.d0  
      q(1) = 4.d0  
      q(2) = 0.d0  
      do 80 j=2,n-1
        r(1) = 2.d0*q(1)
        r(2) = 2.d0*q(2)
        fj   = float(j)
        do 90 i=3,j+1
          r(i) = 2.d0*(q(i)-fj*p(i-2))
 90     continue
        r(j+2) = -2.d0*float(j)*p(j)
        do 100 i=1,j+2
          p(i) = q(i)
          q(i) = r(i)
 100    continue
 80   continue
c
c     Calculate weights.
c
      do 110 j=1,n
        prod = x(j)
        sum1 = 0.d0
        sum2 = 0.d0
        do 120 i=1,n-1
          sum1 = sum1 + r(n-i)*prod
          sum2 = sum2 + d(n-i)*prod
          prod = prod*x(j)
 120    continue
        sum1 = sum1+r(n)
        sum2 = sum2+d(n)
        w(j) = dsqrt(pi)*sum1/sum2
 110  continue
      return
      end