Skip to content

support OffDiagonal multiplication #289

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
dlfivefifty merged 1 commit into from
Oct 13, 2026
Merged

support OffDiagonal multiplication #289

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
44 changes: 44 additions & 0 deletions src/diagonal.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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



Expand Down
2 changes: 2 additions & 0 deletions src/factorizations.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
31 changes: 31 additions & 0 deletions test/test_layoutarray.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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}()
Expand Down
Loading