I have been translating some Numerical recipies programs to matlab
from C and Fortran.  Fortran is obviously closer to matlab than C.

I have successfully implemented tred2 and have been trying to
implement tqli...

I have been getting almost correct results from tqli (with treds
outputs as its input), but there is just one small error either in
syntax translation by me or I don't know what, and I cannot seem to
get it to function 100% correctly...  I have spent days at this and it
is driving me crazy now...

The fortran (77) which actually translates very nicely, i thought is

  SUBROUTINE tqli(d,e,n,np,z)
      INTEGER n,np
      REAL d(np),e(np),z(np,np)
CU    USES pythag
      INTEGER i,iter,k,l,m
      REAL b,c,dd,f,g,p,r,s,pythag
      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 in tqli'
          iter=iter+1
          g=(d(l+1)-d(l))/(2.*e(l))
          r=pythag(g,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)
            r=pythag(f,g)
            e(i+1)=r
            if(r.eq.0.)then
              d(i+1)=d(i+1)-p
              e(m)=0.
              goto 1
            endif
            s=f/r
            c=g/r
            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     Omit lines from here ...
            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     ... to here when finding only eigenvalues.
14        continue
          d(l)=d(l)-p
          e(l)=g
          e(m)=0.
          goto 1
        endif
15    continue
      return
      END

my matlab is

function [d,z] = sceigen_test (d,e,z)

%input
%d is the diagonal vector of the tridiag matrix
%e is the sub and super diagonal of the tridiag matrix
%z is the transformation matrix from the tridiag program

% output,
%d is the eigenvalue vector
% z is the eigenvector matrix
% d is the eigenvale vector

% defining loop size based upon ip vector/matrix size
n = size(d,2);

%shifting e to the left, so that e(1) is not the zero anymore
for i = 2:n
    e(i-1) = e(i);
end
e(n)=0;
for l = 1:n
    iter=0;
    for m = l:n % not to n-1 for reason expained below
        if m == n
            break
        end
        dd = abs(d(m)) + abs(d(m+1));
        if ((abs(e(m))+ dd) == dd)
            break
        end
    end

   % m = n% this is implemented in the above loop, so the loop goes to
n instead of n-1

    if m ~= l
        if iter == 30
            error('bloody typical..too many iterations');
        end
        iter = iter+1;
        g = (d(l+1)-d(l))/(2.0*e(l));
        r = pythag(g, 1.0);%it works correctly
        g = d(m) - d(l) + (e(l) / (g + SIGNO(r, g)));%SIGNO = SIGN...
definitely
        s = 1.0;
        c = 1.0;
        p = 0.0;
        for i = m-1:-1:l
            f = s*e(i);
            b = c*e(i);
            r = pythag(f, g);
            e(i+1) = r;
            if r == 0.0
                d(i+1)= d(i+1)-p;
                e(m) = 0.0;
                break
            end
            s = f/r;
            c = g/r;
            g = d(i + 1) - p;
            r = (d(i) - g)*s + 2.0*c*b;
            p = s*r;
            d(i + 1) = g + p;
            g = c * r - b;
            %eigenvector part now...
            for 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;
            end
        end
        if r ~= 0.0
            d(l) = d(l) - p;
            e(l) = g;
            e(m) = 0.0;
        end
    end
end

The input (example)

which is (if anyone will test for themselves)

[0.0861,0.0172,0.0466,-0.0186,-0.0384;0.0172,0.0635,-0.0146,-0.0154,0.0399; 0.0466,-0.0146,0.0555,-0.0166,-0.0563;-0.0186,-0.0154,-0.0166,0.0459,-0.000 7;-0.0384,0.0399,-0.0563,-0.0007,0.0704;]

which is a covariance matrix, so it is symmetric etc...

eigenvalues (ordered hi to low)

eig_sort =

  Columns 1 through 3

   0.171639537105982                   0
0
                   0                              0.100853395915414
0
                   0
0                              0.037197771768781
                   0
0                              0
                   0
0                              0

  Columns 4 through 5

                   0                               0
                   0                               0
                   0                               0
   0.011737362737036                    0
                   0                               0.000002142655106

I have eigenvectors also, but i want to focus on theeigenvalues for
now

Tha actual matlab generates eigenvalues (reverse ordered) are:

de =

  Columns 1 through 3

   0.000000000000000                   0
0
                   0                              0.011739410520354
0
                   0
0                              0.037197816074304
                   0
0                              0
                   0
0                              0

  Columns 4 through 5

                   0                   0
                   0                   0
                   0                   0
   0.100853446274617        0
                   0                   0.171639537313043

They are similar, but there always a constant error...  I noticed
also, that if I subtract eig_sort(5,5) from itself and add it to
eig_sort(4,4), the results become more accurate...

But this is probably just indicative of some other problem....

It would be great if someone could see what I was doing wrong and
point it out..

Here they are again...

mine are

0.171639537105982      0.100853395915414      0.037197771768781
0.011737362737036      0.000002142655106

matlab's are (in reverse order to my own)