From 0a1f0f1af0e68ae0800f34691efcd84002842198 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Tue, 13 Oct 2026 14:09:29 +0100 Subject: [PATCH] Reapply "support OffDiagonal multiplication" This reverts commit 17f1576241fb6d3996704fec2a947b80b542405b. --- src/diagonal.jl | 44 ++++++++++++++++++++++++++++++++++++++++ src/factorizations.jl | 2 ++ test/test_layoutarray.jl | 31 ++++++++++++++++++++++++++++ 3 files changed, 77 insertions(+) diff --git a/src/diagonal.jl b/src/diagonal.jl index a2c83ea..481a7b3 100644 --- a/src/diagonal.jl +++ b/src/diagonal.jl @@ -19,6 +19,45 @@ copy(M::Lmul{<:DiagonalLayout,<:DiagonalLayout}) = diagonal(diagonaldata(M.A) .* copy(M::Lmul{<:DiagonalLayout}) = diagonaldata(M.A) .* M.B copy(M::Rmul{<:Any,<:DiagonalLayout}) = M.A .* permutedims(diagonaldata(M.B)) +### +# OffDiagonal multiplication shifts and rescales rows (Lmul) or columns (Rmul) +### + +mulreduce(M::Mul{<:OffDiagonalLayout}) = Lmul(M) +mulreduce(M::Mul{<:Any,<:OffDiagonalLayout}) = Rmul(M) +mulreduce(M::Mul{<:OffDiagonalLayout,<:OffDiagonalLayout}) = Lmul(M) +# other structured matrices already have Lmul/Rmul implementations +mulreduce(M::Mul{<:DiagonalLayout,<:OffDiagonalLayout}) = Lmul(M) +mulreduce(M::Mul{<:OffDiagonalLayout,<:DiagonalLayout}) = Rmul(M) +mulreduce(M::Mul{<:TriangularLayout,<:OffDiagonalLayout}) = Lmul(M) +mulreduce(M::Mul{<:OffDiagonalLayout,<:TriangularLayout}) = Rmul(M) + +# The entries of the band of an OffDiagonal A lie in rows _offdiagonalrows(A) and columns _offdiagonalcols(A) +_offdiagonalrows(A) = max(0, -offdiagonalindex(A)) .+ axes(offdiagonaldata(A),1) +_offdiagonalcols(A) = max(0, offdiagonalindex(A)) .+ axes(offdiagonaldata(A),1) + +_rowsview(B::AbstractVector, kr) = view(B, kr) +_rowsview(B::AbstractMatrix, kr) = view(B, kr, :) + +function copyto!(dest::AbstractArray, M::Lmul{<:OffDiagonalLayout}) + A, B = M.A, M.B + axes(dest) == axes(M) || throw(DimensionMismatch(LazyString("destination axes ", axes(dest), " do not match ", axes(M)))) + zero!(dest) + _rowsview(dest, _offdiagonalrows(A)) .= offdiagonaldata(A) .* _rowsview(B, _offdiagonalcols(A)) + dest +end + +function copyto!(dest::AbstractArray, M::Rmul{<:Any,<:OffDiagonalLayout}) + A, B = M.A, M.B + axes(dest) == axes(M) || throw(DimensionMismatch(LazyString("destination axes ", axes(dest), " do not match ", axes(M)))) + zero!(dest) + view(dest, :, _offdiagonalcols(B)) .= view(A, :, _offdiagonalrows(B)) .* permutedims(offdiagonaldata(B)) + dest +end + +copy(M::Lmul{<:OffDiagonalLayout}) = copyto!(similar(M), M) +copy(M::Rmul{<:Any,<:OffDiagonalLayout}) = copyto!(similar(M), M) + dualadjoint(_) = adjoint dualadjoint(::Transpose) = transpose @@ -27,6 +66,11 @@ function copy(M::Rmul{<:DualLayout,<:DiagonalLayout}) adj = dualadjoint(M.A) adj(adj(M.B) * adj(M.A)) end +# keep row vectors as row vectors +function copy(M::Rmul{<:DualLayout,<:OffDiagonalLayout}) + adj = dualadjoint(M.A) + adj(mul(adj(M.B), adj(M.A))) +end diff --git a/src/factorizations.jl b/src/factorizations.jl index 0673695..26ace1f 100644 --- a/src/factorizations.jl +++ b/src/factorizations.jl @@ -136,6 +136,8 @@ mulreduce(M::Mul{<:TriangularLayout,<:AbstractQLayout}) = Rmul(M) mulreduce(M::Mul{<:AbstractQLayout,<:TriangularLayout}) = Lmul(M) mulreduce(M::Mul{<:DiagonalLayout, <:AbstractQLayout}) = Lmul(M) mulreduce(M::Mul{<:AbstractQLayout, <:DiagonalLayout}) = Rmul(M) +mulreduce(M::Mul{<:AbstractQLayout, <:OffDiagonalLayout}) = Lmul(M) +mulreduce(M::Mul{<:OffDiagonalLayout, <:AbstractQLayout}) = Rmul(M) function copyto!(dest::AbstractArray{T}, M::Lmul{<:AbstractQLayout}) where T A,B = M.A,M.B diff --git a/test/test_layoutarray.jl b/test/test_layoutarray.jl index dac1dc2..7bfddd8 100644 --- a/test/test_layoutarray.jl +++ b/test/test_layoutarray.jl @@ -313,6 +313,37 @@ BroadcastStyle(::Type{<:MyVector2{T}}) where {T} = MyBroadcastStyle{1}() @test Diagonal(view(MyVector(x), 2:4))[1:2,2:3] == Diagonal(x[2:4])[1:2,2:3] end + @testset "OffDiagonal mul" begin + for O in (FillArrays.OffDiagonal(randn(4), 1), FillArrays.OffDiagonal(randn(3), -2), + FillArrays.OffDiagonal(randn(3), 2, 3, 5), FillArrays.OffDiagonal(randn(2), -1, 3, 2), + FillArrays.OffDiagonal(Fill(2.0, 3), 0, 3, 4), FillArrays.OffDiagonal(Float64[], 4, 2, 3), + view(FillArrays.OffDiagonal(randn(5), 1), 2:4, 1:6), view(Diagonal(randn(5)), 1:3, 2:5)) + m, n = size(O) + Om = Matrix(O) + b, B, C = randn(n), randn(n, 3), randn(2, m) + @test ArrayLayouts.mulreduce(Mul(O, B)) isa ArrayLayouts.Lmul + @test ArrayLayouts.mulreduce(Mul(C, O)) isa ArrayLayouts.Rmul + @test ArrayLayouts.mul(O, b) ≈ Om * b + @test ArrayLayouts.mul(O, B) ≈ Om * B + @test ArrayLayouts.mul(C, O) ≈ C * Om + @test ArrayLayouts.mul(b', O') ≈ b' * Om' + @test O * MyVector(b) ≈ Om * b + @test O * MyMatrix(B) ≈ Om * B + @test MyMatrix(C) * O ≈ C * Om + @test transpose(MyVector(randn(m))) * O isa Transpose + # mixing with other structured matrices + Dm, Dn = Diagonal(randn(m)), Diagonal(randn(n)) + Um, Un = UpperTriangular(randn(m,m)), UpperTriangular(randn(n,n)) + @test ArrayLayouts.mul(Dm, O) ≈ Dm * Om + @test ArrayLayouts.mul(O, Dn) ≈ Om * Dn + @test ArrayLayouts.mul(Um, O) ≈ Um * Om + @test ArrayLayouts.mul(O, Un) ≈ Om * Un + O2 = FillArrays.OffDiagonal(randn(max(0,n-1)), 1, n, n) + @test ArrayLayouts.mul(O, O2) ≈ Om * Matrix(O2) + end + @test_throws DimensionMismatch copyto!(zeros(2,2), ArrayLayouts.Lmul(FillArrays.OffDiagonal(randn(2), 1), randn(3,3))) + end + @testset "BroadcastStyle" begin A = MyMatrix2(randn(3,6)) @test BroadcastStyle(typeof(reshape(A, (9, 2)))) == MyBroadcastStyle{2}()