`ldiv!` assumes `StridedArray` has unit stride?

[mcabbott]

The following bug turned up in #40504:

julia> using LinearAlgebra
julia> U = UpperTriangular(reshape(collect(3:27.0),5,5));
julia> B = reshape(collect(1:14.0),2,7); B2 = copy(B); B3 = copy(B)
2×7 Matrix{Float64}:
 1.0  3.0  5.0  7.0   9.0  11.0  13.0
 2.0  4.0  6.0  8.0  10.0  12.0  14.0

julia> R1 = ldiv!(U, B[1, 1:5])
5-element Vector{Float64}:
 -0.689790319419949
 -0.20693709582598488
 -0.11640211640211637
 -0.07936507936507933
  0.3333333333333333

julia> R2 = ldiv!(U, view(B2, 1, 1:5))  # not right, StridedVecOrMat -> LAPACK.trtrs!
5-element view(::Matrix{Float64}, 1, 1:5) with eltype Float64:
 -0.3372308228275303
 -0.056907701352145804
  0.18518518518518517
  7.0
  9.0

julia> R3 = ldiv!(U, view(B3', 1:5, 1))  # fallback naivesub!(U, v), correct
5-element view(adjoint(::Matrix{Float64}), 1:5, 1) with eltype Float64:
 -0.689790319419949
 -0.20693709582598488
 -0.11640211640211637
 -0.07936507936507933
  0.3333333333333333

julia> B2  # written into outside the view
2×7 Matrix{Float64}:
 -0.337231  -0.0569077  0.185185  7.0   9.0  11.0  13.0
 -0.101169  -0.0388007  6.0       8.0  10.0  12.0  14.0

julia> B3
2×7 Matrix{Float64}:
 -0.68979  -0.206937  -0.116402  -0.0793651   0.333333  11.0  13.0
  2.0       4.0        6.0        8.0        10.0       12.0  14.0

julia> view(B2, 1, 1:5) isa StridedArray
true

julia> strides(view(B2, 1, 1:5)) == (2,)  # thus not really acceptable to LAPACK function?
true

It appears that #39301 worked around this, by introducing enough adjoints to hit the fallback method. Which #39467 then broke by simplifying some views of adjoint matrices. The minimal work-around would be for #39301 to directly call naivesub!, pending a solution to StridedArray dispatch?

cc. @schneiderfelipe

[schneiderfelipe]

Hi, @mcabbott thanks for looking into that. I think using naivesub! is reasonable.

@andreasnoack, what do you think?

[andreasnoack]

It's unfortunate that we don't have proper dispatch for this. We probably need a fallback call to naivesub! to handle this generally but the vector case can actually be handled by BLAS so if split https://github.com/JuliaLang/julia/blob/728aa90906a8712668120bb6f0912f963678ec63/stdlib/LinearAlgebra/src/triangular.jl#L739-L740 into a method for StridedVector and one for StridedMatrix then the vector version can call BLAS.trsv! since it allows for non-unit stride. That wouldn't solve the problem for something like view(randn(4,4), 1:2:3, 1:2:3) so we probably need to introduce some ugly branching here.

This also reveals that a stride check is missing from https://github.com/JuliaLang/julia/blob/728aa90906a8712668120bb6f0912f963678ec63/stdlib/LinearAlgebra/src/lapack.jl#L3412-L3433

[mcabbott]

Has anyone looked into correcting #39301, so that its test do not block otherwise correct behaviour?

[mcabbott]

Still broken on 1.9, thus still blocks https://github.com/JuliaLang/julia/pull/43719

Maybe https://github.com/JuliaLang/julia/pull/39301 should work around this in a less fragile way?