WIP: SymTridiagonal eigensolver

[andreasnoack]

Just to open the discussion about how we should do the SymTridiagonal. The eigen solver works well when it only uses QL steps. The first is LAPACK and the second is the new julia implementation.

julia> A=SymTridiagonal(randn(5000),randn(4999));

julia> @time eigvals(A);
elapsed time: 0.912759657 seconds (160368 bytes allocated)

julia> @time Base.LinAlg.eigvals!(copy(A));
elapsed time: 1.91420947 seconds (200408 bytes allocated)

LAPACK tries to be clever and switch between QR and QL, but for me, the QR doesn't give fast convergence.

I didn't use the Givens type because I wanted to see how fast I could make it and also because I couldn't use the A_mul_B methods anyway. There is no room for the bulge in a SymTridiagonal.

We also have to consider the vectors. There are at least two options for handling the vectors. We could do the same as LAPACK and apply the rotations in each step, or we could accumulate Givens matrices in a Rotation. We should probably also introduce a SymmetricTridiagonalFactorization immutable to keep track of all the pieces.

[andreasnoack]

Actually, the first timings were not fair to julia because eigvals(SymTridiagonal) uses a fast square root free solver in LAPACK. Comparing to the same algorithm in LAPACK, julia is actually quite a bit faster. See the new timings. Of course the LAPACK version might spend the time on some good things that gives better results, but it doesn't seem that the julia implementation gives less precise eigenvalues.

julia> A=SymTridiagonal(randn(20000),randn(19999));

julia> @time Base.LinAlg.eigvals!(copy(A)); # Julia based
elapsed time: 27.444170306 seconds (640264 bytes allocated)

julia> @time My.steqr!('N',copy(A.dv),copy(A.ev),zeros(0,0)); # Normal QR/QL in LAPACK
elapsed time: 49.393142503 seconds (640352 bytes allocated)

julia> @time eigvals(A);
elapsed time: 12.431734872 seconds (640368 bytes allocated) # Square root free QR/QL in LAPACK

[stevengj]

You might want to try some of the test cases described here and downloadable in the stetester package.

[andreasnoack]

@stevengj Thank you. I'll look into it.

[tkelman]

This is a feature, removing from 0.4.x milestone.

[andreasnoack]

This one is now being developed at https://github.com/andreasnoack/LinearAlgebra.jl