From bb5a8747b109c4c42aae26470a0f4755511eed93 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Tue, 19 Nov 2024 09:34:36 +0000 Subject: [PATCH 1/5] Bump codecov/codecov-action from 4 to 5 (#189) Bumps [codecov/codecov-action](https://github.com/codecov/codecov-action) from 4 to 5. - [Release notes](https://github.com/codecov/codecov-action/releases) - [Changelog](https://github.com/codecov/codecov-action/blob/main/CHANGELOG.md) - [Commits](https://github.com/codecov/codecov-action/compare/v4...v5) --- updated-dependencies: - dependency-name: codecov/codecov-action dependency-type: direct:production update-type: version-update:semver-major ... Signed-off-by: dependabot[bot] Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> --- .github/workflows/ci.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index af919fc..34d7e40 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -63,7 +63,7 @@ jobs: - uses: julia-actions/julia-buildpkg@v1 - uses: julia-actions/julia-runtest@v1 - uses: julia-actions/julia-processcoverage@v1 - - uses: codecov/codecov-action@v4 + - uses: codecov/codecov-action@v5 with: token: ${{ secrets.CODECOV_TOKEN }} file: lcov.info From ad4af83ec549f78098ef770b611eb6961245367f Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Wed, 27 Nov 2024 15:23:20 +0000 Subject: [PATCH 2/5] overload simplifiable for special cases of Mul (#188) * overload simplifiable for special cases of Mul * increase codecov --- Project.toml | 2 +- src/banded/infbanded.jl | 3 +++ test/test_infbanded.jl | 16 ++++++++++++++++ 3 files changed, 20 insertions(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 7d88418..6f939df 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "InfiniteLinearAlgebra" uuid = "cde9dba0-b1de-11e9-2c62-0bab9446c55c" -version = "0.8.4" +version = "0.8.5" [deps] ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a" diff --git a/src/banded/infbanded.jl b/src/banded/infbanded.jl index 215104a..67121ce 100644 --- a/src/banded/infbanded.jl +++ b/src/banded/infbanded.jl @@ -428,6 +428,7 @@ similar(M::MulAdd{<:AbstractBandedLayout,<:AbstractBandedLayout}, ::Type{T}, axe # 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)) @@ -456,6 +457,8 @@ 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) diff --git a/test/test_infbanded.jl b/test/test_infbanded.jl index 804e17e..5a383e0 100644 --- a/test/test_infbanded.jl +++ b/test/test_infbanded.jl @@ -2,6 +2,7 @@ using InfiniteLinearAlgebra, ArrayLayouts, InfiniteArrays, BandedMatrices, FillA import BandedMatrices: _BandedMatrix, bandeddata using InfiniteLinearAlgebra: InfToeplitz, ConstRows, BandedToeplitzLayout, TridiagonalToeplitzLayout, BidiagonalToeplitzLayout, PertToeplitz, PertToeplitzLayout using Base: oneto +using LazyArrays: simplifiable @testset "∞-banded" begin @testset "Diagonal and BandedMatrix" begin @@ -55,6 +56,9 @@ using Base: oneto @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 @@ -87,6 +91,18 @@ using Base: oneto @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 end From 24f9718beb734300de033441fc690c5a7fdb5b04 Mon Sep 17 00:00:00 2001 From: Daniel VandenHeuvel <95613936+DanielVandH@users.noreply.github.com> Date: Wed, 27 Nov 2024 15:23:49 +0000 Subject: [PATCH 3/5] `findblock` for `RealInfinity()` (#190) * findblock for realinfinity * ci --- .github/workflows/ci.yml | 3 ++- src/InfiniteLinearAlgebra.jl | 2 +- src/blockbanded/blockbanded.jl | 1 + test/runtests.jl | 6 +++++- 4 files changed, 9 insertions(+), 3 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 34d7e40..e2bb954 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -37,7 +37,8 @@ jobs: fail-fast: false matrix: version: - - '1.10' + - '1' + - 'lts' os: - ubuntu-latest - macOS-latest diff --git a/src/InfiniteLinearAlgebra.jl b/src/InfiniteLinearAlgebra.jl index 175f1c2..0003bba 100644 --- a/src/InfiniteLinearAlgebra.jl +++ b/src/InfiniteLinearAlgebra.jl @@ -34,7 +34,7 @@ import FillArrays: AbstractFill, AbstractFillMatrix, axes_print_matrix_row, geti import InfiniteArrays: AbstractInfUnitRange, InfAxes, InfRanges, InfStepRange, InfUnitRange, OneToInf, PosInfinity, InfIndexRanges -import Infinities: InfiniteCardinal, Infinity +import Infinities: InfiniteCardinal, Infinity, RealInfinity import LazyArrays: AbstractCachedMatrix, AbstractCachedVector, AbstractLazyLayout, ApplyArray, ApplyLayout, ApplyMatrix, CachedArray, CachedLayout, CachedMatrix, CachedVector, LazyArrayStyle, LazyLayout, diff --git a/src/blockbanded/blockbanded.jl b/src/blockbanded/blockbanded.jl index 68383d2..619a154 100644 --- a/src/blockbanded/blockbanded.jl +++ b/src/blockbanded/blockbanded.jl @@ -17,6 +17,7 @@ BlockBandedMatrices.blockbanded_rowstop(A, x::InfiniteCardinal{0}) = x BlockArrays.blocklasts(a::InfRanges) = Fill(length(a),1) +BlockArrays.findblock(::BlockedOneTo, ::RealInfinity) = Block(ℵ₀) function BlockArrays.sortedunion(a::Vcat{Int,1,<:Tuple{Union{Int,AbstractVector{Int}},<:AbstractRange}}, b::Vcat{Int,1,<:Tuple{Union{Int,AbstractVector{Int}},<:AbstractRange}}) diff --git a/test/runtests.jl b/test/runtests.jl index 4481369..c76adc4 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -6,7 +6,7 @@ import InfiniteLinearAlgebra: qltail, toeptail, tailiterate, tailiterate!, tail_ BandedToeplitzLayout, PertToeplitzLayout, TridiagonalToeplitzLayout, BidiagonalToeplitzLayout, BidiagonalConjugation import Base: BroadcastStyle, oneto -import BlockArrays: _BlockArray, blockcolsupport +import BlockArrays: _BlockArray, blockcolsupport, findblock import BlockBandedMatrices: isblockbanded, _BlockBandedMatrix import MatrixFactorizations: QLPackedQ import BandedMatrices: bandeddata, _BandedMatrix, BandedStyle @@ -472,6 +472,10 @@ end end end +@testset "findblock at +∞, HarmonicOrthogonalPolynomials#88" begin + @test findblock(blockedrange(1:2:∞), RealInfinity()) == Block(ℵ₀) +end + include("test_hessenbergq.jl") include("test_infql.jl") include("test_infqr.jl") From 48f4e0ad4eb82127b55b584d3c13e1df3a3f224c Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Wed, 27 Nov 2024 16:38:29 +0000 Subject: [PATCH 4/5] Drop Julia Date: Thu, 5 Dec 2024 23:40:54 +0000 Subject: [PATCH 5/5] Move out BandedMatrices code (#192) * Drop Julia 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) - -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 - -## -# 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) - -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) - - -function _bidiag_forwardsub!(M::Ldiv{<:Any,<:PaddedColumns}) - A, b_in = M.A, M.B - dv = diagonaldata(A) - ev = subdiagonaldata(A) - b = paddeddata(b_in) - N = length(b) - b[1] = bj1 = dv[1]\b[1] - @inbounds for j = 2:N - bj = b[j] - bj -= ev[j - 1] * bj1 - dvj = dv[j] - if iszero(dvj) - throw(SingularEbception(j)) - end - bj = dvj\bj - b[j] = bj1 = bj - end - - @inbounds for j = N+1:length(b_in) - iszero(bj1) && break - bj = -ev[j - 1] * bj1 - dvj = dv[j] - if iszero(dvj) - throw(SingularEbception(j)) - end - bj = dvj\bj - b_in[j] = bj1 = bj - end - - b_in -end - -### -# Inf-Toeplitz layout -# this could possibly be avoided via an InfFillLayout -### - -const InfFill = AbstractFill{<:Any,1,<:Tuple{OneToInf}} - -for Typ in (:(LinearAlgebra.Tridiagonal{<:Any,<:InfFill}), - :(LinearAlgebra.SymTridiagonal{<:Any,<:InfFill}), - :(LazyBandedMatrices.Tridiagonal{<:Any,<:InfFill,<:InfFill,<:InfFill}), - :(LazyBandedMatrices.SymTridiagonal{<:Any,<:InfFill,<:InfFill})) - @eval begin - MemoryLayout(::Type{<:$Typ}) = TridiagonalToeplitzLayout() - BroadcastStyle(::Type{<:$Typ}) = LazyArrayStyle{2}() - end -end - -for Typ in (:(LinearAlgebra.Bidiagonal{<:Any,<:InfFill}), - :(LazyBandedMatrices.Bidiagonal{<:Any,<:InfFill,<:InfFill})) - @eval begin - MemoryLayout(::Type{<:$Typ}) = BidiagonalToeplitzLayout() - BroadcastStyle(::Type{<:$Typ}) = LazyArrayStyle{2}() - end -end - -*(A::LinearAlgebra.Bidiagonal{<:Any,<:InfFill}, B::LinearAlgebra.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() - - -#### -# Banded, TODO: move to extension in InfiniteArrays.jl -### - -sublayout(::ApplyBandedLayout, ::Type{<:Tuple{KR,Integer}}) where {KR<:InfAxes} = - sublayout(PaddedColumns{UnknownLayout}(), Tuple{KR}) \ No newline at end of file diff --git a/src/blockbanded/blockbanded.jl b/src/blockbanded/blockbanded.jl deleted file mode 100644 index 619a154..0000000 --- a/src/blockbanded/blockbanded.jl +++ /dev/null @@ -1,92 +0,0 @@ -const OneToInfCumsum = RangeCumsum{Int,OneToInf{Int}} - -BlockArrays.sortedunion(::AbstractVector{<:PosInfinity}, ::AbstractVector{<:PosInfinity}) = [∞] -function BlockArrays.sortedunion(::AbstractVector{<:PosInfinity}, b) - @assert isinf(length(b)) - b -end - -function BlockArrays.sortedunion(b, ::AbstractVector{<:PosInfinity}) - @assert isinf(length(b)) - b -end -BlockArrays.sortedunion(a::OneToInfCumsum, ::OneToInfCumsum) = a - -BlockBandedMatrices.blockbanded_colstop(A, x::InfiniteCardinal{0}) = x -BlockBandedMatrices.blockbanded_rowstop(A, x::InfiniteCardinal{0}) = x - -BlockArrays.blocklasts(a::InfRanges) = Fill(length(a),1) - -BlockArrays.findblock(::BlockedOneTo, ::RealInfinity) = Block(ℵ₀) - -function BlockArrays.sortedunion(a::Vcat{Int,1,<:Tuple{Union{Int,AbstractVector{Int}},<:AbstractRange}}, - b::Vcat{Int,1,<:Tuple{Union{Int,AbstractVector{Int}},<:AbstractRange}}) - @assert a == b # TODO: generailse? Not sure how to do so in a type stable fashion - a -end - -sizes_from_blocks(A::AbstractVector, ::Tuple{OneToInf{Int}}) = (map(length,A),) -length(::BlockedOneTo{Int,<:InfRanges}) = ℵ₀ - -const OneToInfBlocks = BlockedOneTo{Int,OneToInfCumsum} -const OneToBlocks = BlockedOneTo{Int,OneToCumsum} - -axes(a::OneToInfBlocks) = (a,) -axes(a::OneToBlocks) = (a,) - -LazyBandedMatrices.unitblocks(a::OneToInf) = blockedrange(Ones{Int}(length(a))) - -BlockArrays.dimlength(start, ::Infinity) = ℵ₀ - -function copy(bc::Broadcasted{<:BroadcastStyle,<:Any,typeof(*),<:Tuple{Ones{T,1,Tuple{OneToInfBlocks}},AbstractArray{V,N}}}) where {N,T,V} - a,b = bc.args - @assert bc.axes == axes(b) - convert(AbstractArray{promote_type(T,V),N}, b) -end - -function copy(bc::Broadcasted{<:BroadcastStyle,<:Any,typeof(*),<:Tuple{AbstractArray{T,N},Ones{V,1,Tuple{OneToInfBlocks}}}}) where {N,T,V} - a,b = bc.args - @assert bc.axes == axes(a) - convert(AbstractArray{promote_type(T,V),N}, a) -end - -_block_interlace_axes(::Int, ax::Tuple{BlockedOneTo{Int,OneToInf{Int}}}...) = (blockedrange(Fill(length(ax), ∞)),) - -_block_interlace_axes(nbc::Int, ax::NTuple{2,BlockedOneTo{Int,OneToInf{Int}}}...) = - (blockedrange(Fill(length(ax) ÷ nbc, ∞)),blockedrange(Fill(mod1(length(ax),nbc), ∞))) - - -include("infblocktridiagonal.jl") - - -####### -# block broadcasted -###### - - -BroadcastStyle(::Type{<:SubArray{T,N,Arr,<:NTuple{N,BlockSlice{BlockRange{1,Tuple{II}}}},false}}) where {T,N,Arr<:BlockArray,II<:InfRanges} = - LazyArrayStyle{N}() - -# TODO: generalise following -BroadcastStyle(::Type{<:BlockArray{T,N,<:AbstractArray{<:AbstractArray{T,N},N},<:NTuple{N,BlockedOneTo{Int,<:InfRanges}}}}) where {T,N} = LazyArrayStyle{N}() -# BroadcastStyle(::Type{<:BlockedArray{T,N,<:AbstractArray{T,N},<:NTuple{N,BlockedOneTo{Int,<:InfRanges}}}}) where {T,N} = LazyArrayStyle{N}() -BroadcastStyle(::Type{<:BlockArray{T,N,<:AbstractArray{<:AbstractArray{T,N},N},<:NTuple{N,BlockedOneTo{Int,<:RangeCumsum{Int,<:InfRanges}}}}}) where {T,N} = LazyArrayStyle{N}() -# BroadcastStyle(::Type{<:BlockedArray{T,N,<:AbstractArray{T,N},<:NTuple{N,BlockedOneTo{Int,<:RangeCumsum{Int,<:InfRanges}}}}}) where {T,N} = LazyArrayStyle{N}() - - -### -# KronTrav -### - -_krontrav_axes(A::OneToInf{Int}, B::OneToInf{Int}) = blockedrange(oneto(length(A))) - - -struct InfKronTravBandedBlockBandedLayout <: AbstractLazyBandedBlockBandedLayout end -MemoryLayout(::Type{<:KronTrav{<:Any,2,<:Any,NTuple{2,BlockedOneTo{Int,OneToInfCumsum}}}}) = InfKronTravBandedBlockBandedLayout() - -sublayout(::InfKronTravBandedBlockBandedLayout, ::Type{<:NTuple{2,BlockSlice1}}) = BroadcastBandedLayout{typeof(*)}() -sublayout(::InfKronTravBandedBlockBandedLayout, ::Type{<:NTuple{2,BlockSlice{BlockRange{1,Tuple{OneTo{Int}}}}}}) = KronTravBandedBlockBandedLayout() - -copy(M::Mul{InfKronTravBandedBlockBandedLayout, InfKronTravBandedBlockBandedLayout}) = KronTrav((krontravargs(M.A) .* krontravargs(M.B))...) - -_broadcast_sub_arguments(::InfKronTravBandedBlockBandedLayout, M, V) = _broadcast_sub_arguments(KronTravBandedBlockBandedLayout(), M, V) diff --git a/src/blockbanded/infblocktridiagonal.jl b/src/blockbanded/infblocktridiagonal.jl deleted file mode 100644 index 9786921..0000000 --- a/src/blockbanded/infblocktridiagonal.jl +++ /dev/null @@ -1,89 +0,0 @@ -const BlockTriPertToeplitz{T} = BlockMatrix{T,Tridiagonal{Matrix{T},Vcat{Matrix{T},1,Tuple{Vector{Matrix{T}},Fill{Matrix{T},1,Tuple{OneToInf{Int}}}}}}, - NTuple{2,BlockedOneTo{Int,Vcat{Int,1,Tuple{Vector{Int},InfStepRange{Int,Int}}}}}} - -const BlockTridiagonalToeplitzLayout = BlockLayout{TridiagonalToeplitzLayout,DenseColumnMajor} - -function BlockTridiagonal(adjA::Adjoint{T,BlockTriPertToeplitz{T}}) where T - A = parent(adjA) - BlockTridiagonal(Matrix.(adjoint.(A.blocks.du)), - Matrix.(adjoint.(A.blocks.d)), - Matrix.(adjoint.(A.blocks.dl))) -end - -for op in (:-, :+) - @eval begin - function $op(A::BlockTriPertToeplitz{T}, λ::UniformScaling) where T - TV = promote_type(T,eltype(λ)) - BlockTridiagonal(Vcat(convert.(AbstractVector{Matrix{TV}}, A.blocks.dl.args)...), - Vcat(convert.(AbstractVector{Matrix{TV}}, broadcast($op, A.blocks.d, Ref(λ)).args)...), - Vcat(convert.(AbstractVector{Matrix{TV}}, A.blocks.du.args)...)) - end - function $op(λ::UniformScaling, A::BlockTriPertToeplitz{V}) where V - TV = promote_type(eltype(λ),V) - BlockTridiagonal(Vcat(convert.(AbstractVector{Matrix{TV}}, broadcast($op, A.blocks.dl.args))...), - Vcat(convert.(AbstractVector{Matrix{TV}}, broadcast($op, Ref(λ), A.blocks.d).args)...), - Vcat(convert.(AbstractVector{Matrix{TV}}, broadcast($op, A.blocks.du.args))...)) - end - $op(adjA::Adjoint{T,BlockTriPertToeplitz{T}}, λ::UniformScaling) where T = $op(BlockTridiagonal(adjA), λ) - $op(λ::UniformScaling, adjA::Adjoint{T,BlockTriPertToeplitz{T}}) where T = $op(λ, BlockTridiagonal(adjA)) - end -end - -*(a::AbstractVector, b::AbstractFill{<:Any,2,Tuple{OneTo{Int},OneToInf{Int}}}) = ApplyArray(*,a,b) - -sizes_from_blocks(A::Diagonal, ::NTuple{2,OneToInf{Int}}) = size.(A.diag, 1), size.(A.diag,2) -sizes_from_blocks(A::Tridiagonal, ::NTuple{2,OneToInf{Int}}) = size.(A.d, 1), size.(A.d,2) -sizes_from_blocks(A::LazyBandedMatrices.Tridiagonal, ::NTuple{2,OneToInf{Int}}) = size.(A.d, 1), size.(A.d,2) -sizes_from_blocks(A::Bidiagonal, ::NTuple{2,OneToInf{Int}}) = size.(A.dv, 1), size.(A.dv,2) -sizes_from_blocks(A::LazyBandedMatrices.Bidiagonal, ::NTuple{2,OneToInf{Int}}) = size.(A.dv, 1), size.(A.dv,2) - -axes_print_matrix_row(_, io, X, A, i, ::AbstractVector{<:PosInfinity}, sep) = nothing -axes_print_matrix_row(::NTuple{2,AbstractBlockedUnitRange}, io, X, A, i, ::AbstractVector{<:PosInfinity}, sep) = nothing - - -function BlockSkylineSizes(A::BlockTriPertToeplitz, (l,u)::NTuple{2,Int}) - N = max(length(A.blocks.du.args[1])+1,length(A.blocks.d.args[1]),length(A.blocks.dl.args[1])) - block_sizes = Vector{Int}(undef, N) # assume square - block_starts = BandedMatrix{Int}(undef, (N+l,N), (l,u)) - block_strides = Vector{Int}(undef, N) - for J=1:N - block_starts[max(1,J-u),J] = J == 1 ? 1 : - block_starts[max(1,J-1-u),J-1]+block_sizes[J-1]*block_strides[J-1] - - for K=max(1,J-u)+1:J+l - block_starts[K,J] = block_starts[K-1,J]+size(A[Block(K-1,J)],1) - end - block_strides[J] = block_starts[J+l,J] + size(A[Block(J+l,J)],1) - block_starts[max(1,J-u),J] - block_sizes[J] = size(A[Block(J,J)],2) - end - - block_stride∞ = 0 - for K=max(1,N+1-u):N+1+l - block_stride∞ += size(A[Block(K,N+1)],1) - end - block_size∞ = size(A[Block(N+1,N+1)],2) - - bs∞ = fill(block_starts[max(1,N-u),N]+block_strides[N]*size(A[Block(N,N)],2):block_stride∞*block_size∞:∞, l+u+1) - for k=2:l+u+1 - bs∞[k] = bs∞[k-1] .+ size(A[Block(N+1-u+k-1,N+1)],1) - end - - BlockSkylineSizes(axes(A), - _BandedMatrix(Hcat(block_starts.data, Vcat(adjoint.(bs∞)...)), ℵ₀, l, u), - Vcat(block_strides, Fill(block_stride∞,∞)), - Fill(l,∞),Fill(u,∞)) -end - -function BlockBandedMatrix(A::BlockTriPertToeplitz{T}, (l,u)::NTuple{2,Int}) where T - data = T[] - append!(data,vec(A[Block.(1:1+l),Block(1)])) - N = max(length(A.blocks.du.args[1])+1,length(A.blocks.d.args[1]),length(A.blocks.dl.args[1])) - for J=2:N - append!(data, vec(A[Block.(max(1,J-u):J+l),Block(J)])) - end - tl = mortar(Fill(vec(A[Block.(max(1,N+1-u):N+1+l),Block(N+1)]),∞)) - - B = _BlockSkylineMatrix(Vcat(data,tl), BlockSkylineSizes(A, (l,u))) -end - -BlockBandedMatrix(A::BlockTriPertToeplitz) = BlockBandedMatrix(A, blockbandwidths(A)) diff --git a/test/runtests.jl b/test/runtests.jl index c76adc4..df60a8e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,9 +1,9 @@ using InfiniteLinearAlgebra, BlockBandedMatrices, BlockArrays, BandedMatrices, InfiniteArrays, FillArrays, LazyArrays, Test, MatrixFactorizations, ArrayLayouts, LinearAlgebra, Random, LazyBandedMatrices, StaticArrays import InfiniteLinearAlgebra: qltail, toeptail, tailiterate, tailiterate!, tail_de, ql_X!, - InfToeplitz, PertToeplitz, TriToeplitz, InfBandedMatrix, InfBandCartesianIndices, - rightasymptotics, QLHessenberg, ConstRows, PertConstRows, chop, chop!, pad, - BandedToeplitzLayout, PertToeplitzLayout, TridiagonalToeplitzLayout, BidiagonalToeplitzLayout, + InfToeplitz, PertToeplitz, TriToeplitz, InfBandedMatrix, + rightasymptotics, QLHessenberg, PertConstRows, chop, chop!, pad, + PertToeplitzLayout, TridiagonalToeplitzLayout, BidiagonalConjugation import Base: BroadcastStyle, oneto import BlockArrays: _BlockArray, blockcolsupport, findblock @@ -52,429 +52,7 @@ end @test pad(BlockVec(transpose(X)), blockedrange(Fill(3,∞))) isa BlockVec{Int,<:Transpose} end -include("test_infbanded.jl") -@testset "∞-block arrays" begin - @testset "fixed block size" begin - k = Base.OneTo.(oneto(∞)) - n = Fill.(oneto(∞), oneto(∞)) - @test broadcast(length, k) ≡ map(length, k) ≡ OneToInf() - @test broadcast(length, n) ≡ map(length, n) ≡ OneToInf() - - b = mortar(Fill([1, 2], ∞)) - @test blockaxes(b, 1) ≡ Block.(OneToInf()) - @test b[Block(5)] == [1, 2] - @test b[Block.(2:∞)][Block.(2:10)] == b[Block.(3:11)] - @test exp.(b)[Block.(2:∞)][Block.(2:10)] == exp.(b[Block.(3:11)]) - - @test blockedrange(Vcat(2, Fill(3, ∞))) isa BlockedOneTo{<:Any,<:InfiniteArrays.InfStepRange} - - c = BlockedArray(1:∞, Vcat(2, Fill(3, ∞))) - @test c[Block.(2:∞)][Block.(2:10)] == c[Block.(3:11)] - - @test length(axes(b, 1)) ≡ ℵ₀ - @test last(axes(b, 1)) ≡ ℵ₀ - @test Base.BroadcastStyle(typeof(b)) isa LazyArrayStyle{1} - - @test unitblocks(oneto(∞)) ≡ blockedrange(Ones{Int}(∞)) - @test unitblocks(2:∞) == 2:∞ - - @test unitblocks(oneto(∞))[Block.(2:∞)] == 2:∞ - end - - @testset "1:∞ blocks" begin - a = blockedrange(oneto(∞)) - @test axes(a, 1) == a - o = Ones((a,)) - @test Base.BroadcastStyle(typeof(a)) isa LazyArrayStyle{1} - b = exp.(a) - @test axes(b, 1) == a - @test o .* b isa typeof(b) - @test b .* o isa typeof(b) - end - - @testset "padded" begin - c = BlockedArray([1; zeros(∞)], Vcat(2, Fill(3, ∞))) - @test c + c isa BlockedVector - end - - @testset "concat" begin - a = unitblocks(1:∞) - b = exp.(a) - c = BlockBroadcastArray(vcat, a, b) - @test length(c) == ∞ - @test blocksize(c) == (∞,) - @test c[Block(5)] == [a[5], b[5]] - - A = unitblocks(BandedMatrix(0 => 1:∞, 1 => Fill(2.0, ∞), -1 => Fill(3.0, ∞))) - B = BlockBroadcastArray(hvcat, 2, A, Zeros(axes(A)), Zeros(axes(A)), A) - @test B[Block(3, 3)] == [A[3, 3] 0; 0 A[3, 3]] - @test B[Block(3, 4)] == [A[3, 4] 0; 0 A[3, 4]] - @test B[Block(3, 5)] == [A[3, 5] 0; 0 A[3, 5]] - @test blockbandwidths(B) == (1, 1) - @test subblockbandwidths(B) == (0, 0) - @test B[Block.(1:10), Block.(1:10)] isa BlockSkylineMatrix - - C = BlockBroadcastArray(hvcat, 2, A, A, A, A) - @test C[Block(3, 3)] == fill(A[3, 3], 2, 2) - @test C[Block(3, 4)] == fill(A[3, 4], 2, 2) - @test C[Block(3, 5)] == fill(A[3, 5], 2, 2) - @test blockbandwidths(C) == (1, 1) - @test subblockbandwidths(C) == (1, 1) - @test B[Block.(1:10), Block.(1:10)] isa BlockSkylineMatrix - end - - @testset "DiagTrav" begin - C = zeros(∞,∞); - C[1:3,1:4] .= [1 2 3 4; 5 6 7 8; 9 10 11 12] - A = DiagTrav(C) - @test blockcolsupport(A) == Block.(1:6) - @test A[Block.(1:7)] == [1; 5; 2; 9; 6; 3; 0; 10; 7; 4; 0; 0; 11; 8; 0; 0; 0; 0; 12; zeros(9)] - - C = zeros(∞,4); - C[1:3,1:4] .= [1 2 3 4; 5 6 7 8; 9 10 11 12] - A = DiagTrav(C) - @test A[Block.(1:7)] == [1; 5; 2; 9; 6; 3; 0; 10; 7; 4; 0; 0; 11; 8; 0; 0; 0; 12; zeros(4)] - - C = zeros(3,∞); - C[1:3,1:4] .= [1 2 3 4; 5 6 7 8; 9 10 11 12] - A = DiagTrav(C) - @test A[Block.(1:7)] == [1; 5; 2; 9; 6; 3; 10; 7; 4; 11; 8; 0; 12; zeros(5)] - end - - @testset "KronTrav" begin - Δ = BandedMatrix(1 => Ones(∞), -1 => Ones(∞)) / 2 - A = KronTrav(Δ - 2I, Eye(∞)) - @test axes(A, 1) isa InfiniteLinearAlgebra.OneToInfBlocks - V = view(A, Block.(Base.OneTo(3)), Block.(Base.OneTo(3))) - - @test MemoryLayout(A) isa InfiniteLinearAlgebra.InfKronTravBandedBlockBandedLayout - @test MemoryLayout(V) isa LazyBandedMatrices.KronTravBandedBlockBandedLayout - - @test A[Block.(Base.OneTo(3)), Block.(Base.OneTo(3))] isa KronTrav - - u = A * [1; zeros(∞)] - @test u[1:3] == A[1:3, 1] - @test bandwidths(view(A, Block(1, 1))) == (1, 1) - - @test A*A isa KronTrav - @test (A*A)[Block.(Base.OneTo(3)), Block.(Base.OneTo(3))] ≈ A[Block.(1:3), Block.(1:4)]A[Block.(1:4), Block.(1:3)] - end - - @testset "triangle recurrences" begin - @testset "n and k" begin - n = mortar(Fill.(oneto(∞), oneto(∞))) - k = mortar(Base.OneTo.(oneto(∞))) - - @test n[Block(5)] ≡ layout_getindex(n, Block(5)) ≡ view(n, Block(5)) ≡ Fill(5, 5) - @test k[Block(5)] ≡ layout_getindex(k, Block(5)) ≡ view(k, Block(5)) ≡ Base.OneTo(5) - @test Base.BroadcastStyle(typeof(n)) isa LazyArrays.LazyArrayStyle{1} - @test Base.BroadcastStyle(typeof(k)) isa LazyArrays.LazyArrayStyle{1} - - N = 1000 - v = view(n, Block.(Base.OneTo(N))) - @test view(v, Block(2)) ≡ Fill(2, 2) - @test axes(v) isa Tuple{BlockedOneTo{Int,ArrayLayouts.RangeCumsum{Int64,Base.OneTo{Int64}}}} - @test @allocated(axes(v)) ≤ 40 - - dest = BlockedArray{Float64}(undef, axes(v)) - @test copyto!(dest, v) == v - @test @allocated(copyto!(dest, v)) ≤ 40 - - v = view(k, Block.(Base.OneTo(N))) - @test view(v, Block(2)) ≡ Base.OneTo(2) - @test axes(v) isa Tuple{BlockedOneTo{Int,ArrayLayouts.RangeCumsum{Int64,Base.OneTo{Int64}}}} - @test @allocated(axes(v)) ≤ 40 - @test copyto!(dest, v) == v - - @testset "stack overflow" begin - i = Base.to_indices(k, (Block.(2:∞),))[1].indices - last(i) - end - - v = view(k, Block.(2:∞)) - @test Base.BroadcastStyle(typeof(v)) isa LazyArrayStyle{1} - @test v[Block(1)] == 1:2 - @test v[Block(1)] ≡ k[Block(2)] ≡ Base.OneTo(2) - - @test axes(n, 1) isa BlockedOneTo{Int,ArrayLayouts.RangeCumsum{Int64,OneToInf{Int64}}} - end - - @testset "BlockHcat copyto!" begin - n = mortar(Fill.(oneto(∞), oneto(∞))) - k = mortar(Base.OneTo.(oneto(∞))) - - a = b = c = 0.0 - dat = BlockHcat( - BroadcastArray((n, k, b, bc1) -> (k + b - 1) * (n + k + bc1) / (2k + bc1), n, k, b, b + c - 1), - BroadcastArray((n, k, abc, bc, bc1) -> (n + k + abc) * (k + bc) / (2k + bc1), n, k, a + b + c, b + c, b + c - 1) - ) - N = 1000 - KR = Block.(Base.OneTo(N)) - V = view(dat, Block.(Base.OneTo(N)), :) - @test MemoryLayout(V) isa LazyArrays.ApplyLayout{typeof(hcat)} - @test BlockedArray(V)[Block.(1:5), :] == dat[Block.(1:5), :] - V = view(dat', :, Block.(Base.OneTo(N))) - @test MemoryLayout(V) isa LazyArrays.ApplyLayout{typeof(vcat)} - a = dat.arrays[1]' - N = 100 - KR = Block.(Base.OneTo(N)) - v = view(a, :, KR) - @time r = BlockedArray(v) - @test v == r - end - - @testset "BlockBanded" begin - a = b = c = 0.0 - n = mortar(Fill.(oneto(∞), oneto(∞))) - k = mortar(Base.OneTo.(oneto(∞))) - Dy = BlockBandedMatrices._BandedBlockBandedMatrix((k .+ (b + c))', axes(k, 1), (-1, 1), (-1, 1)) - N = 100 - @test Dy[Block.(1:N), Block.(1:N)] == BlockBandedMatrices._BandedBlockBandedMatrix((k.+(b+c))[Block.(1:N)]', axes(k, 1)[Block.(1:N)], (-1, 1), (-1, 1)) - @test colsupport(Dy, axes(Dy,2)) == 1:∞ - @test rowsupport(Dy, axes(Dy,1)) == 2:∞ - end - - @testset "Symmetric" begin - k = mortar(Base.OneTo.(oneto(∞))) - n = mortar(Fill.(oneto(∞), oneto(∞))) - - dat = BlockHcat( - BlockBroadcastArray(hcat, float.(k), Zeros((axes(n, 1),)), float.(n)), - Zeros((axes(n, 1), Base.OneTo(3))), - Zeros((axes(n, 1), Base.OneTo(3)))) - M = BlockBandedMatrices._BandedBlockBandedMatrix(dat', axes(k, 1), (1, 1), (1, 1)) - Ms = Symmetric(M) - @test blockbandwidths(M) == (1, 1) - @test blockbandwidths(Ms) == (1, 1) - @test Ms[Block.(1:5), Block.(1:5)] == Symmetric(M[Block.(1:5), Block.(1:5)]) - @test Ms[Block.(1:5), Block.(1:5)] isa BandedBlockBandedMatrix - - b = [ones(10); zeros(∞)] - @test (Ms * b)[Block.(1:6)] == Ms[Block.(1:6), Block.(1:4)]*ones(10) - @test ((Ms * Ms) * b)[Block.(1:6)] == (Ms * (Ms * b))[Block.(1:6)] - @test ((Ms + Ms) * b)[Block.(1:6)] == (2*(Ms * b))[Block.(1:6)] - - dat = BlockBroadcastArray(hcat, float.(k), Zeros((axes(n, 1),)), float.(n)) - M = BlockBandedMatrices._BandedBlockBandedMatrix(dat', axes(k, 1), (-1, 1), (1, 1)) - Ms = Symmetric(M) - @test Symmetric((M+M)[Block.(1:10), Block.(1:10)]) == (Ms+Ms)[Block.(1:10), Block.(1:10)] - end - end - - @testset "blockdiag" begin - D = Diagonal(mortar(Fill.((-(0:∞) - (0:∞) .^ 2), 1:2:∞))) - x = [randn(5); zeros(∞)] - x̃ = BlockedArray(x, (axes(D, 1),)) - @test (D*x)[1:10] == (D*x̃)[1:10] - end - - @testset "sortedunion" begin - a = cumsum(1:2:∞) - @test BlockArrays.sortedunion(a, a) ≡ a - @test BlockArrays.sortedunion([∞], a) ≡ BlockArrays.sortedunion(a, [∞]) ≡ a - @test BlockArrays.sortedunion([∞], [∞]) == [∞] - - b = Vcat([1, 2], 3:∞) - c = Vcat(1, 3:∞) - @test BlockArrays.sortedunion(b, b) ≡ b - @test BlockArrays.sortedunion(c, c) ≡ c - end -end - - - -@testset "Algebra" begin - @testset "BandedMatrix" 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(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 "BlockTridiagonal" begin - A = BlockTridiagonal(Vcat([fill(1.0, 2, 1), Matrix(1.0I, 2, 2), Matrix(1.0I, 2, 2), Matrix(1.0I, 2, 2)], Fill(Matrix(1.0I, 2, 2), ∞)), - Vcat([zeros(1, 1)], Fill(zeros(2, 2), ∞)), - Vcat([fill(1.0, 1, 2), Matrix(1.0I, 2, 2)], Fill(Matrix(1.0I, 2, 2), ∞))) - - @test A isa InfiniteLinearAlgebra.BlockTriPertToeplitz - @test isblockbanded(A) - - @test A[Block.(1:2), Block(1)] == A[1:3, 1:1] == reshape([0.0, 1.0, 1.0], 3, 1) - - @test BlockBandedMatrix(A)[1:100, 1:100] == BlockBandedMatrix(A, (2, 1))[1:100, 1:100] == BlockBandedMatrix(A, (1, 1))[1:100, 1:100] == A[1:100, 1:100] - - @test (A-I)[1:100, 1:100] == A[1:100, 1:100] - I - @test (A+I)[1:100, 1:100] == A[1:100, 1:100] + I - @test (I+A)[1:100, 1:100] == I + A[1:100, 1:100] - @test (I-A)[1:100, 1:100] == I - A[1:100, 1:100] - - @test (A-im*I)[1:100, 1:100] == A[1:100, 1:100] - im * I - @test (A+im*I)[1:100, 1:100] == A[1:100, 1:100] + im * I - @test (im*I+A)[1:100, 1:100] == im * I + A[1:100, 1:100] - @test (im*I-A)[1:100, 1:100] == im * I - A[1:100, 1:100] - - T = mortar(LazyBandedMatrices.Tridiagonal(Fill([1 2; 3 4], ∞), Fill([1 2; 3 4], ∞), Fill([1 2; 3 4], ∞))); - #TODO: copy BlockBidiagonal code from BlockBandedMatrices to LazyBandedMatrices - @test T[Block(2, 2)] == [1 2; 3 4] - @test_broken T[Block(1, 3)] == Zeros(2, 2) - end - - @testset "BlockBidiagonal" begin - B = mortar(LazyBandedMatrices.Bidiagonal(Fill([1 2; 3 4], ∞), Fill([1 2; 3 4], ∞), :U)); - #TODO: copy BlockBidiagonal code from BlockBandedMatrices to LazyBandedMatrices - @test B[Block(2, 3)] == [1 2; 3 4] - @test_broken B[Block(1, 3)] == Zeros(2, 2) - 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 - - @testset "Triangle OP recurrences" begin - k = mortar(Base.OneTo.(1:∞)) - n = mortar(Fill.(1:∞, 1:∞)) - @test k[Block.(2:3)] isa BlockArray - @test n[Block.(2:3)] isa BlockArray - @test k[Block.(2:3)] == [1, 2, 1, 2, 3] - @test n[Block.(2:3)] == [2, 2, 3, 3, 3] - @test blocksize(BroadcastVector(exp, k)) == (ℵ₀,) - @test BroadcastVector(exp, k)[Block.(2:3)] == exp.([1, 2, 1, 2, 3]) - # BroadcastVector(+,k,n) - end - # Multivariate OPs Corollary (3) - # n = 5 - # BlockTridiagonal(Zeros.(1:∞,2:∞), - # (n -> Diagonal(((n+2).+(0:n)))/ (2n + 2)).(0:∞), - # Zeros.(2:∞,1:∞)) - - @testset "KronTrav" begin - Δ = BandedMatrix(1 => Ones(∞) / 2, -1 => Ones(∞)) - A = KronTrav(Δ, Eye(∞)) - @test A[Block(100, 101)] isa BandedMatrix - @test A[Block(100, 100)] isa BandedMatrix - @test A[Block.(1:5), Block.(1:5)] isa BandedBlockBandedMatrix - B = KronTrav(Eye(∞), Δ) - @test B[Block(100, 101)] isa BandedMatrix - @test B[Block(100, 100)] isa BandedMatrix - V = view(A + B, Block.(1:5), Block.(1:5)) - @test MemoryLayout(typeof(V)) isa BroadcastBandedBlockBandedLayout{typeof(+)} - @test arguments(V) == (A[Block.(1:5), Block.(1:5)], B[Block.(1:5), Block.(1:5)]) - @test (A+B)[Block.(1:5), Block.(1:5)] == A[Block.(1:5), Block.(1:5)] + B[Block.(1:5), Block.(1:5)] - - @test blockbandwidths(A + B) == (1, 1) - @test blockbandwidths(2A) == (1, 1) - @test blockbandwidths(2 * (A + B)) == (1, 1) - - @test subblockbandwidths(A + B) == (1, 1) - @test subblockbandwidths(2A) == (1, 1) - @test subblockbandwidths(2 * (A + B)) == (1, 1) - end -end - -@testset "findblock at +∞, HarmonicOrthogonalPolynomials#88" begin - @test findblock(blockedrange(1:2:∞), RealInfinity()) == Block(ℵ₀) -end include("test_hessenbergq.jl") include("test_infql.jl") diff --git a/test/test_infbanded.jl b/test/test_infbanded.jl deleted file mode 100644 index 5a383e0..0000000 --- a/test/test_infbanded.jl +++ /dev/null @@ -1,191 +0,0 @@ -using InfiniteLinearAlgebra, ArrayLayouts, InfiniteArrays, BandedMatrices, FillArrays, LazyBandedMatrices, LazyArrays, Test -import BandedMatrices: _BandedMatrix, bandeddata -using InfiniteLinearAlgebra: InfToeplitz, ConstRows, BandedToeplitzLayout, TridiagonalToeplitzLayout, BidiagonalToeplitzLayout, PertToeplitz, PertToeplitzLayout -using Base: oneto -using LazyArrays: simplifiable - -@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 LazyArrays.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 - @test MemoryLayout(LazyBandedMatrices.Tridiagonal(Fill(1,∞), Zeros(∞), Fill(3,∞))) isa TridiagonalToeplitzLayout - @test MemoryLayout(LazyBandedMatrices.Bidiagonal(Fill(1,∞), Zeros(∞), :U)) isa BidiagonalToeplitzLayout - @test MemoryLayout(LazyBandedMatrices.SymTridiagonal(Fill(1,∞), Zeros(∞))) isa TridiagonalToeplitzLayout - - T = LazyBandedMatrices.Tridiagonal(Fill(1,∞), Zeros(∞), Fill(3,∞)) - @test T[2:∞,3:∞] isa SubArray - @test exp.(T) isa BroadcastMatrix - @test exp.(T)[2:∞,3:∞][1:10,1:10] == exp.(T[2:∞,3:∞])[1:10,1:10] == exp.(T[2:11,3:12]) - @test exp.(T)[2:∞,3:∞] isa BroadcastMatrix - @test exp.(T[2:∞,3:∞]) isa BroadcastMatrix - - B = LazyBandedMatrices.Bidiagonal(Fill(1,∞), Zeros(∞), :U) - @test B[2:∞,3:∞] isa SubArray - @test exp.(B) isa BroadcastMatrix - @test exp.(B)[2:∞,3:∞][1:10,1:10] == exp.(B[2:∞,3:∞])[1:10,1:10] == exp.(B[2:11,3:12]) - @test exp.(B)[2:∞,3:∞] isa BroadcastMatrix - - @testset "algebra" begin - T = Tridiagonal(Fill(1,∞), Fill(2,∞), Fill(3,∞)) - @test T isa InfiniteLinearAlgebra.TriToeplitz - @test (T + 2I)[1:10,1:10] == (2I + T)[1:10,1:10] == T[1:10,1:10] + 2I - 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 - 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 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 InfiniteLinearAlgebra.SymTriPertToeplitz - @test (A + 2I)[1:10,1:10] == (2I + A)[1:10,1:10] == A[1:10,1:10] + 2I - - A = Tridiagonal(Vcat([3.,4.], Fill.(0.5,∞)), Vcat([1,2.], Fill(2.,∞)), Vcat([3.,4.], Fill.(0.5,∞))) - @test A isa InfiniteLinearAlgebra.TriPertToeplitz - @test Adjoint(A) isa InfiniteLinearAlgebra.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 - 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 - 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 "Diagonal{Fill} * Bidiagonal" begin - A, B = Diagonal(Fill(2,∞)) , LazyBandedMatrices.Bidiagonal(exp.(1:∞), exp.(1:∞), :L) - @test (A*B)[1:10,1:10] ≈ (B*A)[1:10,1:10] ≈ 2B[1:10,1:10] - 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 -end