Inaccurate eigenvectors with dstedc

[caliarim]

Dear all,

I think dstedc produces inaccurate eigenvectors for a specific tridiagonal symmetric matrix. In order to show it, I assembled the attached test. It computes eigenpairs of the Jacobi matrix related to Gauss-Hermite quadrature formulas (Golub-Welsch algorithm). Then, it applies the quadrature formula to int_{-infty}^{infty} x^34 * exp(-x^2) dx. The quadrature formula should be exact for 18 or more points. It is numerically exact if the eigenpairs are computed with dsteqr. But it diverges (as the number of points increases) with dstedc. I tested lapack 3.7.1.
tsteig.f.zip

[langou]

Hi @caliarim

Jim Demmel followed up on this. Here is what Jim says:

It seems that Matlab (using dstevd) gets O(macheps) values for
norm(VV'-eye) and norm(JV-V*D) where [V,D]=eig(J), which is
all we promise. The user prefers the result of V from dstev, in
particular the entries of V much smaller than macheps, which
is something that we can't make any general guarantees about
of course, it depends on the special structure of your matrix.
Allowing an option to use dstev instead of dstevd sounds like
the best option for those expert users who know they have such
special structure. We can also refer the user to some literature
about what special structures we have studied, but that certainly
won't cover all cases.

One more comment on the special structure of this problem.
Since the matrix is tridiagonal with 0 diagonal, it is really
equivalent to the bidiagonal SVD, which we have extensively
studied. This is an "easy" problem in the sense that the
condition number (for m=100) is 121, and the minimum
relative gap between eigenvalues/singular values is .031,
so all the eigen/singular values are determined to high
relative accuracy, and so are the eigen/singular vectors,
in a norm sense. But the the user wants componentwise
high relative accuracy, even for the tiniest components,
which we at least did not analyze/promise.

Jim

Cheers,
Julien.