eigenmatrix of complex floating-point matrix is wrong

[dkrenn]

sage: m = Matrix(CDF, 8, [[-1, -1, -1, -1, 1, -3, -1, -1], [1, 1, 1, 1, -1, -1, 1, -3], [-1, 3, -1, -1, 1, 1, -1, -1], [-1, -1, -1, 3, 1, 1, -1, -1], [1, 1, -3, 1, -1, -1, 1, 1], [1, 1, 1, 1, -1, -1, -3, 1], [3, -1, -1, -1, 1, 1, -1, -1], [1, 1, 1, 1, 3, -1, 1, 1]])
sage: d, p = m.eigenmatrix_left()
sage: (p[1] * m)[0] / p[1][0]
1.2360679775 - 3.80422606518*I
sage: d[1][1]
1.2360679775 + 3.80422606518*I

Sage seems to return the complex conjugate of d or something of the sort. Perhaps d is simply wrongly permuted (but real eigenvalues seem to be correct).

p.inverse() * d * p should be at least approximately equal to m.

This was reported on ​the public bug reports from the notebook interface by <david+bugs@…> on 1/24/2012.

[dsm]

Somewhat smaller example, and evidence that something's wrong in RDF too:

sage: M = [[1,-1],[1, 1]]
sage: rings = SR, ZZ, QQ, RDF, CDF
sage: for ring in rings:
....:     m = Matrix(ring, M)
....:     d, p = m.left_eigenmatrix()
....:     print m- (~p)*d*p
....: 
[0 0]
[0 0]
[0 0]
[0 0]
[0 0]
[0 0]
[ 0.0 -2.0]
[ 2.0  0.0]
[2.22044604925e-16 + 1.96068977308e-16*I              -2.0 - 2.22044604925e-16*I]
[              2.0 + 8.14867788484e-17*I                      -2.22044604925e-16]

[dsm]

[Of course even if m is over RDF d and p can be over CDF.]

[rbeezer]

Appears the computed results in both of the above examples are correct for the transpose of the matrix. Digging deeper.......

[dkrenn]

Still getting

sage: norm(p.inverse() * d * p - m)
7.608452130361234

on Sage 6.6.

FYI, doing the same but over QQbar gives 1.3060006046871026e-37.