tred2 has been giving me no joy for the past day. Since I think it's safe to assume that the code itself is not at error, i'm probably doing something wrong with the inputs i'm giving tred2.

In the comments for the C version of tred2, the function takes in a matrix a[1..n][1..n] and two arrays d[1..n] and e[1..n]. Am I correct in interpreting that to mean that the 0th row and column of the matrix a are ignored? For example, if I wanted to tridiagonalise the matrix

1 2 3
1 2 3
1 2 3

a[][] would be

0 0 0 0
0 1 2 3
0 1 2 3
0 1 2 3

Likewise, i'm assuming that d[] and e[] also ignore the 0th element and so are really d[0,1..n] and e[0,1..n]. For this example, i'm taking n to be 3.

PS: On second thought, that's a horrible example since the matrix is not real symmetric. Let's just ignore that for the moment. ;)

Ok try in that way,

void foo(void)
{
float **a,*d,*e;

a = matrix(1,3,1,3);
d = vector(1,3);
e = vector(1,3);

/*filling a, the values in this example are random.
Use your data */

a[1][1]=1;
a[2][1]=2;
a[3][1]=3;
a[1][2]=2;
a[2][2]=2;
a[3][2]=2;
a[1][3]=3;
a[2][3]=3;
a[3][3]=3;

/* d */

d[1] =1;
d[2] =2;
d[3] =3;

/* e */

e[1] =1;
e[2]=2;
e[3]=3;

/* now do what you want */

.
.
.
.
.

/* before exiting you have to free the memory */

free_vector(e,1,3);
free_vector(d,1,3);
free_matrix(a,1,3,1,3);

return;

}



There is a dedicate section in NR to describe the allocation and
deallocation steps. Read it!

Basically any element (i,j) of the matrix a allocated with

a = matrix(n,m,l,o) can be used as

a[i][j]

with i from n to m and j from l to o.


and any element i of the vector v allocated with

v = vector(l,m) can be used as

v[i]

with i from l to m.

It's trivial to convert the float routines to double.


Hope this helps

Regards

Max

Thanks for your reply gaspari. I read through the section in NR on the allocation of vectors and matrices. I seem to be doing exactly the same thing, though with my own code rather than 'matrix' and 'vector'. When I substitute 'matrix' and 'vector' in for me own code I get exactly the same (incorrect?) result.

Can I get someone to try and run tred2 two on the following matrix:

4.0 2.0 2.0 1.0
2.0 -3.0 1.0 1.0
2.0 1.0 3.0 1.0
1.0 1.0 1.0 2.0

If my reference and matlab are not wrong, the result should be

4 -3 0 0
-3 2 3.1623 0
0 3.1623 -1.4 -0.2
0 0 -0.2 1.4

However tred2 gives me the following for the diagonal and the subdiagonal:

diagonal
1.423077 -2.089745 4.666668 2.000000

subdiagonal
0.000000 0.599556 -3.399347 -1.732051

Either tred2 is doing something other than what I think it does or i'm doing something wrong.