Skip to content

Commit

Permalink
Move Tridiagonal Toeplitz code out of extension (#201)
Browse files Browse the repository at this point in the history
* Move Tridiagonal Toeplitz code out of extension

* tests pass

* move ConstRows

* Update InfiniteArrays.jl

* increase cov

* Update test_infbanded.jl
  • Loading branch information
dlfivefifty authored Nov 30, 2024
1 parent 210d087 commit 1a4f215
Show file tree
Hide file tree
Showing 5 changed files with 151 additions and 125 deletions.
2 changes: 1 addition & 1 deletion Project.toml
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
name = "InfiniteArrays"
uuid = "4858937d-0d70-526a-a4dd-2d5cb5dd786c"
version = "0.15.0"
version = "0.15.0-dev"

[deps]
ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
Expand Down
110 changes: 7 additions & 103 deletions ext/InfiniteArraysBandedMatricesExt.jl
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@ using InfiniteArrays.LazyArrays, InfiniteArrays.ArrayLayouts, InfiniteArrays.Fil

import Base: BroadcastStyle, size, getindex, similar, copy, *, +, -, /, \, materialize!, copyto!, OneTo
import Base.Broadcast: Broadcasted
import InfiniteArrays: InfIndexRanges, Infinity, PosInfinity, OneToInf, InfAxes, AbstractInfUnitRange, InfRanges
import InfiniteArrays: InfIndexRanges, Infinity, PosInfinity, OneToInf, InfAxes, AbstractInfUnitRange, InfRanges, InfBaseToeplitzLayouts, ConstRowMatrix, PertConstRowMatrix, SymTriPertToeplitz, TriPertToeplitz, ConstRows, PertConstRows
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
Expand Down Expand Up @@ -51,17 +51,10 @@ sub_materialize(_, V::SubArray{<:Any,1,<:AbstractMatrix,Tuple{InfBandCartesianIn




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}}
const InfToeplitz{T} = InfBandedMatrix{T,<:ConstRowMatrix{T}}
const PertToeplitz{T} = InfBandedMatrix{T,<:PertConstRowMatrix{T}}


_prepad(p, a) = Vcat(Zeros{eltype(a)}(max(p,0)), a)
_prepad(p, a::Zeros{T,1}) where T = Zeros{T}(length(a)+p)
Expand Down Expand Up @@ -157,50 +150,6 @@ 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)
Expand Down Expand Up @@ -339,10 +288,6 @@ ConstRowMatrix(A::AbstractMatrix{T}) where T = ApplyMatrix(*, A[:,1], Ones{T}(1,
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...)
Expand All @@ -355,29 +300,17 @@ for Typ in (:ConstRows, :PertConstRows)
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}
const InfBandedToeplitzLayouts = Union{BandedToeplitzLayout, PertToeplitzLayout, PertTriangularToeplitzLayout}
const InfToeplitzLayouts = Union{InfBaseToeplitzLayouts, InfBandedToeplitzLayouts}

subdiagonalconstant(A) = getindex_value(subdiagonaldata(A))
diagonalconstant(A) = getindex_value(diagonaldata(A))
supdiagonalconstant(A) = getindex_value(supdiagonaldata(A))


islazy_layout(::InfToeplitzLayouts) = Val(true)
islazy_layout(::InfBandedToeplitzLayouts) = Val(true)
islazy(::BandedMatrix{<:Any,<:Any,OneToInf{Int}}) = Val(true)


Expand All @@ -399,8 +332,6 @@ _BandedMatrix(::PertToeplitzLayout, A::AbstractMatrix) =
# 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)

Expand Down Expand Up @@ -473,33 +404,6 @@ 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)))
Expand Down
17 changes: 9 additions & 8 deletions src/InfiniteArrays.jl
Original file line number Diff line number Diff line change
Expand Up @@ -14,10 +14,7 @@ import Base: *, +, -, /, <, ==, >, \, ≤, ≥, (:), @propagate_inbounds,
searchsortedfirst, searchsortedlast, setindex!, show, show_circular, show_delim_array, sign,
signbit, similar, size, sort, sort!, step, sum, tail,
to_shape, transpose, unaliascopy, union, unitrange_last, unsafe_convert, unsafe_indices, unsafe_length,
vcat, zeros


import Base: range_start_step_length
vcat, zeros, copyto!, range_start_step_length

if VERSION v"1.11.0-DEV.21"
using LinearAlgebra: UpperOrLowerTriangular
Expand All @@ -31,16 +28,18 @@ end

using Base.Broadcast
import ArrayLayouts: AbstractBandedLayout, LayoutMatrix, LayoutVecOrMat, LayoutVecOrMats, LayoutVector, MemoryLayout,
RangeCumsum, UnknownLayout, reshapedlayout, sub_materialize, sublayout
RangeCumsum, UnknownLayout, reshapedlayout, sub_materialize, materialize!, sublayout, MatLdivVec,
subdiagonaldata, diagonaldata, supdiagonaldata, triangularlayout

import Base.Broadcast: BroadcastStyle, Broadcasted, DefaultArrayStyle, axistype, broadcasted

import FillArrays: AbstractFill, Eye, Fill, Ones, RectDiagonal, Zeros, fill_reshape, getindex_value
import FillArrays: AbstractFill, Eye, Fill, Ones, RectDiagonal, Zeros, fill_reshape, getindex_value, AbstractFillMatrix

import Infinities: InfiniteCardinal, Infinity, ∞

import LazyArrays: AbstractCachedVector, ApplyLayout, CachedArray, CachedVector, InvColumnLayout,
LazyArrayStyle, LazyLayout, LazyMatrix, PaddedColumns, _padded_sub_materialize, sub_paddeddata
import LazyArrays: AbstractLazyLayout, AbstractCachedVector, ApplyLayout, CachedArray, CachedVector, InvColumnLayout, AbstractLazyBandedLayout,
LazyArrayStyle, LazyLayout, LazyMatrix, PaddedColumns, _padded_sub_materialize, sub_paddeddata,
ApplyBandedLayout, BroadcastBandedLayout, islazy_layout

import LinearAlgebra: AdjOrTrans, HermOrSym, diag, norm, norm1, norm2, normp

Expand Down Expand Up @@ -221,5 +220,7 @@ function ArrayLayouts._power_by_squaring(_, ::NTuple{2,InfiniteCardinal{0}}, A::
end
end

include("inftoeplitz.jl")


end # module
116 changes: 116 additions & 0 deletions src/inftoeplitz.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,116 @@
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}}}

struct ConstRows <: AbstractLazyLayout end
struct PertConstRows <: AbstractLazyLayout end
MemoryLayout(::Type{<:ConstRowMatrix}) = ConstRows()
MemoryLayout(::Type{<:PertConstRowMatrix}) = PertConstRows()


const TriToeplitz{T} = Tridiagonal{T,Fill{T,1,Tuple{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}}}}}}}


"""
TridiagonalToeplitzLayout
represents a matrix which is tridiagonal and toeplitz. Must support
`subdiagonalconstant`, `diagonalconstant`, `supdiagonalconstant`.
"""
struct TridiagonalToeplitzLayout <: AbstractLazyBandedLayout end

struct BidiagonalToeplitzLayout <: AbstractLazyBandedLayout end
struct PertBidiagonalToeplitzLayout <: AbstractLazyBandedLayout end
struct PertTridiagonalToeplitzLayout <: AbstractLazyBandedLayout end

const InfBaseToeplitzLayouts = Union{TridiagonalToeplitzLayout, BidiagonalToeplitzLayout, PertBidiagonalToeplitzLayout, PertTridiagonalToeplitzLayout}


subdiagonalconstant(A) = getindex_value(subdiagonaldata(A))
diagonalconstant(A) = getindex_value(diagonaldata(A))
supdiagonalconstant(A) = getindex_value(supdiagonaldata(A))

islazy_layout(::InfBaseToeplitzLayouts) = Val(true)

@inline sub_materialize(::ApplyBandedLayout{typeof(*)}, V, ::Tuple{InfAxes,InfAxes}) = V
@inline sub_materialize(::BroadcastBandedLayout, V, ::Tuple{InfAxes,InfAxes}) = V


###
# 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))



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
end
end
Loading

0 comments on commit 1a4f215

Please sign in to comment.