I have a symmetric matrix that I am I trying to determine the eigen-values of, and the NR routines return the incorrect eigen-values.

For example, I'm trying to use the Householder reduction method on this matrix:


1 2 0
2 1 2
0 2 1


My implementation of this in code is:


MatDoub e(3,3);
e[0][0]=1;
e[0][1]=2;
e[0][2]=0;
e[1][0]=2;
e[1][1]=1;
e[1][2]=2;
e[2][0]=0;
e[2][1]=2;
e[2][2]=1;

MatrixPrint(e,1);

Symmeig evs(e);

printf("%f\n%f\n%f\n",evs.d[0],evs.d[1],evs.d[2]);


This returns the correctly formatted matrix and the eigen-values:
5.488767, -0.138240, -2.350527.

These values don't checkout as eigenvalues.

Mathcad returns the eigen-values: -1.828, 1, 3.828, which are correct.

Switching to the Unsymmeig functions, I get 3.788183, 0.000000, 1.020177 as the values.

Any advice, am I doing something wrong?

Thanks

I've found a partial cause/solution to the problem.

For reference, I'm using Visual C++ 6.0 SP5. There seems to be a problem with the C library abs function when used on a MatDoub.

For example, in this code:


sm=0.0;
for (ip=0;ip<n-1;ip++) {
for (iq=ip+1;iq<n;iq++)
{
sm += abs(a[ip][iq]);
}
}
printf("%i sm: %f\n",i,sm);
if (sm == 0.0) {
eigsrt(d,&v);
return;
}


'sm' would always be zero and the function will return without performing any rotations. If I declare a simple myabs function and replace all occurances of abs with myabs, the rotations are performed and eigenvalues are returned!


double myabs(double a)
{
if(a<0)
return -a;

return a;
}


This solution works just find for Jacobi and Symmeig, but there seems to be some issue with my fix and the Unsymmeig method. In certain cases I only get two of the eigenvalues with the other one being 0.

Hope somebody finds this useful.

Ok, another piece of follow-up. The eigenvectors are correctly calculated in Unsymmeig, but the sorting function Unsymmeig::sortvecs always seems to zero the smallest eigenvalue.

Hi Aaron,

I am not able to reproduce your problem with the g++ compiler. If you typed in the code, that is usually the cause of the problem in a case like this. Unsymmeig is a very long program, and it is difficult not to make mistakes. It is possible that this is a compiler problem otherwise. Try compiling with all optimization turned off.

Best wishes,
Saul Teukolsky