Try getting degenerate UL working (#46)

* Try getting degenerate UL working

* Update Project.toml

* inf band getindex

* Update infbanded.jl

* Update Project.toml

* Update Project.toml

* Update Project.toml

* v0.4.3

Author: dlfivefifty

Project.toml

@@ -1,6 +1,6 @@
 name = "InfiniteLinearAlgebra"
 uuid = "cde9dba0-b1de-11e9-2c62-0bab9446c55c"
-version = "0.4.2"
+version = "0.4.3"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -16,16 +16,16 @@ MatrixFactorizations = "a3b82374-2e81-5b9e-98ce-41277c0e4c87"
 SemiseparableMatrices = "f8ebbe35-cbfb-4060-bf7f-b10e4670cf57"
 
 [compat]
-ArrayLayouts = "0.4.4"
-BandedMatrices = "0.15.19"
+ArrayLayouts = "0.4.10"
+BandedMatrices = "0.15.23"
 BlockArrays = "0.12.12"
 BlockBandedMatrices = "0.9"
-FillArrays = "0.9.5, 0.10"
-InfiniteArrays = "0.8"
-LazyArrays = "0.17.5, 0.18"
+FillArrays = "0.10"
+InfiniteArrays = "0.8.3"
+LazyArrays = "0.19"
 LazyBandedMatrices = "0.3.2"
 MatrixFactorizations = "0.6"
-SemiseparableMatrices = "0.1.3"
+SemiseparableMatrices = "0.2"
 julia = "1.5"
 
 [extras]

src/banded/infbanded.jl

@@ -1,4 +1,36 @@
 ###
+# BandIndexing
+###
+
+struct InfBandCartesianIndices <: AbstractVector{CartesianIndex{2}}
+    b::Int
+end
+
+InfBandCartesianIndices(b::Band) = InfBandCartesianIndices(b.i)
+
+size(::InfBandCartesianIndices) = (∞,)
+getindex(B::InfBandCartesianIndices, k::Int) = B.b ≥ 0 ? CartesianIndex(k, k+B.b) : CartesianIndex(k-B.b, k)
+
+Base.checkindex(::Type{Bool}, ::NTuple{2,OneToInf{Int}}, ::InfBandCartesianIndices) = true
+BandedMatrices.band_to_indices(_, ::NTuple{2,OneToInf{Int}}, b) = (InfBandCartesianIndices(b),)
+Base.BroadcastStyle(::Type{<:SubArray{<:Any,1,<:Any,Tuple{InfBandCartesianIndices}}}) = LazyArrayStyle{1}()
+
+_inf_banded_sub_materialize(_, V) = V
+function _inf_banded_sub_materialize(::BandedColumns, V)
+    A = parent(V)
+    b = parentindices(V)[1].b
+    data = bandeddata(A)
+    l,u = bandwidths(A)
+    if -l ≤ b ≤ u
+        data[u+1-b, max(1,b+1):end]
+    else
+        Zeros{eltype(V)}(∞) # Not type stable
+    end
+end
+
+sub_materialize(_, V::SubArray{<:Any,1,<:AbstractMatrix,Tuple{InfBandCartesianIndices}}, ::Tuple{InfAxes}) =
+    _inf_banded_sub_materialize(MemoryLayout(parent(V)), V)
+###
 # Diagonal
 ###
 
@@ -423,4 +455,5 @@ function _bidiag_forwardsub!(M::Ldiv{<:Any,<:PaddedLayout})
     end
 
     b_in
-end
+end
+

src/inful.jl

@@ -22,11 +22,14 @@ _ultailL1(c::Number, a::Number, b::Number) = (a + sign(a)*sqrt(a^2-4b*c))/2
 
 function _ultailL1(C::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix)
     d = size(A,1)
-    λs = filter!(λ -> abs(λ) ≤ 1, eigvals([zeros(2,2) -B; -C -A], [B zeros(2,2); zeros(2,2) -B]))
+    λs = filter!(λ -> abs(λ) ≤ 1, eigvals([zeros(d,d) -B; -C -A], [B zeros(d,d); zeros(d,d) -B]))
     @assert length(λs) == d
     V = Matrix{eltype(λs)}(undef, d, d)
-    for (j,λ) in enumerate(λs)
-        V[:,j] = svd(A-λ*B - C/λ).V[:,end] # nullspace(A-λ*B - C/λ)
+    j = 1
+    for λ in union(λs)
+        c = count(==(λ), λs) # multiplicities
+        V[:,j:j+c-1] = svd(A-λ*B - C/λ).V[:,end-c+1:end] # nullspace(A-λ*B - C/λ)
+        j += c
     end
     C*(V*Diagonal(inv.(λs))/V)
 end
@@ -72,4 +75,7 @@ _inf_getL(::BlockTridiagonalToeplitzLayout, F::UL) = mortar(Bidiagonal(F.factors
 
 
 getU(F::UL, ::NTuple{2,Infinity}) = _inf_getU(MemoryLayout(F.factors), F)
-getL(F::UL, ::NTuple{2,Infinity}) = _inf_getL(MemoryLayout(F.factors), F)
+getL(F::UL, ::NTuple{2,Infinity}) = _inf_getL(MemoryLayout(F.factors), F)
+
+getU(F::UL{T,<:Tridiagonal}, ::NTuple{2,Infinity}) where T = _inf_getU(MemoryLayout(F.factors), F)
+getL(F::UL{T,<:Tridiagonal}, ::NTuple{2,Infinity}) where T = _inf_getL(MemoryLayout(F.factors), F)

test/runtests.jl

@@ -1,7 +1,7 @@
 using InfiniteLinearAlgebra, BlockBandedMatrices, BlockArrays, BandedMatrices, InfiniteArrays, FillArrays, LazyArrays, Test, 
         MatrixFactorizations, ArrayLayouts, LinearAlgebra, Random, LazyBandedMatrices
 import InfiniteLinearAlgebra: qltail, toeptail, tailiterate , tailiterate!, tail_de, ql_X!,
-                    InfToeplitz, PertToeplitz, TriToeplitz, InfBandedMatrix, 
+                    InfToeplitz, PertToeplitz, TriToeplitz, InfBandedMatrix, InfBandCartesianIndices,
                     rightasymptotics, QLHessenberg, ConstRows, PertConstRows, BandedToeplitzLayout, PertToeplitzLayout
 import Base: BroadcastStyle
 import BlockArrays: _BlockArray
@@ -12,8 +12,11 @@ import LazyArrays: colsupport, ApplyStyle, MemoryLayout, ApplyLayout, LazyArrayS
 import InfiniteArrays: OneToInf
 import LazyBandedMatrices: BroadcastBandedBlockBandedLayout, BroadcastBandedLayout
 
+
+
 @testset "∞-banded" begin
     D = Diagonal(Fill(2,∞))
+
     B = D[1:∞,2:∞]
     @test B isa BandedMatrix
     @test B[1:10,1:10] == diagm(-1 => Fill(2,9))
@@ -23,6 +26,14 @@ import LazyBandedMatrices: BroadcastBandedBlockBandedLayout, BroadcastBandedLayo
     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 D[band(0)] ≡ Fill(2,∞)
+    @test D[band(1)] ≡ Fill(0,∞)
+    @test A[band(0)][2:10] == 2:10
 end
 
 @testset "∞-block arrays" begin
@@ -133,6 +144,11 @@ end
         @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]
     end
 
     @testset "Fill" begin