Support for compact operators (#88)

* Support for compact operators

* Support general padded

* don't use .data for PaddedLayout

* Don't grow cols so fast

* v0.5.14

* v0.5.15

* chop! results to avoid too much data

* Update infbanded.jl

* Update Project.toml

* fix resizing

* Drop Julia v1.5

* increase coverage

* fix resizedata_chop!

* Update runtests.jl

Author: dlfivefifty

.github/workflows/ci.yml

@@ -10,8 +10,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.5'
-          - '1'
+          - '1.6'
         os:
           - ubuntu-latest
           - macOS-latest

Project.toml

@@ -1,6 +1,6 @@
 name = "InfiniteLinearAlgebra"
 uuid = "cde9dba0-b1de-11e9-2c62-0bab9446c55c"
-version = "0.5.14"
+version = "0.6.0"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -17,18 +17,18 @@ MatrixFactorizations = "a3b82374-2e81-5b9e-98ce-41277c0e4c87"
 SemiseparableMatrices = "f8ebbe35-cbfb-4060-bf7f-b10e4670cf57"
 
 [compat]
-ArrayLayouts = "0.7.2"
+ArrayLayouts = "0.7.5"
 BandedMatrices = "0.16.9"
-BlockArrays = "0.15, 0.16"
-BlockBandedMatrices = "0.10"
-DSP = "0.6, 0.7"
+BlockArrays = "0.16"
+BlockBandedMatrices = "0.11"
+DSP = "0.7"
 FillArrays = "0.11, 0.12"
-InfiniteArrays = "0.11, 0.12"
+InfiniteArrays = "0.12"
 LazyArrays = "0.21.8"
-LazyBandedMatrices = "0.6"
+LazyBandedMatrices = "0.7"
 MatrixFactorizations = "0.8"
 SemiseparableMatrices = "0.2.7"
-julia = "1.5"
+julia = "1.6"
 
 [extras]
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

examples/blocktridiagonal.jl

@@ -1,9 +1,8 @@
-using Revise, InfiniteLinearAlgebra, BlockBandedMatrices, BandedMatrices, BlockArrays, InfiniteArrays, FillArrays, LazyArrays, Test
+using InfiniteLinearAlgebra, BlockBandedMatrices, BandedMatrices, BlockArrays, InfiniteArrays, FillArrays, LazyArrays, Test
 
 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), ∞)))
 
 
-BlockBandedMatrix(A)
-isb
+                 

src/InfiniteLinearAlgebra.jl

@@ -59,6 +59,36 @@ function ArrayLayouts._power_by_squaring(_, ::NTuple{2,InfiniteCardinal{0}}, A::
     end
 end
 
+
+function chop!(c::AbstractVector, tol::Real)
+    @assert tol >= 0
+
+    @inbounds for k=length(c):-1:1
+        if abs(c[k]) > tol
+            resize!(c,k)
+            return c
+        end
+    end
+    resize!(c,0)
+    c
+end
+
+function chop(A::AbstractMatrix, tol)
+    for k = size(A,1):-1:1
+        if norm(view(A,k,:))>tol
+            A=A[1:k,:]
+            break
+        end
+    end
+    for j = size(A,2):-1:1
+        if norm(view(A,:,j))>tol
+            A=A[:,1:j]
+            break
+        end
+    end
+    return A
+end
+
 export Vcat, Fill, ql, ql!, ∞, ContinuousSpectrumError, BlockTridiagonal
 
 include("infconv.jl")

src/banded/infbanded.jl

@@ -367,6 +367,10 @@ _BandedMatrix(::PertToeplitzLayout, A::AbstractMatrix) =
 
 sub_materialize(_, V, ::Tuple{<:BlockedUnitRange{<:InfRanges}}) = V
 sub_materialize(::AbstractBlockLayout, V, ::Tuple{<:BlockedUnitRange{<:InfRanges}}) = V
+function sub_materialize(::PaddedLayout, v::AbstractVector{T}, ax::Tuple{<:BlockedUnitRange{<:InfRanges}}) where T
+    dat = paddeddata(v)
+    PseudoBlockVector(Vcat(sub_materialize(dat), Zeros{T}(∞)), ax)
+end
 
 ##
 # UniformScaling
@@ -428,12 +432,15 @@ _bandedfill_mul(M::MulAdd, ::Tuple{Any,InfAxes}, ::Tuple{InfAxes,InfAxes}) = App
 _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{<:BandedColumns{<:AbstractFillLayout}, <:PertToeplitzLayout}) = ApplyArray(M)
+mulreduce(M::Mul{BandedToeplitzLayout, BandedToeplitzLayout}) = ApplyArray(M)
+mulreduce(M::Mul{BandedToeplitzLayout}) = ApplyArray(M)
+mulreduce(M::Mul{<:Any, BandedToeplitzLayout}) = ApplyArray(M)
+mulreduce(M::Mul{<:BandedColumns{<:AbstractFillLayout}, PertToeplitzLayout}) = ApplyArray(M)
 mulreduce(M::Mul{<:PertToeplitzLayout, <:BandedColumns{<:AbstractFillLayout}}) = ApplyArray(M)
-mulreduce(M::Mul{<:BandedColumns{<:AbstractFillLayout}, <:BandedToeplitzLayout}) = ApplyArray(M)
-mulreduce(M::Mul{<:BandedToeplitzLayout, <:BandedColumns{<:AbstractFillLayout}}) = ApplyArray(M)
-mulreduce(M::Mul{<:AbstractQLayout, <:BandedToeplitzLayout}) = ApplyArray(M)
-mulreduce(M::Mul{<:AbstractQLayout, <:PertToeplitzLayout}) = ApplyArray(M)
+mulreduce(M::Mul{<:BandedColumns{<:AbstractFillLayout}, BandedToeplitzLayout}) = ApplyArray(M)
+mulreduce(M::Mul{BandedToeplitzLayout, <:BandedColumns{<:AbstractFillLayout}}) = ApplyArray(M)
+mulreduce(M::Mul{<:AbstractQLayout, BandedToeplitzLayout}) = ApplyArray(M)
+mulreduce(M::Mul{<:AbstractQLayout, PertToeplitzLayout}) = ApplyArray(M)
 
 
 function _bidiag_forwardsub!(M::Ldiv{<:Any,<:PaddedLayout})

src/infqr.jl

@@ -164,7 +164,7 @@ partialqr!(F::AdaptiveQRFactors, n) = partialqr!(F.data, n)
 _view_QRPackedQ(A, kr, jr) = QRPackedQ(view(A.factors.data.data.data,kr,jr), view(A.τ.data.τ,jr))
 function materialize!(M::MatLmulVec{<:QRPackedQLayout{<:AdaptiveLayout},<:PaddedLayout})
     A,B = M.A,M.B
-    sB = B.datasize[1]
+    sB = size(paddeddata(B),1)
     partialqr!(A.factors.data,sB)
     jr = oneto(sB)
     m = maximum(colsupport(A,jr))
@@ -175,6 +175,39 @@ function materialize!(M::MatLmulVec{<:QRPackedQLayout{<:AdaptiveLayout},<:Padded
     B
 end
 
+function resizedata_chop!(v::CachedVector, tol)
+    c = paddeddata(v)
+    n = length(c)
+    k_tol = n
+    for k = n:-1:1
+        if abs(c[k]) > tol
+            v.datasize = (k_tol,)
+            return v
+        end
+    end
+    v.datasize = (0,)
+    v
+end
+
+function resizedata_chop!(v::PseudoBlockVector, tol)
+    c = paddeddata(v.blocks)
+    n = length(c)
+    k_tol = n
+    for k = n:-1:1
+        if abs(c[k]) > tol
+            k_tol = k
+            break
+        end
+    end
+    ax = axes(v,1)
+    K = findblock(ax,k_tol)
+    n2 = last(ax[K])
+    resize!(v.blocks.data, n2)
+    zero!(view(v.blocks.data, n+1:n2))
+    v.blocks.datasize = (n2,)
+    v
+end
+
 _norm(x::Number) = abs(x)
 
 function materialize!(M::MatLmulVec{<:AdjQRPackedQLayout{<:AdaptiveLayout},<:PaddedLayout})
@@ -190,7 +223,8 @@ function materialize!(M::MatLmulVec{<:AdjQRPackedQLayout{<:AdaptiveLayout},<:Pad
     if mA != mB
         throw(DimensionMismatch("matrix A has dimensions ($mA,$nA) but B has dimensions ($mB, $nB)"))
     end
-    sB = B.datasize[1]
+    Bdata = paddeddata(B)
+    sB = size(Bdata,1)
     l,u = bandwidths(A.factors)
     if l == 0 # diagonal special case
         return B
@@ -205,16 +239,17 @@ function materialize!(M::MatLmulVec{<:AdjQRPackedQLayout{<:AdaptiveLayout},<:Pad
             cs_max = maximum(cs)
             kr = j:cs_max
             resizedata!(B, min(cs_max,mB))
-            if (j > sB && maximum(_norm,view(B.data,j:last(colsupport(A.factors,j)))) ≤ tol)
+            Bdata = paddeddata(B)
+            if (j > sB && maximum(_norm,view(Bdata,j:last(colsupport(A.factors,j)))) ≤ tol)
                 break
             end
             partialqr!(A.factors.data, min(cs_max,nA))
             Q_N = _view_QRPackedQ(A, kr, jr)
-            lmul!(Q_N', view(B.data, kr))
+            lmul!(Q_N', view(Bdata, kr))
             jr = last(jr)+1:min(last(jr)+COLGROWTH,nA)
         end
     end
-    B
+    resizedata_chop!(B, tol)
 end
 
 function resizedata!(B::PseudoBlockVector, M::Block{1})
@@ -230,10 +265,10 @@ end
 
 function materialize!(M::MatLmulVec{<:QRPackedQLayout{<:AdaptiveLayout{<:AbstractBlockBandedLayout}},<:PaddedLayout})
     A,B_in = M.A,M.B
-    sB = B_in.datasize[1]
+    sB = length(paddeddata(B_in))
     ax1,ax2 = axes(A.factors.data.data)
     B = PseudoBlockVector(B_in, (ax2,))
-    SB = findblock(ax2, length(B_in.data))
+    SB = findblock(ax2, sB)
     partialqr!(A.factors.data,SB)
     JR = Block(1):SB
     M = maximum(blockcolsupport(A.factors,JR))
@@ -248,12 +283,12 @@ function materialize!(M::MatLmulVec{<:AdjQRPackedQLayout{<:AdaptiveLayout{<:Abst
     adjA,B_in = M.A,M.B
     A = adjA.parent
     T = eltype(M)
-    COLGROWTH = 1000 # rate to grow columns
+    COLGROWTH = 300 # rate to grow columns
     tol = 1E-30
     ax1 = axes(A.factors.data.data,1)
     B = PseudoBlockVector(B_in, (ax1,))
 
-    SB = findblock(ax1, length(B_in.data))
+    SB = findblock(ax1, length(paddeddata(B_in)))
     MA, NA = blocksize(A.factors.data.data.array)
     JR = Block(1):findblock(ax1,COLGROWTH)
 
@@ -278,7 +313,7 @@ function materialize!(M::MatLmulVec{<:AdjQRPackedQLayout{<:AdaptiveLayout{<:Abst
             JR = last(JR)+1:findblock(ax1,last(jr)+COLGROWTH)
         end
     end
-    B
+    resizedata_chop!(B, tol)
 end
 
 

test/runtests.jl

@@ -13,6 +13,16 @@ import LazyArrays: colsupport, MemoryLayout, ApplyLayout, LazyArrayStyle, argume
 import InfiniteArrays: OneToInf, oneto, RealInfinity
 import LazyBandedMatrices: BroadcastBandedBlockBandedLayout, BroadcastBandedLayout, LazyBandedLayout
 
+@testset "chop" begin
+    a = randn(5)
+    b = [a; zeros(5)]
+    InfiniteLinearAlgebra.chop!(b, eps())
+    @test b == a
+
+    A = randn(5,5)
+    @test InfiniteLinearAlgebra.chop([A zeros(5,2); zeros(2,5) zeros(2,2)],eps()) == A
+end
+
 include("test_infconv.jl")
 include("test_infbanded.jl")
 

test/test_compact.jl

@@ -0,0 +1,11 @@
+using InfiniteLinearAlgebra, FillArrays, LazyArrays, ArrayLayouts, Test
+
+@testset "Compact" begin
+    A = ApplyMatrix(hvcat, 2, randn(5,5), Zeros(5,∞), Zeros(∞,5), Zeros(∞,∞))
+    b = [randn(10); zeros(∞)];
+    @test ((I + A) \ b)[1:10] ≈ (I+A)[1:10,1:10] \ b[1:10]
+
+    C = zeros(∞,∞);
+    C[1:5,1:5] .= randn.()
+    @test_skip ((I + C) \ b)[1:10] ≈ (I+C)[1:10,1:10] \ b[1:10]
+end

test/test_infbanded.jl

@@ -41,6 +41,9 @@ import BandedMatrices: _BandedMatrix
         @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 MemoryLayout(Tridiagonal(Fill(1,∞), Fill(2,∞), Fill(3,∞))) isa TridiagonalToeplitzLayout
         @test MemoryLayout(Bidiagonal(Fill(1,∞), Fill(2,∞), :U)) isa BidiagonalToeplitzLayout
@@ -64,6 +67,14 @@ import BandedMatrices: _BandedMatrix
             @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]
+        end
     end
 
     @testset "Pert-Toeplitz" begin