From 897c8a7636c83e3ba0c31de81a5e7bce5d43c76c Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Wed, 19 Feb 2025 13:10:15 +0530 Subject: [PATCH 1/3] Fix banded-diagonal `Lmul`/`Rmul` copy --- Project.toml | 2 +- src/diagonal.jl | 26 ++++++++++++++++++++------ 2 files changed, 21 insertions(+), 7 deletions(-) diff --git a/Project.toml b/Project.toml index e5c83b2..f729a83 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ArrayLayouts" uuid = "4c555306-a7a7-4459-81d9-ec55ddd5c99a" authors = ["Sheehan Olver "] -version = "1.11.0" +version = "1.11.1" [deps] FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b" diff --git a/src/diagonal.jl b/src/diagonal.jl index 8a8b2be..fd0e6dc 100644 --- a/src/diagonal.jl +++ b/src/diagonal.jl @@ -63,28 +63,42 @@ _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) + if iszero(view(A, diagind(A, -1))) + Bidiagonal(A, :U) + elseif iszero(view(A, diagind(A, 1))) + Bidiagonal(A, :L) + else + throw(InexactError(:Bidiagonal, A)) + end +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 +_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 +_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 From 17c7dc80c84536f2b64abe1387eef379c7770435 Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Wed, 19 Feb 2025 13:38:12 +0530 Subject: [PATCH 2/3] Add tests for lower Bidiagonal --- src/diagonal.jl | 8 ++++---- test/test_layoutarray.jl | 20 +++++++++++++------- 2 files changed, 17 insertions(+), 11 deletions(-) diff --git a/src/diagonal.jl b/src/diagonal.jl index fd0e6dc..94b0830 100644 --- a/src/diagonal.jl +++ b/src/diagonal.jl @@ -65,13 +65,13 @@ _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 if iszero(view(A, diagind(A, -1))) - Bidiagonal(A, :U) - elseif iszero(view(A, diagind(A, 1))) - Bidiagonal(A, :L) + uplo = :U else - throw(InexactError(:Bidiagonal, A)) + uplo = :L end + Bidiagonal(A, uplo) end function copy(M::Rmul{<:BidiagonalLayout,<:DiagonalLayout}) A = _bidiagonal(M.A) diff --git a/test/test_layoutarray.jl b/test/test_layoutarray.jl index 3655aa0..60a664a 100644 --- a/test/test_layoutarray.jl +++ b/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 From abdb8b9ff01d0d72bbee23601c317b489ac9fa63 Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Wed, 19 Feb 2025 13:40:25 +0530 Subject: [PATCH 3/3] Add comments --- src/diagonal.jl | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/diagonal.jl b/src/diagonal.jl index 94b0830..a2c83ea 100644 --- a/src/diagonal.jl +++ b/src/diagonal.jl @@ -65,7 +65,8 @@ _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 + # we assume that the matrix is indeed bidiagonal, + # so that the conversion is lossless if iszero(view(A, diagind(A, -1))) uplo = :U else @@ -81,6 +82,8 @@ function copy(M::Lmul{<:DiagonalLayout,<:BidiagonalLayout}) 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}) @@ -91,6 +94,8 @@ function copy(M::Lmul{<:DiagonalLayout,<:TridiagonalLayout}) 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})