'****************************************************
'*  Calculate the eigenvalues and the eigenvectors  *
'*      of a real symmetric tridiagonal matrix      *
'* ------------------------------------------------ *
'* Ref.: "NUMERICAL RECIPES, Cambridge University   *
'*        Press, 1986" [BIBLI 08].                  *
'*                                                  *
'* SAMPLE RUN:                                      *
'* (find eigenvalues and eigenvectors of matrix:    *
'*          1  2  0  0     0                        *
'*          2  4  7  0     0                        *
'*     A =  0  7 10  8     0                        *
'*          0  0  8 -0.75 -9                        *
'*          0  0  0 -9    10  )                     *
'*                                                  *
'* Eigenvalues 1:   -0.78071442                     *
'* Eigenvector:                                     *
'* 1.000000 -0.890357  0.322363  0.344649  0.287721 *
'*                                                  *
'* Eigenvalues 2:    2.53046878                     *
'* Eigenvector:                                     *
'* 1.000000  0.765234 -0.446362 -0.252815 -0.304616 *
'*                                                  *
'* Eigenvalues 3:   -9.15659229                     *
'* Eigenvector:                                     *
'*-0.056408  0.286458 -0.522285  1.000000  0.469812 *
'*                                                  *
'* Eigenvalues 4:   12.53989066                     *
'* Eigenvector:                                     *
'* 0.097161  0.560616  0.656182 -0.282210  1.000000 *
'*                                                  *
'* Eigenvalues 5:   19.11694726                     *
'* Eigenvector:                                     *
'* 0.051876  0.469920  1.000000  0.728439 -0.719095 *
'*                                                  *
'*              Basic version by J-P Moreau, Paris  *
'*                       (www.jpmoreau.fr)          *
'****************************************************
'PROGRAM TEST_TQLI
defdbl a-h,o-z
defint i-n

NMAX=20
dim D(NMAX), E(NMAX), Z(NMAX,NMAX)

cls
print
input " Input size of matrix: ", n

gosub 1000  'read data

'initialize Z to identity matrix
for i=1 to n
  for j=1 to n
    if i<>j then
      Z(i,j) = 0#
    else
      Z(i,j) = 1#
    end if
  next j
next i

gosub 2000 'call TQLI(D,E,n,Z)

'print results
for j=1 to n
  print
  print " Eigenvalue ";j;": ";
  print using "####.######"; D(j)
  print " Eigenvector:"
'normalize vector before printing
  xmax=Z(1,j)
  for i=2 to n
    if ABS(Z(i,j))>ABS(xmax) then xmax=Z(i,j)
  next i
  for i=1 to n
    print using "####.######"; Z(i,j)/xmax;
  next i
  print
next j

print
print

END 'of main program

'read vectors D and E from screen
'Subroutine Read_data(n,D,E)
1000  print
  print " Input "; n-1; " elements of subdiagonal:"
  E(1)=1#     'arbitrary value
  for i=1 to n-1
    print "  Element ";i; ": "; : input "", E(i+1)
  next i
  print " Input "; n; " elements of main diagonal:"
  for i=1 to n
    print "  Element ";i; ": "; : input "", D(i)
  next i
  print
return

1500 'y=Sign(r,g)
  if g>=0 then
    y=ABS(r)
  else
    y=-ABS(r)
  end if
return

'****************************************************
'* This subroutine implements the QL algorithm with *
'* implicit shifts to determine the eigenvalues and *
'* eigenvectors of a real symmetric tridiagonal     *
'* matrix.                                          *
'* ------------------------------------------------ *
'* INPUTS:                                          *
'*         D    Elements of main diagonal           *
'*         E    Elements of subdiagonal (from E(2)  *
'*              to E(n))                            *
'*         n    size of matrix                      *
'*         Z    identity matrix                     *
'* OUTPUTS                                          *
'*         D    vector storing the n eigenvalues    *
'*         Z    matrix storing the n eigenvectors   *
'*              in columns.                         *
'****************************************************
'Subroutine TQLI(D,E,n,Z)
2000 if n>1 then
  for i=2 to n
    E(i-1)=E(i)
  next i
  E(n)=0#
  for l=1 to n
    iter=0
1   for m=l to n-1
      dd=ABS(D(m))+ABS(D(m+1))
      if ABS(E(m))+dd=dd then goto 2
    next m
    m=n
2   if m<>l then
      if iter=30 then
        print " 30 iterations."
        return
      end if
      iter=iter+1
      g=(D(l+1)-D(l))/(2#*E(l))
      r=SQR(g*g+1#)
      gosub 1500 'y=Sign(r,g)
      g=D(m)-D(l)+E(l)/(g+y)
      s=1# : c=1# : p=0#
      for i=m-1 to l step -1
        f=s*E(i) : b=c*E(i)
        if ABS(f)>=ABS(g) then
          c=g/f
	  r=SQR(c*c+1#)
	  E(i+1)=f*r
	  s=1#/r
	  c=c*s
        else
	  s=f/g
	  r=SQR(s*s+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
        for k=1 to 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
        next k
      next i
      D(l)=D(l)-p
      E(l)=g : E(m)=0#
      goto 1
    end if
  next l
end if
return

'end of file elprotd.bas