Make (Sym)Tridiag-Diagonal solves return TriDiagonal

[dkarrasch]

Closes #42684.

[dkarrasch]

This is ready for review. All tests pass, even though I'm requiring == instead of approximate equality. I guess that's (kind of?) safe due to the many structural zeros in the game.

[dkarrasch]

Thanks for your comments. I was hoping to get the matrix-eltype running as well, but that turns would require making more copies by default. Not sure if it's worth it the code bloat. Maybe if someone requests it, we can get back to it.

[N5N3]

Based on my limited testing, the current implementation is able to hande matrix-eltype:

julia> Revise.track(LinearAlgebra)
julia> D = Diagonal([I(2) for i in 1:4]);
julia> S = SymTridiagonal([randn(2,2) for i in 1:4], [randn(2,2) for i in 1:3]);
julia> D \ S;
julia> (D \ S) |> typeof
Tridiagonal{Matrix{Float64}, Vector{Matrix{Float64}}}
julia> (D \ S) == S
true
julia> (D \ S) === S
false

Edit: My fault, just notice that == between Tridiagonal and SymTridiagonal is broken for matrix eltype.
And Matrix(::SymTridiagonal) also gives incorrect result for matrix eltype.

[dkarrasch]

At the expense of "code bloat" I made the matrix-eltype-case also work, because I thought it was working before... turns out it wasn't, though. If we feel that's not quite necessary, I can revert the last commit and we go with the more concise version.

[N5N3]

The "code bloat" could be solved if we use Tr[j, j-1], Tr[j-1, j] and Tr[j, j]to get the (Sym)Tridiag's elements.
Based on my local benchmarks, llvm is clever enough to make these code branchless if we "improve" the getindex implementation.
Edit: Although LLVM is clever enough to make getindex branchless (some modification to getindex's code is needed). SymTridiag's getindex might call copy, so I think it might not be the best choise.
Instead, we can define some help function like _updiag(T, i::Int) _lowdiag(T, i::Int) and fuse the main loopbody.

BTW, we'd better add 1*1 case to the testset, as:

julia> Tridiagonal(1:0,1:1,1:0) / Diagonal(1:1)
ERROR: BoundsError: attempt to access 0-element UnitRange{Int64} at index [1]

[dkarrasch]

Finally (thanks to @N5N3's comments especially on the edge cases) this seems ready.