!*******************eivtris*****************************************
!                                                     B Bunk 11/2001
!                                                     rev     5/2014
!     eigenvalues of a real, symmetric, tridiagonal matrix
!           - derived from NR program "tqli"
!           - real(8) version
!
!           n           : matrix size
!           alpha(1:n)  : diagonal elements (input, preserved)
!           beta(1:n-1) : sub-diagonal elements (input, preserved)
!           d(1:n)      : (unsorted) eigenvalues (output)

!     for a check of the call format, include
!           use eivmod
!     in the calling program 

module eivmod
  interface
     subroutine eivtris(n, alpha, beta, d)
       integer :: n
       real(8) :: alpha(n), beta(n-1), d(n)
     end subroutine eivtris
  end interface
end module eivmod


subroutine eivtris(n, alpha, beta, d)

implicit real(8) (a-h,o-z)

dimension alpha(n), beta(n-1), d(n)
dimension e(n)

if (n < 1) stop 'error in eivtris: n < 1'
if (n == 1) then
      d(1) = alpha(1)
      return
endif

d = alpha
e(1:n-1) = beta; e(n) = 0.

!!!      subroutine tqli(d,e,n,np,z)
!!!      dimension d(np),e(np),z(np,np)
!!!      if (n.gt.1) then
!!!        do 11 i=2,n
!!!          e(i-1)=e(i)
!!!11      continue
!!!        e(n)=0.
        do 15 l=1,n
          iter=0
1         do 12 m=l,n-1
            dd=abs(d(m))+abs(d(m+1))
            if (abs(e(m))+dd.eq.dd) go to 2
12        continue
          m=n
2         if(m.ne.l)then
!!!            if(iter.eq.30)pause 'too many iterations'
            if(iter.eq.30)stop 'error: eivtris failed'
            iter=iter+1
            g=(d(l+1)-d(l))/(2.*e(l))
            r=sqrt(g**2+1.)
            g=d(m)-d(l)+e(l)/(g+sign(r,g))
            s=1.
            c=1.
            p=0.
            do 14 i=m-1,l,-1
              f=s*e(i)
              b=c*e(i)
              if(abs(f).ge.abs(g))then
                c=g/f
                r=sqrt(c**2+1.)
                e(i+1)=f*r
                s=1./r
                c=c*s
              else
                s=f/g
                r=sqrt(s**2+1.)
                e(i+1)=g*r
                c=1./r  
                s=s*c
              endif
              g=d(i+1)-p
              r=(d(i)-g)*s+2.*c*b
              p=s*r
              d(i+1)=g+p
              g=c*r-b
!!!              do 13 k=1,n
!!!                f=z(k,i+1)
!!!                z(k,i+1)=s*z(k,i)+c*f
!!!                z(k,i)=c*z(k,i)-s*f
!!!13            continue
14          continue
            d(l)=d(l)-p
            e(l)=g
            e(m)=0.
            go to 1
          endif
15      continue
!!!      endif
      return
end subroutine eivtris
