From b09ca02dd4ce2ddad5fba5ad2ef779ffd4b5f761 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Fri, 28 Apr 2023 16:47:49 +0100 Subject: [PATCH 1/4] Make copying a view of a KronTrav faster --- Project.toml | 2 +- src/LazyBandedMatrices.jl | 10 +++++----- src/blockkron.jl | 11 +++++++++++ 3 files changed, 17 insertions(+), 6 deletions(-) diff --git a/Project.toml b/Project.toml index a7af6a6..e1f9523 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "LazyBandedMatrices" uuid = "d7e5e226-e90b-4449-9968-0f923699bf6f" authors = ["Sheehan Olver "] -version = "0.8.9" +version = "0.8.10" [deps] ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a" diff --git a/src/LazyBandedMatrices.jl b/src/LazyBandedMatrices.jl index c282cdb..caf5df7 100644 --- a/src/LazyBandedMatrices.jl +++ b/src/LazyBandedMatrices.jl @@ -41,7 +41,7 @@ import BlockBandedMatrices: BlockSlice, Block1, AbstractBlockBandedLayout, subblockbandwidths, BandedBlockBandedMatrix, BlockBandedMatrix, BlockBandedLayout, AbstractBandedBlockBandedLayout, BandedBlockBandedLayout, BandedBlockBandedStyle, BlockBandedStyle, blockcolsupport, BlockRange1, blockrowsupport, BlockIndexRange1, - BlockBandedColumnMajor + BlockBandedColumnMajor, blockbanded_copyto! import BlockArrays: BlockSlice1, BlockLayout, AbstractBlockStyle, block, blockindex, BlockKron, viewblock, blocks, BlockSlices, AbstractBlockLayout, blockvec # for bidiag/tridiag @@ -582,15 +582,15 @@ _mulbanded_copyto!(dest::AbstractArray{T}, a, b) where T = muladd!(one(T), a, b, _mulbanded_copyto!(dest::AbstractArray{T}, a, b, c, d...) where T = _mulbanded_copyto!(dest, mul(a,b), c, d...) _mulbanded_BandedMatrix(A, _) = A -_mulbanded_BandedMatrix(A, ::NTuple{2,OneTo{Int}}) = BandedMatrix(A) -_mulbanded_BandedMatrix(A) = _mulbanded_BandedMatrix(A, axes(A)) +_mulbanded_BandedMatrix(A, ::NTuple{2,Int}) = BandedMatrix(A) +_mulbanded_BandedMatrix(A) = _mulbanded_BandedMatrix(A, size(A)) _copyto!(::AbstractBandedLayout, ::ApplyBandedLayout{typeof(*)}, dest::AbstractMatrix, src::AbstractMatrix) = _mulbanded_copyto!(dest, map(_mulbanded_BandedMatrix,arguments(src))...) _mulbanded_BandedBlockBandedMatrix(A, _) = A -_mulbanded_BandedBlockBandedMatrix(A, ::NTuple{2,OneTo{Int}}) = BandedBlockBandedMatrix(A) -_mulbanded_BandedBlockBandedMatrix(A) = _mulbanded_BandedBlockBandedMatrix(A, axes(A)) +_mulbanded_BandedBlockBandedMatrix(A, ::NTuple{2,Int}) = BandedBlockBandedMatrix(A) +_mulbanded_BandedBlockBandedMatrix(A) = _mulbanded_BandedBlockBandedMatrix(A, size(A)) _copyto!(::AbstractBandedBlockBandedLayout, ::ApplyBandedBlockBandedLayout{typeof(*)}, dest::AbstractMatrix, src::AbstractMatrix) = _mulbanded_copyto!(dest, map(_mulbanded_BandedBlockBandedMatrix,arguments(src))...) diff --git a/src/blockkron.jl b/src/blockkron.jl index 8630334..d46e47d 100644 --- a/src/blockkron.jl +++ b/src/blockkron.jl @@ -223,6 +223,7 @@ isbandedblockbanded(A::KronTrav) = isblockbanded(A) && length(A.args) == 2 convert(::Type{B}, A::KronTrav{<:Any,2}) where B<:BandedBlockBandedMatrix = convert(B, BandedBlockBandedMatrix(A)) struct KronTravBandedBlockBandedLayout <: AbstractBandedBlockBandedLayout end +struct SubKronTravBandedBlockBandedLayout <: AbstractBandedBlockBandedLayout end struct KronTravLayout{M} <: AbstractBlockLayout end @@ -243,6 +244,16 @@ sublayout(::KronTravBandedBlockBandedLayout, ::Type{<:NTuple{2,BlockSlice{BlockR sub_materialize(::Union{KronTravLayout,KronTravBandedBlockBandedLayout}, V) = KronTrav(map(sub_materialize, krontravargs(V))...) +_copyto!(_, ::KronTravBandedBlockBandedLayout, dest::AbstractMatrix, src::AbstractMatrix) = blockbanded_copyto!(dest, sub_materialize(src)) +function _copyto!(_, ::SubKronTravBandedBlockBandedLayout, dest::AbstractMatrix, src::AbstractMatrix) + KR, JR = parentindices(src) + # use Base.OneTo range which is a Kron Trav + KR2,JR2 = Block.(Base.OneTo(Int(KR.block[end]))),Block.(Base.OneTo(Int(JR.block[end]))) + # materialize parent with KR2,JR2 ranges then copyto! + copyto!(dest, view(BandedBlockBandedMatrix(view(parent(src),KR2, JR2)), KR, JR)) +end + + krontravargs(K::KronTrav) = K.args function krontravargs(V::SubArray) KR,JR = parentindices(V) From 2665aa54f0e4ee5e692ccc89f09eddf44204da55 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Mon, 1 May 2023 08:44:39 +0100 Subject: [PATCH 2/4] Add Padded * Structured --- src/LazyBandedMatrices.jl | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/LazyBandedMatrices.jl b/src/LazyBandedMatrices.jl index caf5df7..bea9526 100644 --- a/src/LazyBandedMatrices.jl +++ b/src/LazyBandedMatrices.jl @@ -958,9 +958,13 @@ copy(M::Mul{BroadcastLayout{typeof(*)}, <:StructuredLazyLayouts}) = lazymaterial ## padded copy mulreduce(M::Mul{<:StructuredLazyLayouts, <:PaddedLayout}) = MulAdd(M) mulreduce(M::Mul{<:StructuredApplyLayouts{F}, D}) where {F,D<:PaddedLayout} = Mul{ApplyLayout{F},D}(M.A, M.B) +mulreduce(M::Mul{<:PaddedLayout, <:StructuredLazyLayouts}) = MulAdd(M) +mulreduce(M::Mul{D, <:StructuredApplyLayouts{F}}) where {F,D<:PaddedLayout} = Mul{D,ApplyLayout{F}}(M.A, M.B) # need to overload copy due to above copy(M::Mul{<:StructuredLazyLayouts, <:PaddedLayout}) = copy(mulreduce(M)) +copy(M::Mul{<:PaddedLayout, <:StructuredLazyLayouts}) = copy(mulreduce(M)) simplifiable(::Mul{<:StructuredLazyLayouts, <:PaddedLayout}) = Val(true) +simplifiable(::Mul{<:PaddedLayout, <:StructuredLazyLayouts}) = Val(true) copy(L::Ldiv{ApplyBandedLayout{typeof(*)}, Lay}) where Lay = copy(Ldiv{ApplyLayout{typeof(*)},Lay}(L.A, L.B)) From c1047ff613cb168d1ad8e09b77539f883ac5ad83 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Mon, 1 May 2023 21:45:51 +0100 Subject: [PATCH 3/4] Bidiagonal layout --- src/LazyBandedMatrices.jl | 3 ++- src/blockkron.jl | 1 + 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/src/LazyBandedMatrices.jl b/src/LazyBandedMatrices.jl index bea9526..f4312c4 100644 --- a/src/LazyBandedMatrices.jl +++ b/src/LazyBandedMatrices.jl @@ -881,7 +881,8 @@ BandedLazyLayouts = Union{AbstractLazyBandedLayout, BandedColumns{LazyLayout}, B TriangularLayout{UPLO,UNIT,BandedRows{LazyLayout}} where {UPLO,UNIT}, TriangularLayout{UPLO,UNIT,BandedColumns{LazyLayout}} where {UPLO,UNIT}, SymTridiagonalLayout{LazyLayout}, BidiagonalLayout{LazyLayout}, TridiagonalLayout{LazyLayout}, - SymmetricLayout{BandedColumns{LazyLayout}}, HermitianLayout{BandedColumns{LazyLayout}}} + SymmetricLayout{BandedColumns{LazyLayout}}, HermitianLayout{BandedColumns{LazyLayout}}, + BidiagonalLayout{LazyLayout,LazyLayout}} StructuredLazyLayouts = Union{BandedLazyLayouts, BlockBandedColumns{LazyLayout}, BandedBlockBandedColumns{LazyLayout}, diff --git a/src/blockkron.jl b/src/blockkron.jl index d46e47d..15e7459 100644 --- a/src/blockkron.jl +++ b/src/blockkron.jl @@ -201,6 +201,7 @@ function _krontrav_getindex(Kin::Block{2}, A, B, C) end getindex(M::KronTrav{<:Any,2}, K::Block{2}) = _krontrav_getindex(K, M.args...) +# getindex(M::KronTrav{<:Any,2}, K::BlockIndex{2}) = _krontrav_getindex(K, M.args...) getindex(A::KronTrav{<:Any,N}, kj::Vararg{Int,N}) where N = From 60725f361646ec144e0f8ee2cf1c413650764356 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Tue, 2 May 2023 21:42:37 +0100 Subject: [PATCH 4/4] Increase coverage --- src/LazyBandedMatrices.jl | 4 ++-- src/blockkron.jl | 5 ++++- test/test_banded.jl | 9 ++++++++- 3 files changed, 14 insertions(+), 4 deletions(-) diff --git a/src/LazyBandedMatrices.jl b/src/LazyBandedMatrices.jl index f4312c4..ab8e369 100644 --- a/src/LazyBandedMatrices.jl +++ b/src/LazyBandedMatrices.jl @@ -533,10 +533,10 @@ _broadcastarray2broadcasted(::StructuredBroadcastLayouts{F}, A) where F = _broad _broadcastarray2broadcasted(::StructuredBroadcastLayouts{F}, A::BroadcastArray) where F = _broadcastarray2broadcasted(BroadcastLayout{F}(), A) _copyto!(::AbstractBandedLayout, ::BroadcastBandedLayout, dest::AbstractMatrix, bc::AbstractMatrix) = - copyto!(dest, _broadcastarray2broadcasted(bc)) + copyto!(dest, _broadcastarray2broadcasted(MemoryLayout(bc), bc)) _copyto!(_, ::BroadcastBandedLayout, dest::AbstractMatrix, bc::AbstractMatrix) = - copyto!(dest, _broadcastarray2broadcasted(bc)) + copyto!(dest, _broadcastarray2broadcasted(MemoryLayout(bc), bc)) _banded_broadcast!(dest::AbstractMatrix, f, (A,B)::Tuple{AbstractMatrix{T},AbstractMatrix{V}}, _, ::Tuple{<:Any,ApplyBandedLayout{typeof(*)}}) where {T,V} = broadcast!(f, dest, BandedMatrix(A), BandedMatrix(B)) diff --git a/src/blockkron.jl b/src/blockkron.jl index 15e7459..2c41b5e 100644 --- a/src/blockkron.jl +++ b/src/blockkron.jl @@ -64,6 +64,9 @@ function colsupport(A::DiagTrav{<:Any,2}, _) axes(A,1)[bs] end +copy(A::DiagTrav) = DiagTrav(copy(A.array)) +==(A::DiagTrav, B::PseudoBlockVector) = PseudoBlockVector(A) == B +==(A::PseudoBlockVector, B::DiagTrav) = A == PseudoBlockVector(B) function getindex(A::DiagTrav, K::Block{1}) @boundscheck checkbounds(A, K) @@ -249,7 +252,7 @@ _copyto!(_, ::KronTravBandedBlockBandedLayout, dest::AbstractMatrix, src::Abstra function _copyto!(_, ::SubKronTravBandedBlockBandedLayout, dest::AbstractMatrix, src::AbstractMatrix) KR, JR = parentindices(src) # use Base.OneTo range which is a Kron Trav - KR2,JR2 = Block.(Base.OneTo(Int(KR.block[end]))),Block.(Base.OneTo(Int(JR.block[end]))) + KR2,JR2 = Block.(Base.OneTo(Int(last(KR.block)))),Block.(Base.OneTo(Int(last(JR.block)))) # materialize parent with KR2,JR2 ranges then copyto! copyto!(dest, view(BandedBlockBandedMatrix(view(parent(src),KR2, JR2)), KR, JR)) end diff --git a/test/test_banded.jl b/test/test_banded.jl index 3d2f043..0c11623 100644 --- a/test/test_banded.jl +++ b/test/test_banded.jl @@ -1,4 +1,4 @@ - +using LazyBandedMatrices, LazyArrays, BandedMatrices, Test struct PseudoBandedMatrix{T} <: AbstractMatrix{T} data::Array{T} @@ -67,12 +67,19 @@ LinearAlgebra.lmul!(β::Number, A::PseudoBandedMatrix) = (lmul!(β, A.data); A) b = Vcat(randn(2), Zeros(3)) @test A*b ≈ Matrix(A)b @test B*b ≈ Matrix(B)b + + P = Vcat(randn(2,5), Zeros(3,5)) + @test A*P ≈ Matrix(A)P + @test P*A ≈ P*Matrix(A) end @testset "Apply * Banded" begin B = brand(5,5,2,1) A = ApplyArray(*, B, B) @test A * Vcat([1,2], Zeros(3)) ≈ B*B*[1,2,0,0,0] + P = Vcat(randn(2,5), Zeros(3,5)) + @test A*P ≈ Matrix(A)P + @test P*A ≈ P*Matrix(A) end @testset "Banded Perturbed" begin