peace,

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..