Error in hqr2?


glk63
12-06-2011, 04:43 PM
Apparently there is an error in code hqr2
Look at the code fragment:
...
if (anorm! = 0.0) {
for (nn = n-1; nn> = 0; nn --) {
p = real (wri [nn]);
q = imag (wri [nn]);
na = nn-1;
...
If nn = 0, the value na =- 1, but then:
...
} Else if (q <0.0) {
m = na;
if (abs (a [nn] [na])> abs (a [na] [nn])) {
a [na] [na] = q / a [nn] [na];
...
Thus, when nn = 0 there is a reference to a [0] [-1]. For example, when calculating the eigenvalues and eigenvectors of the matrix:
A = {{7,1,2}, {2,3,1}, {-8, -2, -1}} an incorrect addressing error occurs in hqr2 code.

Saul Teukolsky
12-07-2011, 06:39 PM
I'm unable to reproduce an addressing error in hqr2 with the matrix A. The potential reference to a[0][-1] when nn=0 should never be triggered. The reason is that complex conjugate eigenvalues are stored so that the one with negative imaginary part has an index 1 higher than the one with positive imaginary part. This is implemented by the code
z=sqrt(abs(q));
...
wri[nn]=Complex(x+p,-z);
wri[nn-1]=conj(wri[nn]);
This means if there is a complex conjugate pair of eigenvalues with indexes 0 and 1, the code you were worried about
q=imag(wri[nn]);
na=nn-1;
is always executed with q<0 for nn=1 and q>0 for nn=0. But only the q<0 case goes into the block where a[nn][na] is referenced, so there is never a problem.

I must admit this is pretty obscure: it goes back to the original Eispack code.

The matrix A is difficult for hqr2 because it is defective - there is only one independent eigenvector. You can see this, as suggested in the text, by noting that two of the purported eigenvectors that are output are parallel.

You can check for out-of-bounds references by compiling with the corresponding NR flag set, e.g.
g++ -D_CHECKBOUNDS_

Hope this helps.
Saul Teukolsky