Cleanup code

Author: dlfivefifty

Project.toml

@@ -8,19 +8,24 @@ BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
 BlockBandedMatrices = "ffab5731-97b5-5995-9138-79e8c1846df0"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 InfiniteArrays = "4858937d-0d70-526a-a4dd-2d5cb5dd786c"
-IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
 LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MatrixFactorizations = "a3b82374-2e81-5b9e-98ce-41277c0e4c87"
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
 
 [compat]
-BandedMatrices = "0.9"
+BandedMatrices = "0.10"
 BlockArrays = "0.9"
 BlockBandedMatrices = "0.4"
-FillArrays = "≥ 0.5.0"
+FillArrays = "0.6.4"
 InfiniteArrays = "0.1"
-IntervalArithmetic = "≥ 0.15.2"
-LazyArrays = "0.8, 0.9"
-MatrixFactorizations = "0.0.4, 0.1"
-julia = "0.7, 1"
+LazyArrays = "0.10"
+MatrixFactorizations = "0.1"
+julia = "1"
+
+[extras]
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[targets]
+test = ["Test", "Random"]

src/InfiniteLinearAlgebra.jl

@@ -9,7 +9,7 @@ import FillArrays: AbstractFill
 import BandedMatrices: BandedMatrix, _BandedMatrix, bandeddata
 import LinearAlgebra: lmul!,  ldiv!, matprod, qr, QRPackedQ, AbstractTriangular, AbstractQ, adjoint, transpose,
                         QR
-import LazyArrays: CachedArray, VecMulMat, DenseColumnMajor, FillLayout, ApplyMatrix
+import LazyArrays: CachedArray, DenseColumnMajor, FillLayout, ApplyMatrix
 import MatrixFactorizations: ql, ql!, QLPackedQ, getL, reflector!, reflectorApply!,
                             QL
 
@@ -20,8 +20,6 @@ import BandedMatrices: BandedMatrix, bandwidths
 import BlockBandedMatrices: _BlockSkylineMatrix, _BandedMatrix, AbstractBlockSizes, cumulsizes, _BlockSkylineMatrix, BlockSizes, blockstart, blockstride,
         BlockSkylineSizes, BlockSkylineMatrix, BlockBandedMatrix, _BlockBandedMatrix, BlockTridiagonal
 
-import IntervalArithmetic
-import IntervalArithmetic: Interval, emptyinterval
 
 if VERSION < v"1.2-"
     import Base: has_offset_axes
@@ -38,6 +36,5 @@ include("banded/infbanded.jl")
 include("blockbanded/infblocktridiagonal.jl")
 include("banded/infqltoeplitz.jl")
 include("infql.jl")
-include("banded/rigorous.jl")
 
 end # module

src/banded/hessenbergq.jl

@@ -87,6 +87,10 @@ bandwidths(Q::LowerHessenbergQ) = (size(Q,1)-1,1)
 adjoint(Q::UpperHessenbergQ) = LowerHessenbergQ(adjoint.(Q.q))
 adjoint(Q::LowerHessenbergQ) = UpperHessenbergQ(adjoint.(Q.q))
 
+check_mul_axes(A::AbstractHessenbergQ, B, C...) =
+    axes(A,2) == axes(B,1) || throw(DimensionMismatch("Second axis of A, $(axes(A,2)), and first axis of B, $(axes(B,1)) must match"))
+
+
 function lmul!(Q::LowerHessenbergQ{T}, x::AbstractVector) where T
     t = Array{T}(undef, 2)
     for n = 1:length(Q.q)

src/banded/infbanded.jl

@@ -1,6 +1,7 @@
 const TriToeplitz{T} = Tridiagonal{T,Fill{T,1,Tuple{OneToInf{Int}}}}
-const InfToeplitz{T} = BandedMatrix{T,<:ApplyMatrix{T,<:VecMulMat{DenseColumnMajor,FillLayout,T,T}},OneToInf{Int}}
-const PertToeplitz{T} = BandedMatrix{T,<:Hcat{T,<:Tuple{Matrix{T},ApplyMatrix{T,<:VecMulMat{DenseColumnMajor,FillLayout,T,T}}}},OneToInf{Int}}
+const ConstRows{T} = ApplyMatrix{T,typeof(*),<:Tuple{<:AbstractVector,<:AbstractFill}}
+const InfToeplitz{T} = BandedMatrix{T,<:ConstRows{T},OneToInf{Int}}
+const PertToeplitz{T} = BandedMatrix{T,<:Hcat{T,<:Tuple{Matrix{T},<:ConstRows{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}}}}}}

src/banded/reducetoeplitzql.jl

@@ -28,7 +28,7 @@ for op in (:-, :+)
     end
 end
 
-*(a::AbstractVector, b::AbstractFill{<:Any,2,Tuple{OneTo{Int},OneToInf{Int}}}) = MulArray(a,b)
+*(a::AbstractVector, b::AbstractFill{<:Any,2,Tuple{OneTo{Int},OneToInf{Int}}}) = ApplyArray(*,a,b)
 
 
 sizes_from_blocks(A::AbstractVector, ::Tuple{OneToInf{Int}}) = BlockSizes((Vcat(1, 1 .+ cumsum(length.(A))),))

src/banded/rigorous.jl

@@ -1,6 +1,6 @@
-using InfiniteLinearAlgebra, BlockBandedMatrices, BlockArrays, BandedMatrices, InfiniteArrays, FillArrays, LazyArrays, Test, MatrixFactorizations, LinearAlgebra
+using InfiniteLinearAlgebra, BlockBandedMatrices, BlockArrays, BandedMatrices, InfiniteArrays, FillArrays, LazyArrays, Test, MatrixFactorizations, LinearAlgebra, Random
 import InfiniteLinearAlgebra: qltail, toeptail, tailiterate , tailiterate!, tail_de, ql_X!,
-                    InfToeplitz, PertToeplitz, TriToeplitz, InfBandedMatrix, householderparams, combine_two_Q, periodic_combine_two_Q, householderparams,
+                    InfToeplitz, PertToeplitz, TriToeplitz, InfBandedMatrix, 
                     rightasymptotics, QLHessenberg
 import BlockBandedMatrices: isblockbanded, _BlockBandedMatrix
 import MatrixFactorizations: QLPackedQ
@@ -44,7 +44,6 @@ end
 
 include("test_hessenbergq.jl")
 
-
 @testset "PertTriToeplitz QL" begin
     A = Tridiagonal(Vcat(Float64[], Fill(2.0,∞)), 
                     Vcat(Float64[2.0], Fill(0.0,∞)), 
@@ -55,7 +54,6 @@ include("test_hessenbergq.jl")
     end
 end
 
-
 @testset "Pert Hessenberg Toeplitz" begin
     a = [1,2,5,0.5]
     Random.seed!(0)
@@ -100,7 +98,6 @@ end
     end
 end
 
-
 @testset "Pert faux-periodic QL" begin
     a = [0.5794879759059747 + 0.0im,0.538107104952824 - 0.951620830938543im,-0.19352887774167749 - 0.3738926065520737im,0.4314153362874331,0.0]
     T = _BandedMatrix(a*Ones{ComplexF64}(1,∞), ∞, 3,1)
@@ -119,86 +116,3 @@ end
     @test Qn[1:10,1:10] * diagm(0 => [Ones(5); -(-1).^(1:5)]) ≈ Q[1:10,1:10]
     @test diagm(0 => [Ones(5); -(-1).^(1:5)]) * Ln[1:10,1:10] ≈ L[1:10,1:10]
 end
-
-
-function qdL(A)
-    l,u = bandwidths(A)
-    H = _BandedMatrix(A.data, ∞, l+u-1, 1)
-    Q1,L1 = ql(H)
-    D1, Q1, L1 = reduceband(A)
-    T2 = _BandedMatrix(Lrightasymptotics(L1), ∞, l, u)
-    l1 = L1[1,1]
-    A2 = [[D1 l1 zeros(1,10-size(D1,2)-1)]; T2[1:10-1,1:10]] # TODO: remove
-    B2 = _BandedMatrix(T2.data, ∞, l+u-2, 2)
-    B2 = _BandedMatrix(T2.data, ∞, l+u-2, 2)
-    D2, Q2, L2 = reduceband(B2)
-    l2 = L2[1,1]
-    # peroidic tail
-    T3 = _BandedMatrix(Lrightasymptotics(L2), ∞, l+1, u-1)
-    A3 = [[D2 l2 zeros(1,10-size(D2,2)-1)]; T3[1:10-1,1:10]] # TODO: remove
-
-    Q3,L3 = ql( [A2[1,1] A2[1:1,2:3]; [Q2[1:3,1:1]' * T2[1:3,1]  A3[1:1,1:2] ]])
-
-    fd_data = hcat([0; L3[:,1]; Q2[1:3,2:3]' * T2[1:3,1]], [L3[:,2]; T3[1:3,1]], [L3[2,3]; T3[1:4,2]])
-    B3 = _BandedMatrix(Hcat(fd_data, T3.data), ∞, l+u-1, 1)
-
-    ql(B3).L
-end
-
-@testset "quick-and-dirty L" begin
-    for λ in (5,1,0.1+0.1im,-0.5-0.1im), A in (BandedMatrix(3 => Fill(7/10,∞), 2 => Fill(1,∞), 0 => Fill(-λ,∞), -1 => Fill(2im,∞)),
-                                        BandedMatrix(3 => Fill(7/10,∞), 2 => Fill(1,∞), 0 => Fill(-conj(λ),∞), -1 => Fill(-2im,∞)))
-        L∞ = qdL(A)[1:10,1:10]
-        Ln = ql(A[1:1000,1:1000]).L[1:10,1:10]
-        @test L∞ .* sign.(diag(L∞)) ≈ Matrix(Ln) .* sign.(diag(Ln))
-    end
-    for λ in (-3-0.1im, 0.0, -1im)
-        A = BandedMatrix(3 => Fill(7/10,∞), 2 => Fill(1,∞), 0 => Fill(-conj(λ),∞), -1 => Fill(-2im,∞))
-        @test abs(qdL(A)[1,1]) ≈ abs(ql(A[1:10000,1:10000]).L[1,1])
-    end
-    for λ in (1+2im,)
-        A = BandedMatrix(3 => Fill(7/10,∞), 2 => Fill(1,∞), 0 => Fill(-λ,∞), -1 => Fill(2im,∞))
-        @test_throws DomainError qdL(A)
-    end
-end
-
-
-
-function tail_de_j(a::AbstractVector{T}, j) where T
-    m = length(a)
-    C = [view(a,m-1:-1:1) Vcat(-a[end]*Eye(m-2), Zeros{T}(1,m-2))]
-    λ, V = eigen(C)::Eigen{T,T,Matrix{T},Vector{T}}
-    n2 = abs2.(λ[j])
-    n2 ≥ abs2(a[end]) || throw(DomainError(a, "QL factorization does not exist. This could indicate that the operator is not Fredholm or that the dimension of the kernel exceeds that of the co-kernel. Try again with the adjoint."))
-    c_abs = sqrt((n2 - abs2(a[end]))/abs2(V[1,j]))
-    c_sgn = -sign(λ[j])/sign(V[1,j]*a[end-1] - V[2,j]*a[end])
-    c_sgn*c_abs*V[end:-1:1,j]    
-end
-
-@testset "non-uniqueness" begin
-    a =    [10.774290245267503 - 2.01600077393353im   , -0.21512211005762638 + 0.4071609685512763im , -1.1744464421530598 + 0.6046364065537878im , 0.9690771351593747 + 0.24407852288135806im,
-                -0.17679826119222275 - 1.0449912257889253im ,  1.350321850620113 + 0.1195877826052787im , -0.7557518148047799 - 0.809927665736972im  , -0.24869467464627973 + 0.06062801043516876im,
-        -0.83619838577036 + 1.053001604590783im  ,1.0 + 0.0im ]
-
-    A = _BandedMatrix(reverse(a) * Ones{ComplexF64}(1,∞), ∞, length(a)-2, 1)    
-    l,u = bandwidths(A)
-    Q,L = ql(A)
-    @test Q[1:10,1:11] * L[1:11,1:10] ≈ A[1:10,1:10]
-
-    de = tail_de(a)
-    de2 = tail_de_j(a, 2)
-
-    @test !(de ≈ de2)
-
-    X = [transpose(a); 0 transpose(de2)]
-    F = ql_X!(X)
-    @test X[1,1:end-1] ≈ de2
-    factors = _BandedMatrix(Hcat([zero(T); X[1,end-1]; X[2,end-1:-1:1]], [0; X[2,end:-1:1]] * Ones{T}(1,∞)), ∞, l+u, 1)
-    Q2,L2 = QLHessenberg(factors, Fill(F.Q,∞))
-    @test Q2[1:10,1:11] * L2[1:11,1:10] ≈ A[1:10,1:10]
-    @test Q2[1:10,1:11] * (Q2')[1:11,1:10] ≈ I
-    @test (Q2')[1:10,1:100] * (Q2)[1:100,1:10] ≈ I
-    @test (Q')[1:10,1:100] * (Q)[1:100,1:10] ≈ I
-
-    @test !(abs(L2[1,1]) ≈ abs(L[1,1]))
-end