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..ab8e369 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 @@ -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)) @@ -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))...) @@ -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}, @@ -958,9 +959,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)) diff --git a/src/blockkron.jl b/src/blockkron.jl index 8630334..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) @@ -201,6 +204,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 = @@ -223,6 +227,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 +248,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(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 + + krontravargs(K::KronTrav) = K.args function krontravargs(V::SubArray) KR,JR = parentindices(V) 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