Unalias off-diagonals in Tridiagonal constructor

[jishnub]

As noted by Samuel Isaacson on slack, currently,

julia> e = ones(4); f = e[1:end-1];

julia> T = Tridiagonal(f, 2e, f)
4×4 Tridiagonal{Float64, Vector{Float64}}:
 2.0  1.0   ⋅    ⋅ 
 1.0  2.0  1.0   ⋅ 
  ⋅   1.0  2.0  1.0
  ⋅    ⋅   1.0  2.0

julia> T ./= 10
4×4 Tridiagonal{Float64, Vector{Float64}}:
 0.2   0.01   ⋅     ⋅ 
 0.01  0.2   0.01   ⋅ 
  ⋅    0.01  0.2   0.01
  ⋅     ⋅    0.01  0.2

This PR resolves such issues by ensuring that the off-diagonals are not aliased. After this,

julia> T ./= 10
4×4 Tridiagonal{Float64, Vector{Float64}}:
 0.2  0.1   ⋅    ⋅ 
 0.1  0.2  0.1   ⋅ 
  ⋅   0.1  0.2  0.1
  ⋅    ⋅   0.1  0.2

This also gets the in-place lu! to work on such an array, which otherwise wouldn't as the off-diagonals are aliased.

It was suggested in JuliaLang/LinearAlgebra.jl#706 that a warning for mutable functions may suffice (instead of an explicit check for dl !== du), but preventing the aliasing seems to be a better solution that would avoid various bugs.

[mbauman]

I like this much better, but it looks like you'll also need to revert #35197 (as those tests added there should fail).

[dkarrasch]

I believe we have tests in the lu! code (and perhaps elsewhere) that check for potential aliasing. Those could be removed. Moreover, this might help to simplify recent BunchKaufman -eigen-related efforts. @aravindh-krishnamoorthy

EDIT: 🤦 That's exactly what @mbauman was talking about.

[dkarrasch]

Should we document which off-diagonal is aliased to the argument of the constructor and which one may be just a copy?

[jishnub]

I feel we should document that the off diagonals shouldn't be aliased, and leave the de-aliasing as an implementation detail.

[aravindh-krishnamoorthy]

Hello, as a Julia beginner and new user, I feel I might prefer to have a warning here along with the unaliasing if I alias dl and du (and not an error as earlier). The reason is that, in Julia, I always expect matrix assignments and calls to reuse the same memory. Silent unaliasing of dl and du would be (from my point of view) unexpected behaviour.

[jishnub]

I've added a comment to the docstring explaining that du is copied if aliasing is detected.

[aravindh-krishnamoorthy]

Thank you very much, @jishnub. I'm not sure if it makes sense, in my comment above, I was actually referring to a runtime warning that would seem most correct to me as a beginner.

[dkarrasch]

I think runtime warnings may ruin performance. Moreover, we have the hint in the documentation, so it's not exactly "silent".

[dkarrasch]

Is this ready to go? I'd say so.

[jishnub]

Yeah, seems ok to me

[PallHaraldsson]

This seems better since you prevent a bug, wrong results.

I'm just curious this makes a copy sometimes, e.g. in that case I presume, and why didn't it throw in the old code?

julia/base/abstractarray.jl

Line 1521 in ea261ce

_unaliascopy(A, C) = throw(ArgumentError("""

Sometimes you do not want allocations, and that is now a possibility, but also still the possibility of throwing (for unusual types).

[dkarrasch]

It didn't throw in the old code because there was no aliasing check included. It's true, it may now imply an allocation that may not be intended by the user, but it appears only once at construction, and afterwards makes no difference except that it avoids errors.