cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c from Numerical Recipes of Press et al. (1996)
c
c a  (input) is the symmetric matrix with n columns and n rows
c n  (input) is the number true (as opposed to declared) dimension of a
c np (input) is the dimension of the matrix a defined in the calling program np >= n
c d  (output) to be used in TQLI
c e  (output) to be used in TQLI
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	subroutine tred2(a,n,np,d,e)
	dimension a(np,np),d(np),e(np)
	if (n .gt. 1) then
		do 18 i=n,2,-1
			l = i-1
			h = 0.
			scale = 0.
			if (l .gt. 1) then
				do 11 k=1,l
					scale = scale+abs(a(i,k))
11				continue
				if (scale .eq. 0.) then
					e(i) = a(i,l)
				else
					do 12 k=1,l
						a(i,k) = a(i,k)/scale
						h = h+a(i,k)**2
12					continue
					f = a(i,l)
					g = -sign(sqrt(h),f)
					e(i) = scale*g
					h = h-f*g
					a(i,l) = f-g
					f = 0.
					do 15 j=1,l 
						a(j,i) = a(i,j)/h
						g = 0.
						do 13 k=1,j
							g = g+a(j,k)*a(i,k)
13						continue
						if (l.gt.j) then
							do 14 k=j+1,l
								g = g+a(k,j)*a(i,k)
14							continue
						end if
						e(j) = g/h
						f =f+e(j)*a(i,j)
15					continue
					hh = f/(h+h)
					do 17 j=1,l
						f = a(i,j)
						g = e(j)-hh*f
						e(j) = g
						do 16 k=1,j
							a(j,k) = a(j,k)-f*e(k)-g*a(i,k)
16						continue
17					continue
				end if
			else
				e(i) = a(i,l)
			end if
			d(i) = h
18		continue
	end if
	d(1) = 0.
	e(1) = 0.
	do 23 i=1,n
		l = i-1
		if (d(i) .ne. 0.) then
			do 21 j=1,l
				g = 0.
				do 19 k=1,l
					g = g+a(i,k)*a(k,j)
19				continue
				do 20 k=1,l
					a(k,j) = a(k,j)-g*a(k,i)
20				continue
21			continue
		end if
		d(i) = a(i,i)
		a(i,i) = 1.
		if (l.ge.1) then
			do 22 j=1,l
				a(i,j) = 0.
				a(j,i) = 0.
22			continue
		end if
23	continue
	return
	end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c from Numerical Recipes of Press et al. (1996)
c
c d  (input) generated by TRED2; contains eigenvalues on output
c e  (input) generated by TRED2
c n  (input) is the number true (as opposed to declared) dimension of d and e
c np (input) is the dimension of d and e as defined in the calling program np >= n
c z  (output) is the matrix containing n eigenvectors each of length n 
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	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) goto 2
12			continue
			m = n
2			if (m .ne. l) then
				if (iter.eq.30) pause 'too many iterations'
				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) .gt. 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
					end if
					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
c form eigenvectors
					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
c end eigenvectors
14				continue
				d(l) = d(l)-p
				e(l) = g
				e(m) = 0.
				goto 1
			end if
15		continue
	end if
	return
	end