symmetrize/transpose in getindex of SymTridiagonal

[dkarrasch]

Fixes JuliaLang/LinearAlgebra.jl#613.

Previously, there was no distinction between whether the second argument in the SymTridiagonal constructor was the super- or subdiagional, and whether the matrices on the diagonal are :U or :L symmetric. With this PR, this is hardcoded to :U symmetry and superdiagonal, so that the subdiagonal gets transposed. Also, the code comment in line 8 should be changed then to "superdiagonal" as a correct reminder for developers.

We could introduce an uplo::Char field, to switch between the two cases, with :U being the default? But I'm not sure if this would be breaking...

[dkarrasch]

Hm, CI is complaining that it can't find a reference for Transpose, but I thought we merged my corresponding PR...? I'm confused. CI failure seems unrelated? I'd say this is ready for review, a careful one. At first, I only adopted the getindex method, and only later realized that there are more methods that require modification, like constructors from matrices, diag. I hope I haven't missed anything.

[dkarrasch]

I think this is ready for review/merge.

[andreasnoack]

I think this one hits the same issue as #32041. I'm wondering if we should materialize the transpose to avoid the type instability.

[dkarrasch]

Do you mean we should make block SymTridiagonals be block Tridiagonals with subdiagonal the materialized transposed superdiagonal? Is there a price one would pay for "giving up" symmetric tridiagonality by type?

[dkarrasch]

I fixed the type instability, at the expense of materializing the subdiagonal elements at every getindex call. We should either add a comment such that this gets revisited once the type instability gets fully resolved, or think about constructing a TriDiagonal matrix right at the start (well, or better leave that as a design option to the user).

BTW, the getindex(::Symmetric, (2,1)) suffers from the same type instability.

[dkarrasch]

After some long detour in between, I think I figured out how to solve the problem. I see two open issues:

1. should we mention/recommend considering to generate a TriDiagonal(transpose.(ev), dv, ev) for block symtridiagonal matrices in the docstring, to avoid repeated materialization of the transposes?
2. In the constructor from matrices, we used to check whether the off-diagonals are the same. I have replaced this by checking that they are pointwise transpose of each other and whether the diagonal elements are symmetric. Should we drop any tests and simply construct SymTridiagonal(diag(A, 0), diag(A, 1)? For instance, in Symmetric, we don't check whether the input is symmetric, but make it symmetric. Similarly, we could not check anything and make SymTridiagonal(A::Matrix) be symmetric tridiagonal.

Otherwise, the tests even check for type inference... and pass.

[dkarrasch]

Is this still of interest? @dlfivefifty @andreasnoack

[andreasnoack]

Yes. Thanks for following up.