diff --git a/Project.toml b/Project.toml index 725c2b7..e2050fb 100644 --- a/Project.toml +++ b/Project.toml @@ -1,38 +1,41 @@ name = "InfiniteArrays" uuid = "4858937d-0d70-526a-a4dd-2d5cb5dd786c" -version = "0.14.3" +version = "0.15.0" [deps] ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a" +BandedMatrices = "aae01518-5342-5314-be14-df237901396f" FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b" Infinities = "e1ba4f0e-776d-440f-acd9-e1d2e9742647" LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" [weakdeps] +BandedMatrices = "aae01518-5342-5314-be14-df237901396f" DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" +[extensions] +InfiniteArraysBandedMatricesExt = "BandedMatrices" +InfiniteArraysDSPExt = "DSP" +InfiniteArraysStatisticsExt = "Statistics" + [compat] Aqua = "0.8" ArrayLayouts = "1.8" -BandedMatrices = "1.0" -BlockArrays = "1.0" +BandedMatrices = "1.7.6" Base64 = "1" +BlockArrays = "1.0" DSP = "0.7" FillArrays = "1.0" Infinities = "0.1.1" -LazyArrays = "2.0.2" +LazyArrays = "2.2.3" LinearAlgebra = "1.6" SparseArrays = "1" Statistics = "1" Test = "1" julia = "1.10" -[extensions] -InfiniteArraysDSPExt = "DSP" -InfiniteArraysStatisticsExt = "Statistics" - [extras] Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" BandedMatrices = "aae01518-5342-5314-be14-df237901396f" diff --git a/ext/InfiniteArraysBandedMatricesExt.jl b/ext/InfiniteArraysBandedMatricesExt.jl new file mode 100644 index 0000000..7be24b8 --- /dev/null +++ b/ext/InfiniteArraysBandedMatricesExt.jl @@ -0,0 +1,526 @@ +module InfiniteArraysBandedMatricesExt +using InfiniteArrays, BandedMatrices, LinearAlgebra +using InfiniteArrays.LazyArrays, InfiniteArrays.ArrayLayouts, InfiniteArrays.FillArrays + +import Base: BroadcastStyle, size, getindex, similar, copy, *, +, -, /, \, materialize!, copyto!, OneTo +import Base.Broadcast: Broadcasted +import InfiniteArrays: InfIndexRanges, Infinity, PosInfinity, OneToInf, InfAxes, AbstractInfUnitRange, InfRanges +import ArrayLayouts: sub_materialize, MemoryLayout, sublayout, mulreduce, triangularlayout, MatLdivVec, subdiagonaldata, diagonaldata, supdiagonaldata +import LazyArrays: applybroadcaststyle, applylayout, islazy, islazy_layout, simplifiable, AbstractLazyLayout, PaddedColumns, LazyArrayStyle, ApplyLayout, AbstractLazyBandedLayout, ApplyBandedLayout, BroadcastBandedLayout +import BandedMatrices: _BandedMatrix, AbstractBandedMatrix, banded_similar, BandedMatrix, bandedcolumns, BandedColumns, bandeddata +import FillArrays: AbstractFillMatrix, AbstractFill, getindex_value + +BroadcastStyle(::Type{<:SubArray{<:Any,2,<:AbstractBandedMatrix,<:Tuple{<:InfIndexRanges,<:InfIndexRanges}}})= LazyArrayStyle{2}() +_BandedMatrix(data::AbstractMatrix{T}, ::Infinity, l, u) where T = _BandedMatrix(data, ℵ₀, l, u) + +### +# BandIndexing +### + +struct InfBandCartesianIndices <: AbstractVector{CartesianIndex{2}} + b::Int +end + +InfBandCartesianIndices(b::Band) = InfBandCartesianIndices(b.i) + +size(::InfBandCartesianIndices) = (∞,) +getindex(B::InfBandCartesianIndices, k::Int) = B.b ≥ 0 ? CartesianIndex(k, k+B.b) : CartesianIndex(k-B.b, k) + +Base.checkindex(::Type{Bool}, ::NTuple{2,OneToInf{Int}}, ::InfBandCartesianIndices) = true +BandedMatrices.band_to_indices(_, ::NTuple{2,OneToInf{Int}}, b) = (InfBandCartesianIndices(b),) +BroadcastStyle(::Type{<:SubArray{<:Any,1,<:Any,Tuple{InfBandCartesianIndices}}}) = LazyArrayStyle{1}() + +_inf_banded_sub_materialize(_, V) = V +function _inf_banded_sub_materialize(::BandedColumns, V) + A = parent(V) + b = parentindices(V)[1].b + data = bandeddata(A) + l,u = bandwidths(A) + if -l ≤ b ≤ u + data[u+1-b, max(1,b+1):end] + else + Zeros{eltype(V)}(∞) # Not type stable + end +end + +sub_materialize(_, V::SubArray{<:Any,1,<:AbstractMatrix,Tuple{InfBandCartesianIndices}}, ::Tuple{InfAxes}) = + _inf_banded_sub_materialize(MemoryLayout(parent(V)), V) + + + + + + + +const TriToeplitz{T} = Tridiagonal{T,Fill{T,1,Tuple{OneToInf{Int}}}} +const ConstRowMatrix{T} = ApplyMatrix{T,typeof(*),<:Tuple{<:AbstractVector,<:AbstractFillMatrix{<:Any,Tuple{OneTo{Int},OneToInf{Int}}}}} +const PertConstRowMatrix{T} = Hcat{T,<:Tuple{Array{T},<:ConstRowMatrix{T}}} +const InfToeplitz{T} = BandedMatrix{T,<:ConstRowMatrix{T},OneToInf{Int}} +const PertToeplitz{T} = BandedMatrix{T,<:PertConstRowMatrix{T},OneToInf{Int}} + +const SymTriPertToeplitz{T} = SymTridiagonal{T,Vcat{T,1,Tuple{Vector{T},Fill{T,1,Tuple{OneToInf{Int}}}}}} +const TriPertToeplitz{T} = Tridiagonal{T,Vcat{T,1,Tuple{Vector{T},Fill{T,1,Tuple{OneToInf{Int}}}}}} +const AdjTriPertToeplitz{T} = Adjoint{T,Tridiagonal{T,Vcat{T,1,Tuple{Vector{T},Fill{T,1,Tuple{OneToInf{Int}}}}}}} +const InfBandedMatrix{T,V<:AbstractMatrix{T}} = BandedMatrix{T,V,OneToInf{Int}} + +_prepad(p, a) = Vcat(Zeros{eltype(a)}(max(p,0)), a) +_prepad(p, a::Zeros{T,1}) where T = Zeros{T}(length(a)+p) +_prepad(p, a::Ones{T,1}) where T = Ones{T}(length(a)+p) +_prepad(p, a::AbstractFill{T,1}) where T = Fill{T}(getindex_value(a), length(a)+p) + +banded_similar(T, (m,n)::Tuple{Int,PosInfinity}, (l,u)::Tuple{Int,Int}) = BandedMatrix{T}(undef, (n,m), (u,l))' + +function BandedMatrix{T}(kv::Tuple{Vararg{Pair{<:Integer,<:AbstractVector}}}, + ::NTuple{2,PosInfinity}, + (l,u)::NTuple{2,Integer}) where T + ks = getproperty.(kv, :first) + rws = Vcat(permutedims.(_prepad.(ks,getproperty.(kv, :second)))...) + c = zeros(T, l+u+1, length(ks)) + for (k,j) in zip(u .- ks .+ 1,1:length(ks)) + c[k,j] = one(T) + end + _BandedMatrix(ApplyArray(*,c,rws), ℵ₀, l, u) +end + +# Construct InfToeplitz +function BandedMatrix{T}(kv::Tuple{Vararg{Pair{<:Integer,<:Fill{<:Any,1,Tuple{OneToInf{Int}}}}}}, + mn::NTuple{2,PosInfinity}, + lu::NTuple{2,Integer}) where T + m,n = mn + @assert isinf(n) + l,u = lu + t = zeros(T, u+l+1) + for (k,v) in kv + p = length(v) + t[u-k+1] = v.value + end + + return _BandedMatrix(t * Ones{T}(1,∞), Integer(m), l, u) +end + +function BandedMatrix{T}(kv::Tuple{Vararg{Pair{<:Integer,<:Vcat{<:Any,1,<:Tuple{<:AbstractVector,Fill{<:Any,1,Tuple{OneToInf{Int}}}}}}}}, + mn::NTuple{2,PosInfinity}, + lu::NTuple{2,Integer}) where T + m,n = mn + @assert isinf(n) + l,u = lu + M = mapreduce(x -> length(x.second.args[1]) + max(0,x.first), max, kv) # number of data rows + data = zeros(T, u+l+1, M) + t = zeros(T, u+l+1) + for (k,v) in kv + a,b = v.args + p = length(a) + t[u-k+1] = b.value + if k ≤ 0 + data[u-k+1,1:p] = a + data[u-k+1,p+1:end] .= b.value + else + data[u-k+1,k+1:k+p] = a + data[u-k+1,k+p+1:end] .= b.value + end + end + + return _BandedMatrix(Hcat(data, t * Ones{T}(1,∞)), Integer(m), l, u) +end + + +function BandedMatrix(Ac::Adjoint{T,<:InfToeplitz}) where T + A = parent(Ac) + l,u = bandwidths(A) + a = A.data.args[1] + _BandedMatrix(reverse(conj(a)) * Ones{T}(1,∞), ℵ₀, u, l) +end + +function BandedMatrix(Ac::Transpose{T,<:InfToeplitz}) where T + A = parent(Ac) + l,u = bandwidths(A) + a = A.data.args[1] + _BandedMatrix(reverse(a) * Ones{T}(1,∞), ℵ₀, u, l) +end + +function BandedMatrix(Ac::Adjoint{T,<:PertToeplitz}) where T + A = parent(Ac) + l,u = bandwidths(A) + a,b = A.data.args + Ac_fd = BandedMatrix(_BandedMatrix(Hcat(a, b[:,1:l+1]), size(a,2)+l, l, u)') + _BandedMatrix(Hcat(Ac_fd.data, reverse(conj(b.args[1])) * Ones{T}(1,∞)), ℵ₀, u, l) +end + +function BandedMatrix(Ac::Transpose{T,<:PertToeplitz}) where T + A = parent(Ac) + l,u = bandwidths(A) + a,b = A.data.args + Ac_fd = BandedMatrix(transpose(_BandedMatrix(Hcat(a, b[:,1:l+1]), size(a,2)+l, l, u))) + _BandedMatrix(Hcat(Ac_fd.data, reverse(b.args[1]) * Ones{T}(1,∞)), ℵ₀, u, l) +end + + +for op in (:-, :+) + @eval begin + function $op(A::SymTriPertToeplitz{T}, λ::UniformScaling) where T + TV = promote_type(T,eltype(λ)) + dv = Vcat(convert.(AbstractVector{TV}, A.dv.args)...) + ev = Vcat(convert.(AbstractVector{TV}, A.ev.args)...) + SymTridiagonal(broadcast($op, dv, Ref(λ.λ)), ev) + end + function $op(λ::UniformScaling, A::SymTriPertToeplitz{V}) where V + TV = promote_type(eltype(λ),V) + SymTridiagonal(convert(AbstractVector{TV}, broadcast($op, Ref(λ.λ), A.dv)), + convert(AbstractVector{TV}, broadcast($op, A.ev))) + end + function $op(A::SymTridiagonal{T,<:AbstractFill}, λ::UniformScaling) where T + TV = promote_type(T,eltype(λ)) + SymTridiagonal(convert(AbstractVector{TV}, broadcast($op, A.dv, Ref(λ.λ))), + convert(AbstractVector{TV}, A.ev)) + end + + function $op(A::TriPertToeplitz{T}, λ::UniformScaling) where T + TV = promote_type(T,eltype(λ)) + Tridiagonal(Vcat(convert.(AbstractVector{TV}, A.dl.args)...), + Vcat(convert.(AbstractVector{TV}, broadcast($op, A.d, λ.λ).args)...), + Vcat(convert.(AbstractVector{TV}, A.du.args)...)) + end + function $op(λ::UniformScaling, A::TriPertToeplitz{V}) where V + TV = promote_type(eltype(λ),V) + Tridiagonal(Vcat(convert.(AbstractVector{TV}, broadcast($op, A.dl.args))...), + Vcat(convert.(AbstractVector{TV}, broadcast($op, λ.λ, A.d).args)...), + Vcat(convert.(AbstractVector{TV}, broadcast($op, A.du.args))...)) + end + function $op(adjA::AdjTriPertToeplitz{T}, λ::UniformScaling) where T + A = parent(adjA) + TV = promote_type(T,eltype(λ)) + Tridiagonal(Vcat(convert.(AbstractVector{TV}, A.du.args)...), + Vcat(convert.(AbstractVector{TV}, broadcast($op, A.d, λ.λ).args)...), + Vcat(convert.(AbstractVector{TV}, A.dl.args)...)) + end + function $op(λ::UniformScaling, adjA::AdjTriPertToeplitz{V}) where V + A = parent(adjA) + TV = promote_type(eltype(λ),V) + Tridiagonal(Vcat(convert.(AbstractVector{TV}, broadcast($op, A.du.args))...), + Vcat(convert.(AbstractVector{TV}, broadcast($op, λ.λ, A.d).args)...), + Vcat(convert.(AbstractVector{TV}, broadcast($op, A.dl.args))...)) + end + + function $op(λ::UniformScaling, A::InfToeplitz{V}) where V + l,u = bandwidths(A) + TV = promote_type(eltype(λ),V) + a = convert(AbstractVector{TV}, $op.(A.data.args[1])) + a[u+1] += λ.λ + _BandedMatrix(a*Ones{TV}(1,∞), ℵ₀, l, u) + end + + function $op(A::InfToeplitz{T}, λ::UniformScaling) where T + l,u = bandwidths(A) + TV = promote_type(T,eltype(λ)) + a = TV[Zeros{TV}(max(-u,0)); A.data.args[1]; Zeros{TV}(max(-l,0))] + a[max(0,u)+1] = $op(a[max(u,0)+1], λ.λ) + _BandedMatrix(a*Ones{TV}(1,∞), ℵ₀, max(l,0), max(u,0)) + end + + function $op(λ::UniformScaling, A::PertToeplitz{V}) where V + l,u = bandwidths(A) + TV = promote_type(eltype(λ),V) + a, t = map(AbstractArray{TV}, A.data.args) + b = $op.(t.args[1]) + a[u+1,:] .+= λ.λ + b[u+1] += λ.λ + _BandedMatrix(Hcat(a, b*Ones{TV}(1,∞)), ℵ₀, l, u) + end + + function $op(A::PertToeplitz{T}, λ::UniformScaling) where T + l,u = bandwidths(A) + TV = promote_type(T,eltype(λ)) + ã, t = A.data.args + a = AbstractArray{TV}(ã) + b = AbstractVector{TV}(t.args[1]) + a[u+1,:] .= $op.(a[u+1,:],λ.λ) + b[u+1] = $op(b[u+1], λ.λ) + _BandedMatrix(Hcat(a, b*Ones{TV}(1,∞)), ℵ₀, l, u) + end + end +end + + + +#### +# Conversions to BandedMatrix +#### + +function BandedMatrix(A::PertToeplitz{T}, (l,u)::Tuple{Int,Int}) where T + @assert A.u == u # Not implemented + a, b = A.data.args + t = b.args[1] # topelitz part + t_pad = vcat(t,Zeros(l-A.l)) + data = Hcat([vcat(a,Zeros{T}(l-A.l,size(a,2))) repeat(t_pad,1,l)], t_pad * Ones{T}(1,∞)) + _BandedMatrix(data, ℵ₀, l, u) +end + +function BandedMatrix(A::SymTriPertToeplitz{T}, (l,u)::Tuple{Int,Int}) where T + a,a∞ = A.dv.args + b,b∞ = A.ev.args + n = max(length(a), length(b)+1) + 1 + data = zeros(T, l+u+1, n) + data[u,2:length(b)+1] .= b + data[u,length(b)+2:end] .= b∞.value + data[u+1,1:length(a)] .= a + data[u+1,length(a)+1:end] .= a∞.value + data[u+2,1:length(b)] .= b + data[u+2,length(b)+1:end] .= b∞.value + _BandedMatrix(Hcat(data, [Zeros{T}(u-1); b∞.value; a∞.value; b∞.value; Zeros{T}(l-1)] * Ones{T}(1,∞)), ℵ₀, l, u) +end + +function BandedMatrix(A::SymTridiagonal{T,Fill{T,1,Tuple{OneToInf{Int}}}}, (l,u)::Tuple{Int,Int}) where T + a∞ = A.dv + b∞ = A.ev + n = 2 + data = zeros(T, l+u+1, n) + data[u,2:end] .= b∞.value + data[u+1,1:end] .= a∞.value + data[u+2,1:end] .= b∞.value + _BandedMatrix(Hcat(data, [Zeros{T}(u-1); b∞.value; a∞.value; b∞.value; Zeros{T}(l-1)] * Ones{T}(1,∞)), ℵ₀, l, u) +end + +function BandedMatrix(A::TriPertToeplitz{T}, (l,u)::Tuple{Int,Int}) where T + a,a∞ = A.d.args + b,b∞ = A.du.args + c,c∞ = A.dl.args + n = max(length(a), length(b)+1, length(c)-1) + 1 + data = zeros(T, l+u+1, n) + data[u,2:length(b)+1] .= b + data[u,length(b)+2:end] .= b∞.value + data[u+1,1:length(a)] .= a + data[u+1,length(a)+1:end] .= a∞.value + data[u+2,1:length(c)] .= c + data[u+2,length(c)+1:end] .= c∞.value + _BandedMatrix(Hcat(data, [Zeros{T}(u-1); b∞.value; a∞.value; c∞.value; Zeros{T}(l-1)] * Ones{T}(1,∞)), ℵ₀, l, u) +end + +function BandedMatrix(A::Tridiagonal{T,Fill{T,1,Tuple{OneToInf{Int}}}}, (l,u)::Tuple{Int,Int}) where T + a∞ = A.d + b∞ = A.du + c∞ = A.dl + n = 2 + data = zeros(T, l+u+1, n) + data[u,2:end] .= b∞.value + data[u+1,1:end] .= a∞.value + data[u+2,1:end] .= c∞.value + _BandedMatrix(Hcat(data, [Zeros{T}(u-1); b∞.value; a∞.value; c∞.value; Zeros{T}(l-1)] * Ones{T}(1,∞)), ℵ₀, l, u) +end + +function InfToeplitz(A::Tridiagonal{T,Fill{T,1,Tuple{OneToInf{Int}}}}, (l,u)::Tuple{Int,Int}) where T + a∞ = A.d + b∞ = A.du + c∞ = A.dl + _BandedMatrix([Zeros{T}(u-1); b∞.value; a∞.value; c∞.value; Zeros{T}(l-1)] * Ones{T}(1,∞), ℵ₀, l, u) +end + +InfToeplitz(A::Tridiagonal{T,Fill{T,1,Tuple{OneToInf{Int}}}}) where T = InfToeplitz(A, bandwidths(A)) + + +#### +# Toeplitz layout +#### + +_pertdata(A::ConstRowMatrix{T}) where T = Array{T}(undef,size(A,1),0) +_pertdata(A::Hcat{T,<:Tuple{Vector{T},<:ConstRowMatrix{T}}}) where T = 1 +_pertdata(A::PertConstRowMatrix) = A.args[1] +function _pertdata(A::SubArray) + P = parent(A) + kr,jr = parentindices(A) + dat = _pertdata(P) + dat[kr,jr ∩ axes(dat,2)] +end + +_constrows(A::ConstRowMatrix) = A.args[1]*getindex_value(A.args[2]) +_constrows(A::PertConstRowMatrix) = _constrows(A.args[2]) +_constrows(A::SubArray) = _constrows(parent(A))[parentindices(A)[1]] + +ConstRowMatrix(A::AbstractMatrix{T}) where T = ApplyMatrix(*, A[:,1], Ones{T}(1,size(A,2))) +PertConstRowMatrix(A::AbstractMatrix{T}) where T = + Hcat(_pertdata(A), ApplyMatrix(*, _constrows(A), Ones{T}(1,size(A,2)))) + +struct ConstRows <: AbstractLazyLayout end +struct PertConstRows <: AbstractLazyLayout end +MemoryLayout(::Type{<:ConstRowMatrix}) = ConstRows() +MemoryLayout(::Type{<:PertConstRowMatrix}) = PertConstRows() +bandedcolumns(::ConstRows) = BandedToeplitzLayout() +bandedcolumns(::PertConstRows) = PertToeplitzLayout() +sublayout(::ConstRows, inds...) = sublayout(ApplyLayout{typeof(*)}(), inds...) +sublayout(::PertConstRows, inds...) = sublayout(ApplyLayout{typeof(hcat)}(), inds...) +for Typ in (:ConstRows, :PertConstRows) + @eval begin + sublayout(::$Typ, ::Type{<:Tuple{Any,AbstractInfUnitRange{Int}}}) = $Typ() # no way to lose const rows + applybroadcaststyle(::Type{<:AbstractMatrix}, ::$Typ) = LazyArrayStyle{2}() + applylayout(::Type, ::$Typ, _...) = LazyLayout() + end +end + +""" + TridiagonalToeplitzLayout + +represents a matrix which is tridiagonal and toeplitz. Must support +`subdiagonalconstant`, `diagonalconstant`, `supdiagonalconstant`. +""" +struct TridiagonalToeplitzLayout <: AbstractLazyBandedLayout end +const BandedToeplitzLayout = BandedColumns{ConstRows} +const PertToeplitzLayout = BandedColumns{PertConstRows} +const PertTriangularToeplitzLayout{UPLO,UNIT} = TriangularLayout{UPLO,UNIT,BandedColumns{PertConstRows}} +struct BidiagonalToeplitzLayout <: AbstractLazyBandedLayout end +struct PertBidiagonalToeplitzLayout <: AbstractLazyBandedLayout end +struct PertTridiagonalToeplitzLayout <: AbstractLazyBandedLayout end + +const InfToeplitzLayouts = Union{TridiagonalToeplitzLayout, BandedToeplitzLayout, BidiagonalToeplitzLayout, + PertToeplitzLayout, PertTriangularToeplitzLayout, PertBidiagonalToeplitzLayout, PertTridiagonalToeplitzLayout} + +subdiagonalconstant(A) = getindex_value(subdiagonaldata(A)) +diagonalconstant(A) = getindex_value(diagonaldata(A)) +supdiagonalconstant(A) = getindex_value(supdiagonaldata(A)) + + +islazy_layout(::InfToeplitzLayouts) = Val(true) +islazy(::BandedMatrix{<:Any,<:Any,OneToInf{Int}}) = Val(true) + + + +_BandedMatrix(::BandedToeplitzLayout, A::AbstractMatrix) = + _BandedMatrix(ConstRowMatrix(bandeddata(A)), size(A,1), bandwidths(A)...) +_BandedMatrix(::PertToeplitzLayout, A::AbstractMatrix) = + _BandedMatrix(PertConstRowMatrix(bandeddata(A)), size(A,1), bandwidths(A)...) + +# for Lay in (:BandedToeplitzLayout, :PertToeplitzLayout) +# @eval begin +# sublayout(::$Lay, ::Type{<:Tuple{AbstractInfUnitRange{Int},AbstractInfUnitRange{Int}}}) = $Lay() +# sublayout(::$Lay, ::Type{<:Tuple{Slice,AbstractInfUnitRange{Int}}}) = $Lay() +# sublayout(::$Lay, ::Type{<:Tuple{AbstractInfUnitRange{Int},Slice}}) = $Lay() +# sublayout(::$Lay, ::Type{<:Tuple{Slice,Slice}}) = $Lay() + +# sub_materialize(::$Lay, V) = BandedMatrix(V) +# end +# end + + +@inline sub_materialize(::ApplyBandedLayout{typeof(*)}, V, ::Tuple{InfAxes,InfAxes}) = V +@inline sub_materialize(::BroadcastBandedLayout, V, ::Tuple{InfAxes,InfAxes}) = V +@inline sub_materialize(::BandedColumns, V, ::Tuple{InfAxes,InfAxes}) = BandedMatrix(V) +@inline sub_materialize(::BandedColumns, V, ::Tuple{InfAxes,OneTo{Int}}) = BandedMatrix(V) + +## +# UniformScaling +## + +# for op in (:+, :-), Typ in (:(BandedMatrix{<:Any,<:Any,OneToInf{Int}}), +# :(Adjoint{<:Any,<:BandedMatrix{<:Any,<:Any,OneToInf{Int}}}), +# :(Transpose{<:Any,<:BandedMatrix{<:Any,<:Any,OneToInf{Int}}})) +# @eval begin +# $op(A::$Typ, λ::UniformScaling) = $op(A, Diagonal(Fill(λ.λ,∞))) +# $op(λ::UniformScaling, A::$Typ) = $op(Diagonal(Fill(λ.λ,∞)), A) +# end +# end + + +_default_banded_broadcast(bc::Broadcasted, ::Tuple{<:OneToInf,<:Any}) = copy(Broadcasted{LazyArrayStyle{2}}(bc.f, bc.args)) + +### +# Banded * Banded +### + +BandedMatrix{T}(::UndefInitializer, axes::Tuple{OneToInf{Int},OneTo{Int}}, lu::NTuple{2,Integer}) where T = BandedMatrix{T}(undef, map(length,axes), lu) +BandedMatrix{T}(::UndefInitializer, (m,n)::Tuple{Infinity,Int}, lu::NTuple{2,Integer}) where T = BandedMatrix{T}(undef, (ℵ₀,n), lu) + +similar(M::MulAdd{<:AbstractBandedLayout,<:AbstractBandedLayout}, ::Type{T}, axes::Tuple{OneTo{Int},OneToInf{Int}}) where T = + transpose(BandedMatrix{T}(undef, reverse(axes), reverse(bandwidths(M)))) +similar(M::MulAdd{<:AbstractBandedLayout,<:AbstractBandedLayout}, ::Type{T}, axes::Tuple{OneToInf{Int},OneTo{Int}}) where T = + BandedMatrix{T}(undef, axes, bandwidths(M)) + + +### +# BandedFill * BandedFill +### + +simplifiable(::Mul{<:BandedColumns{<:AbstractFillLayout},<:BandedColumns{<:AbstractFillLayout}}) = Val(true) +copy(M::MulAdd{<:BandedColumns{<:AbstractFillLayout},<:BandedColumns{<:AbstractFillLayout},ZerosLayout}) = + _bandedfill_mul(M, axes(M.A), axes(M.B)) + +_bandedfill_mul(M, _, _) = copyto!(similar(M), M) +function _bandedfill_mul(M::MulAdd, ::Tuple{InfAxes,InfAxes}, ::Tuple{InfAxes,InfAxes}) + A, B = M.A, M.B + Al,Au = bandwidths(A) + Bl,Bu = bandwidths(B) + l,u = Al+Bl,Au+Bu + m = min(Au+Al,Bl+Bu)+1 + λ = getindex_value(bandeddata(A))*getindex_value(bandeddata(B)) + ret = _BandedMatrix(Hcat(Array{typeof(λ)}(undef, l+u+1,max(0,u)), [1:m-1; Fill(m,l+u-2m+3); m-1:-1:1]*Fill(λ,1,∞)), ℵ₀, l, u) + mul!(view(ret, 1:l+u,1:u), view(A,1:l+u,1:u+Bl), view(B,1:u+Bl,1:u)) + ret +end + +_bandedfill_mul(M::MulAdd, ::Tuple{InfAxes,Any}, ::Tuple{Any,Any}) = ApplyArray(*, M.A, M.B) +_bandedfill_mul(M::MulAdd, ::Tuple{InfAxes,Any}, ::Tuple{Any,InfAxes}) = ApplyArray(*, M.A, M.B) +_bandedfill_mul(M::MulAdd, ::Tuple{Any,Any}, ::Tuple{Any,InfAxes}) = ApplyArray(*, M.A, M.B) +_bandedfill_mul(M::MulAdd, ::Tuple{Any,InfAxes}, ::Tuple{InfAxes,InfAxes}) = ApplyArray(*, M.A, M.B) +_bandedfill_mul(M::MulAdd, ::Tuple{InfAxes,InfAxes}, ::Tuple{InfAxes,Any}) = ApplyArray(*, M.A, M.B) +_bandedfill_mul(M::MulAdd, ::Tuple{Any,InfAxes}, ::Tuple{InfAxes,Any}) = ApplyArray(*, M.A, M.B) + +mulreduce(M::Mul{<:InfToeplitzLayouts, <:InfToeplitzLayouts}) = ApplyArray(M) +mulreduce(M::Mul{<:InfToeplitzLayouts}) = ApplyArray(M) +mulreduce(M::Mul{<:InfToeplitzLayouts,<:PaddedColumns}) = MulAdd(M) +mulreduce(M::Mul{<:Any, <:InfToeplitzLayouts}) = ApplyArray(M) +mulreduce(M::Mul{<:AbstractQLayout, <:InfToeplitzLayouts}) = ApplyArray(M) +simplifiable(::Mul{<:DiagonalLayout, <:InfToeplitzLayouts}) = Val(true) +simplifiable(::Mul{<:InfToeplitzLayouts, <:DiagonalLayout}) = Val(true) +mulreduce(M::Mul{<:DiagonalLayout, <:InfToeplitzLayouts}) = Lmul(M) +mulreduce(M::Mul{<:InfToeplitzLayouts, <:DiagonalLayout}) = Rmul(M) + + + +### +# Inf-Toeplitz layout +# this could possibly be avoided via an InfFillLayout +### + +const InfFill = AbstractFill{<:Any,1,<:Tuple{OneToInf}} + +for Typ in (:(Tridiagonal{<:Any,<:InfFill}), + :(SymTridiagonal{<:Any,<:InfFill})) + @eval begin + MemoryLayout(::Type{<:$Typ}) = TridiagonalToeplitzLayout() + BroadcastStyle(::Type{<:$Typ}) = LazyArrayStyle{2}() + end +end + +MemoryLayout(::Type{<:Bidiagonal{<:Any,<:InfFill}}) = BidiagonalToeplitzLayout() +BroadcastStyle(::Type{<:Bidiagonal{<:Any,<:InfFill}}) = LazyArrayStyle{2}() + +*(A::Bidiagonal{<:Any,<:InfFill}, B::Bidiagonal{<:Any,<:InfFill}) = + mul(A, B) + +# fall back for Ldiv +triangularlayout(::Type{<:TriangularLayout{UPLO,'N'}}, ::TridiagonalToeplitzLayout) where UPLO = BidiagonalToeplitzLayout() +materialize!(L::MatLdivVec{BidiagonalToeplitzLayout,Lay}) where Lay = materialize!(Ldiv{BidiagonalLayout{FillLayout,FillLayout},Lay}(L.A, L.B)) +copyto!(dest::AbstractArray, L::Ldiv{BidiagonalToeplitzLayout,Lay}) where Lay = copyto!(dest, Ldiv{BidiagonalLayout{FillLayout,FillLayout},Lay}(L.A, L.B)) + + +# copy for AdjOrTrans +copy(A::Adjoint{T,<:BandedMatrix{T,<:Any,OneToInf{Int}}}) where T = copy(parent(A))' +copy(A::Transpose{T,<:BandedMatrix{T,<:Any,OneToInf{Int}}}) where T = transpose(copy(parent(A))) + + +## +# hcat +## + +Base.typed_hcat(::Type{T}, A::BandedMatrix{<:Any,<:Any,OneToInf{Int}}, B::AbstractVecOrMat...) where T = Hcat{T}(A, B...) + + + +### +# SymTriPertToeplitz +### + +MemoryLayout(::Type{<:SymTriPertToeplitz}) = PertTridiagonalToeplitzLayout() + + +sublayout(::ApplyBandedLayout, ::Type{<:Tuple{KR,Integer}}) where {KR<:InfAxes} = + sublayout(PaddedColumns{UnknownLayout}(), Tuple{KR}) + +end # module \ No newline at end of file diff --git a/ext/InfiniteArraysBlockArraysExt.jl b/ext/InfiniteArraysBlockArraysExt.jl new file mode 100644 index 0000000..5c2388f --- /dev/null +++ b/ext/InfiniteArraysBlockArraysExt.jl @@ -0,0 +1,10 @@ +module InfiniteArraysBlockArraysExt + +sub_materialize(_, V, ::Tuple{BlockedOneTo{Int,<:InfRanges}}) = V +sub_materialize(::AbstractBlockLayout, V, ::Tuple{BlockedOneTo{Int,<:InfRanges}}) = V +function sub_materialize(::PaddedColumns, v::AbstractVector{T}, ax::Tuple{BlockedOneTo{Int,<:InfRanges}}) where T + dat = paddeddata(v) + BlockedVector(Vcat(sub_materialize(dat), Zeros{T}(∞)), ax) +end + +end # module \ No newline at end of file diff --git a/src/InfiniteArrays.jl b/src/InfiniteArrays.jl index 7a86b57..b3bb8f2 100644 --- a/src/InfiniteArrays.jl +++ b/src/InfiniteArrays.jl @@ -208,5 +208,18 @@ end collect(G::Base.Generator{<:InfRanges}) = BroadcastArray(G.f, G.iter) +# Support PowerBySquaring w/ ∞-arrays +function ArrayLayouts._power_by_squaring(_, ::NTuple{2,InfiniteCardinal{0}}, A::AbstractMatrix{T}, p::Integer) where T + if p < 0 + inv(A)^(-p) + elseif p == 0 + Eye{T}(∞) + elseif p == 1 + copy(A) + else + A*A^(p-1) + end +end + end # module diff --git a/src/infarrays.jl b/src/infarrays.jl index 541b30b..29668a2 100644 --- a/src/infarrays.jl +++ b/src/infarrays.jl @@ -429,3 +429,5 @@ Base.checkindex(::Type{Bool}, inds::AbstractUnitRange, I::AbstractFill) = Base.c LazyArrays.cache_getindex(::InfiniteCardinal{0}, A::AbstractVector, I, J...) = layout_getindex(A, I, J...) LazyArrays.cache_getindex(::InfiniteCardinal{0}, A::CachedVector{<:Any,<:AbstractVector,<:AbstractFill{<:Any,1}}, I::AbstractVector) = LazyArrays.cache_getindex(nothing, A, I) + +*(a::AbstractVector, b::AbstractFill{<:Any,2,Tuple{OneTo{Int},OneToInf{Int}}}) = ApplyArray(*,a,b) \ No newline at end of file diff --git a/src/infrange.jl b/src/infrange.jl index 62610f2..72dca0a 100644 --- a/src/infrange.jl +++ b/src/infrange.jl @@ -618,4 +618,8 @@ end # bounds-checking function Base.checkindex(::Type{Bool}, inds::NTuple{N, AbstractInfUnitRange}, i::AbstractRange{CartesianIndex{N}}) where {N} isempty(i) | checkindex(Bool, inds, first(i)) -end \ No newline at end of file +end + + +LinearAlgebra._cut_B(x::AbstractVector, ::InfUnitRange) = x +LinearAlgebra._cut_B(X::AbstractMatrix, ::InfUnitRange) = X diff --git a/src/reshapedarray.jl b/src/reshapedarray.jl index 2178fd0..989d824 100644 --- a/src/reshapedarray.jl +++ b/src/reshapedarray.jl @@ -2,7 +2,7 @@ using Base.MultiplicativeInverses: SignedMultiplicativeInverse -struct ReshapedArray{T,N,P<:AbstractArray,DIMS<:Tuple,MI<:Tuple{Vararg{SignedMultiplicativeInverse{Int}}}} <: AbstractArray{T,N} +struct ReshapedArray{T,N,P<:AbstractArray,DIMS<:Tuple,MI<:Tuple{Vararg{SignedMultiplicativeInverse{Int}}}} <: LayoutArray{T,N} parent::P dims::DIMS mi::MI diff --git a/test/runtests.jl b/test/runtests.jl index 2a8cfb2..e359cc4 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1254,3 +1254,14 @@ include("test_block.jl") @test checkbounds(Bool, D, diagind(D, IndexCartesian())) end end + +@testset "Vector * ∞ matrix" begin + a = [1+im,2+im] + A = a * Ones{Complex{Int}}(1,∞) + @test A[:,1:5] == a * ones(1,5) + + @test (a*permutedims(1:∞))[:,1:5] == a*(1:5)' + @test (a*Hcat(Zeros(1,2), permutedims(1:∞)))[1,1:5] == (a*Vcat(Hcat(Zeros(1,2), permutedims(1:∞))))[1,1:5] +end + +include("test_infbanded.jl") \ No newline at end of file diff --git a/test/test_infbanded.jl b/test/test_infbanded.jl new file mode 100644 index 0000000..dd26622 --- /dev/null +++ b/test/test_infbanded.jl @@ -0,0 +1,348 @@ +using ArrayLayouts, InfiniteArrays, BandedMatrices, FillArrays, LazyArrays, Test +import BandedMatrices: _BandedMatrix, bandeddata + +const InfiniteArraysBandedMatricesExt = Base.get_extension(InfiniteArrays, :InfiniteArraysBandedMatricesExt) +const InfToeplitz = InfiniteArraysBandedMatricesExt.InfToeplitz +const TriToeplitz = InfiniteArraysBandedMatricesExt.TriToeplitz +const SymTriPertToeplitz = InfiniteArraysBandedMatricesExt.SymTriPertToeplitz +const TriPertToeplitz = InfiniteArraysBandedMatricesExt.TriPertToeplitz +const AdjTriPertToeplitz = InfiniteArraysBandedMatricesExt.AdjTriPertToeplitz +const ConstRows = InfiniteArraysBandedMatricesExt.ConstRows +const BandedToeplitzLayout = InfiniteArraysBandedMatricesExt.BandedToeplitzLayout +const TridiagonalToeplitzLayout = InfiniteArraysBandedMatricesExt.TridiagonalToeplitzLayout +const BidiagonalToeplitzLayout = InfiniteArraysBandedMatricesExt.BidiagonalToeplitzLayout +const PertToeplitz = InfiniteArraysBandedMatricesExt.PertToeplitz +const PertToeplitzLayout = InfiniteArraysBandedMatricesExt.PertToeplitzLayout +const InfBandCartesianIndices = InfiniteArraysBandedMatricesExt.InfBandCartesianIndices +const subdiagonalconstant = InfiniteArraysBandedMatricesExt.subdiagonalconstant +const diagonalconstant = InfiniteArraysBandedMatricesExt.diagonalconstant +const supdiagonalconstant = InfiniteArraysBandedMatricesExt.supdiagonalconstant + +using Base: oneto +using LazyArrays: simplifiable, ApplyLayout, BroadcastBandedLayout, islazy + +@testset "∞-banded" begin + @testset "Diagonal and BandedMatrix" begin + D = Diagonal(Fill(2,∞)) + + B = D[1:∞,2:∞] + @test B isa BandedMatrix + @test B[1:10,1:10] == diagm(-1 => Fill(2,9)) + @test B[1:∞,2:∞] isa BandedMatrix + + A = BandedMatrix(0 => 1:∞, 1=> Fill(2.0,∞), -1 => Fill(3.0,∞)) + x = [1; 2; zeros(∞)] + @test A*x isa Vcat + @test (A*x)[1:10] == A[1:10,1:10]*x[1:10] + + @test InfBandCartesianIndices(0)[1:5] == CartesianIndex.(1:5,1:5) + @test InfBandCartesianIndices(1)[1:5] == CartesianIndex.(1:5,2:6) + @test InfBandCartesianIndices(-1)[1:5] == CartesianIndex.(2:6,1:5) + + @test A[band(0)][2:10] == 2:10 + @test D[band(0)] ≡ Fill(2,∞) + @test D[band(1)] ≡ Fill(0,∞) + + @test B*A*x isa Vcat + @test (B*A*x)[1:10] == [0; 10; 14; 12; zeros(6)] + + @test _BandedMatrix((1:∞)', ∞, -1, 1) isa BandedMatrix + end + + @testset "∞-Toeplitz" begin + A = BandedMatrix(1 => Fill(2im,∞), 2 => Fill(-1,∞), 3 => Fill(2,∞), -2 => Fill(-4,∞), -3 => Fill(-2im,∞)) + @test A isa InfToeplitz + @test MemoryLayout(A.data) == ConstRows() + @test MemoryLayout(A) == BandedToeplitzLayout() + @test islazy(A) == Val(true) + + V = view(A,:,3:∞) + @test MemoryLayout(typeof(bandeddata(V))) == ConstRows() + @test MemoryLayout(typeof(V)) == BandedToeplitzLayout() + @test BandedMatrix(V) isa InfToeplitz + @test A[:,3:end] isa InfToeplitz + + @test (A + 2I)[1:10,1:10] == (2I + A)[1:10,1:10] == A[1:10,1:10] + 2I + @test (A*A)[1:10,1:10] ≈ A[1:10,1:16]A[1:16,1:10] + @test (A * Fill(2,∞))[1:10] ≈ 2A[1:10,1:16]*ones(16) + @test (Fill(2,∞,∞)*A)[1:10,1:10] ≈ fill(2,10,13)A[1:13,1:10] + + @test Eye(∞) * A isa BandedMatrix + @test A * Eye(∞) isa BandedMatrix + + @test A * [1; 2; Zeros(∞)] isa Vcat + @test A * [1; 2; Zeros(∞)] == [A[1:5,1:2] * [1,2]; Zeros(∞)] + + @test A * Vcat([1 2; 3 4], Zeros(∞,2)) isa ApplyMatrix{ComplexF64,typeof(Base.setindex)} + @test simplifiable(*, A, Vcat([1 2; 3 4], Zeros(∞,2))) == Val(true) + + @test MemoryLayout(Tridiagonal(Fill(1,∞), Fill(2,∞), Fill(3,∞))) isa TridiagonalToeplitzLayout + @test MemoryLayout(Bidiagonal(Fill(1,∞), Fill(2,∞), :U)) isa BidiagonalToeplitzLayout + @test MemoryLayout(SymTridiagonal(Fill(1,∞), Fill(2,∞))) isa TridiagonalToeplitzLayout + + + @testset "algebra" begin + T = Tridiagonal(Fill(1,∞), Fill(2,∞), Fill(3,∞)) + @test T isa TriToeplitz + @test (T + 2I)[1:10,1:10] == (2I + T)[1:10,1:10] == T[1:10,1:10] + 2I + @test BandedMatrix(T, (2,3))[1:10,1:10] == InfToeplitz(T,(3,2))[1:10,1:10] == InfToeplitz(T)[1:10,1:10] == T[1:10,1:10] + @test subdiagonalconstant(T) == 1 + @test diagonalconstant(T) == 2 + @test supdiagonalconstant(T) == 3 + + S = SymTridiagonal(Fill(1,∞), Fill(2,∞)) + @test (S + 2I)[1:10,1:10] == (2I + S)[1:10,1:10] == S[1:10,1:10] + 2I + @test BandedMatrix(S, (2,3))[1:10,1:10] == S[1:10,1:10] + end + + @testset "constant data" begin + A = BandedMatrix(1 => Fill(2im,∞), 2 => Fill(-1,∞), 3 => Fill(2,∞), -2 => Fill(-4,∞), -3 => Fill(-2im,∞)) + B = _BandedMatrix(Fill(2,4,∞), ∞, 1,2) + @test (B*B)[1:10,1:10] ≈ B[1:10,1:14]B[1:14,1:10] + @test (A*B)[1:10,1:10] ≈ A[1:10,1:14]B[1:14,1:10] + @test (B*A)[1:10,1:10] ≈ B[1:10,1:14]A[1:14,1:10] + @test simplifiable(*, B, B) == Val(true) + @test length((A*B*B).args) == 2 + @test length((B*B*A).args) == 2 + end + + @testset "Toep * Diag" begin + A = BandedMatrix(1 => Fill(2im,∞), 2 => Fill(-1,∞), 3 => Fill(2,∞), -2 => Fill(-4,∞), -3 => Fill(-2im,∞)) + D = Diagonal(1:∞) + @test D*A isa BroadcastMatrix + @test A*D isa BroadcastMatrix + @test simplifiable(*, D, A) == Val(true) + @test simplifiable(*, A, D) == Val(true) + end + + @testset "change bands" begin + A = SymTridiagonal(Fill(2,∞), Fill(1,∞)) + @test BandedMatrix(A, (3,1))[1:5,1:5] == A[1:5,1:5] + end + end + + @testset "Pert-Toeplitz" begin + @testset "Inf Pert" begin + A = BandedMatrix(-2 => Vcat(Float64[], Fill(1/4,∞)), 0 => Vcat([1.0+im,2,3],Fill(0,∞)), 1 => Vcat(Float64[], Fill(1,∞))) + @test A isa PertToeplitz + @test MemoryLayout(A) isa PertToeplitzLayout + + V = view(A,2:∞,2:∞) + @test MemoryLayout(V) isa PertToeplitzLayout + @test BandedMatrix(V) isa PertToeplitz + @test islazy(V) == Val(true) + @test A[2:∞,2:∞] isa PertToeplitz + + @test (A + 2I)[1:10,1:10] == (2I + A)[1:10,1:10] == A[1:10,1:10] + 2I + + @test Eye(∞) * A isa BandedMatrix + @test A * Eye(∞) isa BandedMatrix + end + + @testset "TriPert" begin + A = SymTridiagonal(Vcat([1,2.], Fill(2.,∞)), Vcat([3.,4.], Fill.(0.5,∞))) + @test A isa SymTriPertToeplitz + @test (A + 2I)[1:10,1:10] == (2I + A)[1:10,1:10] == A[1:10,1:10] + 2I + @test BandedMatrix(A, (2,3))[1:10,1:10] == A[1:10,1:10] + + A = Tridiagonal(Vcat([3.,4.], Fill.(0.5,∞)), Vcat([1,2.], Fill(2.,∞)), Vcat([3.,4.], Fill.(0.5,∞))) + @test A isa TriPertToeplitz + @test Adjoint(A) isa AdjTriPertToeplitz + @test (A + 2I)[1:10,1:10] == (2I + A)[1:10,1:10] == A[1:10,1:10] + 2I + @test (Adjoint(A) + 2I)[1:10,1:10] == (2I + Adjoint(A))[1:10,1:10] == Adjoint(A)[1:10,1:10] + 2I + @test BandedMatrix(A, (2,3))[1:10,1:10] == A[1:10,1:10] + end + + + + @testset "InfBanded" begin + A = _BandedMatrix(Fill(2,4,∞),ℵ₀,2,1) + B = _BandedMatrix(Fill(3,2,∞),ℵ₀,-1,2) + @test mul(A,A) isa PertToeplitz + @test A*A isa PertToeplitz + @test (A*A)[1:20,1:20] == A[1:20,1:23]*A[1:23,1:20] + @test (A*B)[1:20,1:20] == A[1:20,1:23]*B[1:23,1:20] + end + + @testset "change bands" begin + A = BandedMatrix(-2 => Vcat(Float64[], Fill(1/4,∞)), 0 => Vcat([1.0+im,2,3],Fill(0,∞)), 1 => Vcat(Float64[], Fill(1,∞))) + @test BandedMatrix(A, (3,1))[1:5,1:5] == A[1:5,1:5] + end + end + + @testset "adjortrans" begin + A = BandedMatrix(0 => 1:∞, 1=> Fill(2.0+im,∞), -1 => Fill(3.0,∞)) + @test copy(A')[1:10,1:10] == (A')[1:10,1:10] + @test copy(transpose(A))[1:10,1:10] == transpose(A)[1:10,1:10] + end + + @testset "Eye subindex" begin + @test Eye(∞)[:,1:3][1:5,:] == Eye(∞)[Base.Slice(oneto(∞)),1:3][1:5,:] == Eye(5,3) + @test Eye(∞)[1:3,:][:,1:5] == Eye(∞)[1:3,Base.Slice(oneto(∞))][:,1:5] == Eye(3,5) + @test Eye(∞)[:,:][1:5,1:3] == Eye(∞)[Base.Slice(oneto(∞)),Base.Slice(oneto(∞))][1:5,1:3] == Eye(5,3) + end + + @testset "band(0) indexing" begin + D = ApplyArray(*, Diagonal(1:∞), Diagonal(1:∞)) + @test D[band(0)][1:10] == (1:10).^2 + end + + @testset "Fill * Banded" begin + A = _BandedMatrix(Ones(1,∞), ∞, 1,-1) + B = _BandedMatrix(Fill(1.0π,1,∞), ∞, 0,0) + @test (A*B)[1:10,1:10] ≈ BandedMatrix(-1 => Fill(1.0π,9)) + end + + @testset "concat" begin + H = ApplyArray(hvcat, 2, 1, [1 Zeros(1,∞)], [1; Zeros(∞)], Diagonal(1:∞)) + @test bandwidths(H) == (1,1) + H = ApplyArray(hvcat, 2, 1, [0 Zeros(1,∞)], [0; Zeros(∞)], Diagonal(1:∞)) + @test bandwidths(H) == (0,0) + H = ApplyArray(hvcat, (2,2), 1, [1 Zeros(1,∞)], [1; Zeros(∞)], Diagonal(1:∞)) + @test_broken bandwidths(H) == (1,1) + end + + @testset "Banded * PaddedMatrix" begin + A = Eye(∞)[2:∞,:] + @test LazyArrays.islazy(A) == Val(true) + B = PaddedArray(randn(3,3),ℵ₀,ℵ₀) + @test (A*B)[1:10,1:10] ≈ A[1:10,1:10] * B[1:10,1:10] + end + + @testset "SubArray broadcasting" begin + A = BandedMatrix(2 => 1:∞) + @test exp.(A[1:2:∞,1:2:∞])[1:10,1:10] ≈ exp.(A[1:2:20,1:2:20]) + @test A[band(2)][1:5] == 1:5 + @test _BandedMatrix((1:∞)', ∞, -1,1)[band(1)][1:5] == 2:6 + @test exp.(view(A,band(2)))[1:10] ≈ exp.(1:10) + + @test BandedMatrices.banded_similar(Int, (∞,5), (1,1)) isa BandedMatrix + @test BandedMatrices.banded_similar(Int, (5,∞), (1,1)) isa Adjoint{<:Any,<:BandedMatrix} + + A = BandedMatrix{Int}((2 => 1:∞,), (∞,∞), (0,2)) + @test eltype(A) == Int + @test bandwidths(A) == (0,2) + + A = BandedMatrix{Int}((2 => Vcat([1,2], Fill(2,∞)),), (∞,∞), (0,2)) + @test A[band(2)][1:5] == [1; fill(2,4)] + @test A[band(4)] ≡ Zeros{Int}(∞) + end + + @testset "prepad with fill" begin + A = BandedMatrix(2 => Zeros(∞)) + @test A[band(2)][1:10] == Zeros(10) + A = BandedMatrix(2 => Ones(∞)) + @test A[band(2)][1:10] == Ones(10) + end + + @testset "Algebra" begin + A = BandedMatrix(-3 => Fill(7 / 10, ∞), -2 => 1:∞, 1 => Fill(2im, ∞)) + @test A isa BandedMatrix{ComplexF64} + @test A[1:10, 1:10] == diagm(-3 => Fill(7 / 10, 7), -2 => 1:8, 1 => Fill(2im, 9)) + + A = BandedMatrix(0 => Vcat([1, 2, 3], Zeros(∞)), 1 => Vcat(1, Zeros(∞))) + @test A[1, 2] == 1 + + A = BandedMatrix(-3 => Fill(7 / 10, ∞), -2 => Fill(1, ∞), 1 => Fill(2im, ∞)) + Ac = BandedMatrix(A') + At = BandedMatrix(transpose(A)) + @test Ac[1:10, 1:10] ≈ (A')[1:10, 1:10] ≈ A[1:10, 1:10]' + @test At[1:10, 1:10] ≈ transpose(A)[1:10, 1:10] ≈ transpose(A[1:10, 1:10]) + + A = BandedMatrix(-1 => Vcat(Float64[], Fill(1 / 4, ∞)), 0 => Vcat([1.0 + im], Fill(0, ∞)), 1 => Vcat(Float64[], Fill(1, ∞))) + @test MemoryLayout(typeof(view(A.data, :, 1:10))) == ApplyLayout{typeof(hcat)}() + Ac = BandedMatrix(A') + At = BandedMatrix(transpose(A)) + @test Ac[1:10, 1:10] ≈ (A')[1:10, 1:10] ≈ A[1:10, 1:10]' + @test At[1:10, 1:10] ≈ transpose(A)[1:10, 1:10] ≈ transpose(A[1:10, 1:10]) + + A = BandedMatrix(-2 => Vcat(Float64[], Fill(1 / 4, ∞)), 0 => Vcat([1.0 + im, 2, 3], Fill(0, ∞)), 1 => Vcat(Float64[], Fill(1, ∞))) + Ac = BandedMatrix(A') + At = BandedMatrix(transpose(A)) + @test Ac[1:10, 1:10] ≈ (A')[1:10, 1:10] ≈ A[1:10, 1:10]' + @test At[1:10, 1:10] ≈ transpose(A)[1:10, 1:10] ≈ transpose(A[1:10, 1:10]) + + A = _BandedMatrix(Fill(1, 4, ∞), ℵ₀, 1, 2) + @test A^2 isa BandedMatrix + @test (A^2)[1:10, 1:10] == (A*A)[1:10, 1:10] == (A[1:100, 1:100]^2)[1:10, 1:10] + @test A^3 isa ApplyMatrix{<:Any,typeof(*)} + @test (A^3)[1:10, 1:10] == (A*A*A)[1:10, 1:10] == ((A*A)*A)[1:10, 1:10] == (A*(A*A))[1:10, 1:10] == (A[1:100, 1:100]^3)[1:10, 1:10] + + @testset "∞ x finite" begin + A = BandedMatrix(1 => 1:∞) + BandedMatrix(-1 => Fill(2, ∞)) + B = _BandedMatrix(randn(3, 5), ℵ₀, 1, 1) + + @test lmul!(2.0, copy(B)')[:, 1:10] == (2B')[:, 1:10] + + @test_throws ArgumentError BandedMatrix(A) + @test A * B isa MulMatrix + @test B'A isa MulMatrix + + @test all(diag(A[1:6, 1:6]) .=== zeros(Int,6)) + + @test (A*B)[1:7, 1:5] ≈ A[1:7, 1:6] * B[1:6, 1:5] + @test (B'A)[1:5, 1:7] ≈ (B')[1:5, 1:6] * A[1:6, 1:7] + end + end + + @testset "Fill" begin + A = _BandedMatrix(Ones(1, ∞), ℵ₀, -1, 1) + @test 1.0 .* A isa BandedMatrix{Float64,<:Fill} + @test Zeros(∞) .* A ≡ Zeros(∞, ∞) .* A ≡ A .* Zeros(1, ∞) ≡ A .* Zeros(∞, ∞) ≡ Zeros(∞, ∞) + @test Ones(∞) .* A isa BandedMatrix{Float64,<:Ones} + @test A .* Ones(1, ∞) isa BandedMatrix{Float64,<:Ones} + @test 2.0 .* A isa BandedMatrix{Float64,<:Fill} + @test A .* 2.0 isa BandedMatrix{Float64,<:Fill} + @test Eye(∞) * A isa BandedMatrix{Float64,<:Ones} + @test A * Eye(∞) isa BandedMatrix{Float64,<:Ones} + + @test A * A isa BandedMatrix + @test (A*A)[1:10, 1:10] == BandedMatrix(2 => Ones(8)) + + à = _BandedMatrix(Fill(1, 1, ∞), ℵ₀, -1, 1) + @test A * à isa BandedMatrix + @test à * A isa BandedMatrix + @test à * à isa BandedMatrix + + B = _BandedMatrix(Ones(1, 10), ℵ₀, -1, 1) + C = _BandedMatrix(Ones(1, 10), 10, -1, 1) + D = _BandedMatrix(Ones(1, ∞), 10, -1, 1) + + @test (A*B)[1:10, 1:10] == (B*C)[1:10, 1:10] == (D*A)[1:10, 1:10] == D * B == (C*D)[1:10, 1:10] == BandedMatrix(2 => Ones(8)) + end + + @testset "Banded Broadcast" begin + A = _BandedMatrix((1:∞)', ℵ₀, -1, 1) + @test 2.0 .* A isa BandedMatrix{Float64,<:Adjoint} + @test A .* 2.0 isa BandedMatrix{Float64,<:Adjoint} + @test Eye(∞) * A isa BandedMatrix{Float64,<:Adjoint} + @test A * Eye(∞) isa BandedMatrix{Float64,<:Adjoint} + A = _BandedMatrix(Vcat((1:∞)', Ones(1, ∞)), ℵ₀, 0, 1) + @test 2.0 .* A isa BandedMatrix + @test A .* 2.0 isa BandedMatrix + @test Eye(∞) * A isa BandedMatrix + @test A * Eye(∞) isa BandedMatrix + b = 1:∞ + @test bandwidths(b .* A) == (0, 1) + + @test colsupport(b .* A, 1) == 1:1 + @test Base.replace_in_print_matrix(b .* A, 2, 1, "0.0") == " ⋅ " + @test bandwidths(A .* b) == (0, 1) + @test A .* b' isa BroadcastArray + @test bandwidths(A .* b') == bandwidths(A .* b') + @test colsupport(A .* b', 3) == 2:3 + + A = _BandedMatrix(Ones{Int}(1, ∞), ℵ₀, 0, 0)' + B = _BandedMatrix((-2:-2:-∞)', ℵ₀, -1, 1) + C = Diagonal(2 ./ (1:2:∞)) + @test bandwidths(A * (B * C)) == (-1, 1) + @test bandwidths((A * B) * C) == (-1, 1) + + A = _BandedMatrix(Ones{Int}(1, ∞), ℵ₀, 0, 0)' + B = _BandedMatrix((-2:-2:-∞)', ℵ₀, -1, 1) + @test MemoryLayout(A + B) isa BroadcastBandedLayout{typeof(+)} + @test MemoryLayout(2 * (A + B)) isa BroadcastBandedLayout{typeof(*)} + @test bandwidths(A + B) == (0, 1) + @test bandwidths(2 * (A + B)) == (0, 1) + end +end