C  http://www.netlib.org/go/gaussq.f on 7 August 2012.
C  The function was modified for portability (Aug and Sep 2012) and
C  updated to Fortran 77 (28 Aug 2016) by Gordon Smyth.
C  All modified lines are commented out with a capital "C" and the new
C  version follows immediately.

subroutine gausq2(n, d, e, z, ierr)
c
c     this subroutine is a translation of an algol procedure,
c     num. math. 12, 377-383(1968) by martin and wilkinson,
c     as modified in num. math. 15, 450(1970) by dubrulle.
c     handbook for auto. comp., vol.ii-linear algebra, 241-248(1971).
c     this is a modified version of the 'eispack' routine imtql2.
c
c     this subroutine finds the eigenvalues and first components of the
c     eigenvectors of a symmetric tridiagonal matrix by the implicit ql
c     method.
c
c     on input:
c
c        n is the order of the matrix;
c
c        d contains the diagonal elements of the input matrix;
c
c        e contains the subdiagonal elements of the input matrix
c          in its first n-1 positions.  e(n) is arbitrary;
c
c        z contains the first row of the identity matrix.
c
c      on output:
c
c        d contains the eigenvalues in ascending order.  if an
c          error exit is made, the eigenvalues are correct but
c          unordered for indices 1, 2, ..., ierr-1;
c
c        e has been destroyed;
c
c        z contains the first components of the orthonormal eigenvectors
c          of the symmetric tridiagonal matrix.  if an error exit is
c          made, z contains the eigenvectors associated with the stored
c          eigenvalues;
c
c        ierr is set to
c          zero       for normal return,
c          j          if the j-th eigenvalue has not been
c                     determined after 30 iterations.
c
c     ------------------------------------------------------------------
c
integer i, j, k, l, m, n, ii, mml, ierr
C      real*8 d(n), e(n), z(n), b, c, f, g, p, r, s, machep
double precision d(n), e(n), z(n), b, c, f, g, p, r, s, machep
C      real*8 dsqrt, dabs, dsign, d1mach
double precision dsqrt, dabs, dsign
c
C      machep=d1mach(4)
machep = 2.0d0**(-52.0d0)
c
ierr = 0
if (n .eq. 1) go to 1001
c
e(n) = 0.0d0
do 240 l = 1, n
j = 0
c     :::::::::: look for small sub-diagonal element ::::::::::
105    do 110 m = l, n
if (m .eq. n) go to 120
if (dabs(e(m)) .le. machep * (dabs(d(m)) + dabs(d(m+1))))
x         go to 120
110    continue
c
120    p = d(l)
if (m .eq. l) go to 240
if (j .eq. 30) go to 1000
j = j + 1
c     :::::::::: form shift ::::::::::
g = (d(l+1) - p) / (2.0d0 * e(l))
r = dsqrt(g*g+1.0d0)
g = d(m) - p + e(l) / (g + dsign(r, g))
s = 1.0d0
c = 1.0d0
p = 0.0d0
mml = m - l
c
c     :::::::::: for i=m-1 step -1 until l do -- ::::::::::
do 200 ii = 1, mml
i = m - ii
f = s * e(i)
b = c * e(i)
if (dabs(f) .lt. dabs(g)) go to 150
c = g / f
r = dsqrt(c*c+1.0d0)
e(i+1) = f * r
s = 1.0d0 / r
c = c * s
go to 160
150       s = f / g
r = dsqrt(s*s+1.0d0)
e(i+1) = g * r
c = 1.0d0 / r
s = s * c
160       g = d(i+1) - p
r = (d(i) - g) * s + 2.0d0 * c * b
p = s * r
d(i+1) = g + p
g = c * r - b
c     :::::::::: form first component of vector ::::::::::
f = z(i+1)
z(i+1) = s * z(i) + c * f
C  200       z(i) = c * z(i) - s * f
z(i) = c * z(i) - s * f
200    continue
c
d(l) = d(l) - p
e(l) = g
e(m) = 0.0d0
go to 105
240 continue
c
c     :::::::::: order eigenvalues and eigenvectors ::::::::::
do 300 ii = 2, n
i = ii - 1
k = i
p = d(i)
c
do 260 j = ii, n
if (d(j) .ge. p) go to 260
k = j
p = d(j)
260    continue
c
if (k .eq. i) go to 300
d(k) = d(i)
d(i) = p
p = z(i)
z(i) = z(k)
z(k) = p
300 continue
c
go to 1001
c     :::::::::: set error -- no convergence to an
c                eigenvalue after 30 iterations ::::::::::
1000 ierr = l
1001 return
c     :::::::::: last card of gausq2 ::::::::::
end