Resolve stack-overflow in Diagonal*Bidiagonal and (Sym)Tridiagonal

Author: jishnub

src/diagonal.jl

@@ -57,12 +57,12 @@ copy(M::Rdiv{<:DiagonalLayout,<:DiagonalLayout{<:AbstractFillLayout}}) = diagona
 
 
 ## bi/tridiagonal copy
-copy(M::Rmul{<:BidiagonalLayout,<:DiagonalLayout}) = convert(Bidiagonal, M.A) * M.B
-copy(M::Lmul{<:DiagonalLayout,<:BidiagonalLayout}) = M.A * convert(Bidiagonal, M.B)
-copy(M::Rmul{<:TridiagonalLayout,<:DiagonalLayout}) = convert(Tridiagonal, M.A) * M.B
-copy(M::Lmul{<:DiagonalLayout,<:TridiagonalLayout}) = M.A * convert(Tridiagonal, M.B)
-copy(M::Rmul{<:SymTridiagonalLayout,<:DiagonalLayout}) = convert(SymTridiagonal, M.A) * M.B
-copy(M::Lmul{<:DiagonalLayout,<:SymTridiagonalLayout}) = M.A * convert(SymTridiagonal, M.B)
+copy(M::Rmul{<:BidiagonalLayout,<:DiagonalLayout}) = convert(Bidiagonal, M.A) .* permutedims(parent(convert(Diagonal, M.B)))
+copy(M::Lmul{<:DiagonalLayout,<:BidiagonalLayout}) = parent(convert(Diagonal, M.A)) .* convert(Bidiagonal, M.B)
+copy(M::Rmul{<:TridiagonalLayout,<:DiagonalLayout}) = convert(Tridiagonal, M.A) .* permutedims(parent(convert(Diagonal, M.B)))
+copy(M::Lmul{<:DiagonalLayout,<:TridiagonalLayout}) = parent(convert(Diagonal, M.A)) .* convert(Tridiagonal, M.B)
+copy(M::Rmul{<:SymTridiagonalLayout,<:DiagonalLayout}) = convert(SymTridiagonal, M.A) .* permutedims(parent(convert(Diagonal, M.B)))
+copy(M::Lmul{<:DiagonalLayout,<:SymTridiagonalLayout}) = parent(convert(Diagonal, M.A)) .* convert(SymTridiagonal, M.B)
 
 copy(M::Lmul{DiagonalLayout{OnesLayout}}) = _copy_oftype(M.B, eltype(M))
 copy(M::Lmul{DiagonalLayout{OnesLayout},<:DiagonalLayout}) = Diagonal(_copy_oftype(diagonaldata(M.B), eltype(M)))

test/test_layoutarray.jl

@@ -337,6 +337,22 @@ MemoryLayout(::Type{MyVector}) = DenseColumnMajor()
         @test B̃\D ≈ B̃\Matrix(D)
         @test D\D̃ ≈ D̃\D
         @test B̃/D ≈ B̃/Matrix(D)
+
+        @testset "Diagonal * Bidiagonal/Tridiagonal with structured diags" begin
+            n = size(D,1)
+            B = Bidiagonal(map(MyVector, (rand(n), rand(n-1)))..., :U)
+            S = SymTridiagonal(map(MyVector, (rand(n), rand(n-1)))...)
+            T = Tridiagonal(map(MyVector, (rand(n-1), rand(n), rand(n-1)))...)
+            DA, BA, SA, TA = map(Array, (D, B, S, T))
+            @test D * B ≈ DA * BA
+            @test B * D ≈ BA * DA
+            if VERSION >= v"1.12.0-DEV.824"
+                @test D * S ≈ DA * SA
+                @test D * T ≈ DA * TA
+                @test S * D ≈ SA * DA
+                @test T * D ≈ TA * DA
+            end
+        end
     end
 
     @testset "Adj/Trans" begin