Skip to content

Move Tridiagonal Toeplitz code out of extension #201

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 6 commits into from
Nov 30, 2024
Merged

Move Tridiagonal Toeplitz code out of extension #201

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
fa605c3
Move Tridiagonal Toeplitz code out of extension
dlfivefifty Nov 30, 2024
6c4de6e
tests pass
dlfivefifty Nov 30, 2024
5cda2be
move ConstRows
dlfivefifty Nov 30, 2024
dfe2bdb
Update InfiniteArrays.jl
dlfivefifty Nov 30, 2024
4ff3eea
increase cov
dlfivefifty Nov 30, 2024
f73173e
Update test_infbanded.jl
dlfivefifty Nov 30, 2024
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
2 changes: 1 addition & 1 deletion Project.toml
View file
Open in desktop
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
View file
Open in desktop
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
View file
Open in desktop
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
View file
Open in desktop
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