Fix banded-diagonal Lmul/Rmul copy (#255)

* Fix banded-diagonal `Lmul`/`Rmul` copy

* Add tests for lower Bidiagonal

* Add comments

Author: jishnub

Project.toml

@@ -1,7 +1,7 @@
 name = "ArrayLayouts"
 uuid = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "1.11.0"
+version = "1.11.1"
 
 [deps]
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"

src/diagonal.jl

@@ -63,28 +63,47 @@ _similar(A::SymTridiagonal) = similar(Tridiagonal(A.ev, A.dv, A.ev))
 _copy_diag(M::T, ::T) where {T<:Rmul} = copyto!(_similar(M.A), M)
 _copy_diag(M::T, ::T) where {T<:Lmul} = copyto!(_similar(M.B), M)
 _copy_diag(M, _) = copy(M)
+_bidiagonal(A::Bidiagonal) = A
+function _bidiagonal(A)
+    # we assume that the matrix is indeed bidiagonal,
+    # so that the conversion is lossless
+    if iszero(view(A, diagind(A, -1)))
+        uplo = :U
+    else
+        uplo = :L
+    end
+    Bidiagonal(A, uplo)
+end
 function copy(M::Rmul{<:BidiagonalLayout,<:DiagonalLayout})
-    A = convert(Bidiagonal, M.A)
+    A = _bidiagonal(M.A)
     _copy_diag(Rmul(A, M.B), M)
 end
 function copy(M::Lmul{<:DiagonalLayout,<:BidiagonalLayout})
-    B = convert(Bidiagonal, M.B)
+    B = _bidiagonal(M.B)
     _copy_diag(Lmul(M.A, B), M)
 end
+# we assume that the matrix is indeed tridiagonal,
+# so that the conversion is lossless
+_tridiagonal(A::Tridiagonal) = A
+_tridiagonal(A) = Tridiagonal(A)
 function copy(M::Rmul{<:TridiagonalLayout,<:DiagonalLayout})
-    A = convert(Tridiagonal, M.A)
+    A = _tridiagonal(M.A)
     _copy_diag(Rmul(A, M.B), M)
 end
 function copy(M::Lmul{<:DiagonalLayout,<:TridiagonalLayout})
-    B = convert(Tridiagonal, M.B)
+    B = _tridiagonal(M.B)
     _copy_diag(Lmul(M.A, B), M)
 end
+# we assume that the matrix is indeed symmetric tridiagonal,
+# so that the conversion is lossless
+_symtridiagonal(A::SymTridiagonal) = A
+_symtridiagonal(A) = SymTridiagonal(A)
 function copy(M::Rmul{<:SymTridiagonalLayout,<:DiagonalLayout})
-    A = convert(SymTridiagonal, M.A)
+    A = _symtridiagonal(M.A)
     _copy_diag(Rmul(A, M.B), M)
 end
 function copy(M::Lmul{<:DiagonalLayout,<:SymTridiagonalLayout})
-    B = convert(SymTridiagonal, M.B)
+    B = _symtridiagonal(M.B)
     _copy_diag(Lmul(M.A, B), M)
 end
 

test/test_layoutarray.jl

@@ -338,18 +338,24 @@ Base.copy(A::MyVector) = MyVector(copy(A.A))
 
         @testset "Diagonal * Bidiagonal/Tridiagonal with structured diags" begin
             n = size(D,1)
-            B = Bidiagonal(map(MyVector, (rand(n), rand(n-1)))..., :U)
-            MB = MyMatrix(B)
+            BU = Bidiagonal(map(MyVector, (rand(n), rand(n-1)))..., :U)
+            MBU = MyMatrix(BU)
+            BL = Bidiagonal(map(MyVector, (rand(n), rand(n-1)))..., :L)
+            MBL = MyMatrix(BL)
             S = SymTridiagonal(map(MyVector, (rand(n), rand(n-1)))...)
             MS = MyMatrix(S)
             T = Tridiagonal(map(MyVector, (rand(n-1), rand(n), rand(n-1)))...)
             MT = MyMatrix(T)
-            DA, BA, SA, TA = map(Array, (D, B, S, T))
+            DA, BUA, BLA, SA, TA = map(Array, (D, BU, BL, S, T))
             if VERSION >= v"1.11"
-                @test D * B ≈ DA * BA
-                @test B * D ≈ BA * DA
-                @test D * MB ≈ DA * BA
-                @test MB * D ≈ BA * DA
+                @test D * BU ≈ DA * BUA
+                @test BU * D ≈ BUA * DA
+                @test D * MBU ≈ DA * BUA
+                @test MBU * D ≈ BUA * DA
+                @test D * BL ≈ DA * BLA
+                @test BL * D ≈ BLA * DA
+                @test D * MBL ≈ DA * BLA
+                @test MBL * D ≈ BLA * DA
             end
             if VERSION >= v"1.12.0-DEV.824"
                 @test D * S ≈ DA * SA