Add limited memory U and L matrices

Author: Pawel Latawiec

src/lal.jl

@@ -87,9 +87,9 @@ mutable struct LookAheadLanczosDecomp{OpT, OptT, VecT, MatT, ElT, ElRT}
 
     # Eq. 3.9
     # need to keep previous columns of U for G checks
-    U::UpperTriangular{ElT, Matrix{ElT}}
+    U::LimitedMemoryUpperTriangular{ElT, Matrix{ElT}}
     # need to keep previous columns of L for H checks
-    L::UpperHessenberg{ElT, Matrix{ElT}}
+    L::LimitedMemoryUpperHessenberg{ElT, Matrix{ElT}}
 
     # Indices tracking location in block and sequence
     n::Int
@@ -188,8 +188,8 @@ function LookAheadLanczosDecomp(
     Flastrow = Vector{elT}()
     F̃lastcol = Vector{elT}()
 
-    U = UpperTriangular(Matrix{elT}(undef, 0, 0))
-    L = UpperHessenberg(Matrix{elT}(undef, 0, 0))
+    U = LimitedMemoryUpperTriangular{elT, Matrix{elT}}(max_memory)
+    L = LimitedMemoryUpperHessenberg{elT, Matrix{elT}}(max_memory)
 
     # Alg 5.2.0
     n     = 1
@@ -473,12 +473,9 @@ function _update_U!(ld, innerp)
     # U is upper triangular matrix in decomposition of recurrence relation for P-Q sequence
     # updates last column of U
     n, mk, k, kstar = ld.n, ld.mk, ld.k, ld.kstar
-    ld.U  = UpperTriangular(
-        [
-            ld.U fill(0.0, n-1, 1)
-            fill(0.0, 1, n-1) 1.0
-        ]
-    )
+    uvec = fill(0, n)
+    uvec[end] = 1
+    _grow_hcat!(ld.U, uvec)
 
     for i = kstar:k-1
         block_start = mk[i]
@@ -639,17 +636,8 @@ function _update_L!(ld, innerv)
     if !innerv
         Llastcol[nl[l]:end] .= blocks(ld.D)[end] \ ld.Flastcol[nl[l]:end]
     end
-    if isone(n)
-        ld.L = UpperHessenberg(
-            reshape([Llastcol[1]
-                0.0], 2, 1)
-        )
-    else
-        ld.L = UpperHessenberg(
-            [ld.L Llastcol
-            fill(0.0, 1, n)]
-        )
-    end
+    lvec = [Llastcol; 0]
+    _grow_hcat!(ld.L, lvec)
     return ld
 end
 

src/limited_memory_matrices.jl

@@ -4,33 +4,34 @@
     A matrix which keeps only the last `memory` columns. Access outside these columns throws an error. Thus, the underlying storage is a matrix of size `(N, memory)` where `N` is the size of the first dimension
 """
 mutable struct LimitedMemoryMatrix{T, V<:AbstractMatrix{T}} <: AbstractMatrix{T}
-    value::V
+    data::V
     memory::Int
     hsize::Int
 
-    function LimitedMemoryMatrix{T, V}(value::V, memory, hsize) where {T, V<:AbstractMatrix{T}}
-        @assert Base.require_one_based_indexing(value)
-        return new{T, V}(value, memory, hsize)
+    function LimitedMemoryMatrix{T, V}(data::V, memory, hsize) where {T, V<:AbstractMatrix{T}}
+        @assert Base.require_one_based_indexing(data)
+        return new{T, V}(data, memory, hsize)
     end
 end
 
 function LimitedMemoryMatrix(mat::AbstractMatrix{T}, memory) where T
-    value = similar(mat, size(mat, 1), memory)
-    value[:, end-size(mat, 2)+1:end] .= mat
-    return LimitedMemoryMatrix{T, typeof(value)}(value, memory, size(mat, 2))
+    data = similar(mat, size(mat, 1), memory)
+    data[:, end-size(mat, 2)+1:end] .= mat
+    return LimitedMemoryMatrix{T, typeof(data)}(data, memory, size(mat, 2))
 end
 function LimitedMemoryMatrix(vec::AbstractVector{T}, memory) where T
-    value = similar(vec, size(vec, 1), memory)
-    value[:, end] .= vec
-    return LimitedMemoryMatrix{T, typeof(value)}(value, memory, 1)
+    data = similar(vec, size(vec, 1), memory)
+    data[:, end] .= vec
+    return LimitedMemoryMatrix{T, typeof(data)}(data, memory, 1)
 end
 
-Base.size(A::LimitedMemoryMatrix) = (size(A.value, 1), A.hsize)
-Base.similar(A::LimitedMemoryMatrix) = LimitedMemoryMatrix(similar(A.value), A.memory, A.hsize)
-Base.similar(A::LimitedMemoryMatrix, ::Type{T}) where T = LimitedMemoryMatrix(similar(A.value, T), A.memory, A.hsize)
+Base.size(A::LimitedMemoryMatrix) = (size(A.data, 1), A.hsize)
+Base.similar(A::LimitedMemoryMatrix) = LimitedMemoryMatrix(similar(A.data), A.memory, A.hsize)
+Base.similar(A::LimitedMemoryMatrix, ::Type{T}) where T = LimitedMemoryMatrix(similar(A.data, T), A.memory, A.hsize)
 Base.isempty(A::LimitedMemoryMatrix) = size(A, 1) == 0 || size(A, 2) == 0
 
 Base.@propagate_inbounds function Base.setindex!(A::LimitedMemoryMatrix, v, i::Integer, j::Integer)
+    @boundscheck checkbounds(A, i, j)
     lowest_stored_index = size(A, 2) - A.memory + 1
     if j < lowest_stored_index || j > size(A, 2)
         throw(
@@ -40,12 +41,13 @@ Base.@propagate_inbounds function Base.setindex!(A::LimitedMemoryMatrix, v, i::I
         )
     else
         jj = j - lowest_stored_index + 1
-        A.value[i, jj] = v
+        A.data[i, jj] = v
     end
     return v
 end
 
 Base.@propagate_inbounds function Base.getindex(A::LimitedMemoryMatrix, i::Integer, j::Integer)
+    @boundscheck checkbounds(A, i, j)
     lowest_stored_index = size(A, 2) - A.memory + 1
     if j < lowest_stored_index || j > size(A, 2)
         throw(
@@ -54,19 +56,19 @@ Base.@propagate_inbounds function Base.getindex(A::LimitedMemoryMatrix, i::Integ
             )
         )
     else
-        @inbounds v = A.value[i, j - lowest_stored_index + 1]
+        @inbounds v = A.data[i, j - lowest_stored_index + 1]
     end
     return v
 end
 
 function Base.show(io::IO, ::MIME"text/plain", A::LimitedMemoryMatrix) where {T}
     print(io, Base.dims2string(size(A)), " ")
-    Base.showarg(io, A.value, true)
+    Base.showarg(io, A.data, true)
 end
 
 function Base.show(io::IO, A::LimitedMemoryMatrix) where {T}
     print(io, Base.dims2string(size(A)), " ")
-    Base.showarg(io, A.value, true)
+    Base.showarg(io, A.data, true)
 end
 
 """
@@ -78,22 +80,152 @@ size of `A`
 function hcat!(A::LimitedMemoryMatrix, B::AbstractVector)
     @assert size(B, 1) == size(A, 1)
     A.hsize += 1
-    A.value[:, 1:end-1] .= A.value[:, 2:end]
-    A.value[:, end] .= B
+    A.data[:, 1:end-1] .= A.data[:, 2:end]
+    A.data[:, end] .= B
     return A
 end
 
+
 """
-    LimitedMemoryBandedMatrix    
+    LimitedMemoryUpperTriangularMatrix    
 
-A matrix which keeps only the last `memory` columns. Access outside these columns throws an
-  error. Additionally, the matrix is understood to be banded, such that only values
-  `band_offset` from the diagonal are non-zero, where the band is of width `band_width`. Thus, the underlying storage is a matrix of size `(band_width, memory)`
+A matrix which keeps only the last `memory` columns. `setindex!` outside these columns
+throws an error, otherwise entries are set to 0.
 """
-mutable struct LimitedMemoryBandedMatrix{T, V<:AbstractMatrix{T}} <: AbstractMatrix{T}
-    value::V
+mutable struct LimitedMemoryUpperTriangular{T, V<:AbstractMatrix{T}} <: AbstractMatrix{T}
+    data::V
     memory::Int
-    size::Tuple{Int, Int}
-    band_offset::Int
-    band_width::Int
+    hsize::Int
+
+    function LimitedMemoryUpperTriangular{T, V}(data::AbstractMatrix{T}, memory, hsize) where {T, V<:AbstractMatrix{T}}
+        @assert Base.require_one_based_indexing(data)
+        return new{T, typeof(data)}(data, memory, hsize)
+    end
+end
+
+Base.size(A::LimitedMemoryUpperTriangular) = (A.hsize, A.hsize)
+function LimitedMemoryUpperTriangular(mat::AbstractMatrix{T}, memory) where T
+    data = similar(mat, memory, memory)
+    fill!(data, 0)
+    data[:, end-size(mat, 2)+1:end] .= mat
+    return LimitedMemoryUpperTriangular{T, typeof(data)}(data, memory, size(mat, 2))
+end
+function LimitedMemoryUpperTriangular(vec::AbstractVector{T}, memory) where T
+    data = similar(vec, memory, memory)
+    fill!(data, 0)
+    data[:, end] .= vec
+    return LimitedMemoryUpperTriangular{T, typeof(data)}(data, memory, 1)
+end
+function LimitedMemoryUpperTriangular{T, V}(memory) where {T, V<:AbstractMatrix{T}}
+    data = V(undef, memory, memory)
+    fill!(data, 0)
+    return LimitedMemoryUpperTriangular{T, typeof(data)}(data, memory, 0)
+end
+
+Base.@propagate_inbounds function Base.setindex!(A::LimitedMemoryUpperTriangular, v, i::Integer, j::Integer)
+    @boundscheck checkbounds(A, i, j)
+    lowest_stored_index = size(A, 2) - A.memory + 1
+    if j < lowest_stored_index || (i > j) || (j-A.memory) >= i
+        throw(
+            ArgumentError(
+                "Cannot set index ($i, $j) out of memory range, memory kept from $(max(1, lowest_stored_index)) to $(size(A, 2)))"
+            )
+        )
+    else
+        jj = j - lowest_stored_index + 1
+        ii = i - lowest_stored_index + 1 - (size(A, 2) - j)
+        A.data[ii, jj] = v
+    end
+    return v
+end
+
+Base.@propagate_inbounds function Base.getindex(A::LimitedMemoryUpperTriangular, i::Integer, j::Integer)
+    @boundscheck checkbounds(A, i, j)
+    lowest_stored_index_j = size(A, 2) - A.memory + 1
+    if j < lowest_stored_index_j || (i > j) || (j-A.memory) >= i
+        v = zero(eltype(A))
+    else
+        jj = j - lowest_stored_index_j + 1
+        ii = i - lowest_stored_index_j + 1 + (size(A, 2) - j)
+        @inbounds v = A.data[ii, jj]
+    end
+    return v
+end
+
+"""
+    LimitedMemoryUpperHessenberg
+
+A matrix which keeps only the last `memory` columns. `setindex!` outside these columns
+throws an error, otherwise entries are set to 0.
+"""
+mutable struct LimitedMemoryUpperHessenberg{T, V<:AbstractMatrix{T}} <: AbstractMatrix{T}
+    data::V
+    memory::Int
+    hsize::Int
+
+    function LimitedMemoryUpperHessenberg{T, V}(data::AbstractMatrix{T}, memory, hsize) where {T, V<:AbstractMatrix{T}}
+        @assert Base.require_one_based_indexing(data)
+        return new{T, typeof(data)}(data, memory, hsize)
+    end
+end
+
+Base.size(A::LimitedMemoryUpperHessenberg) = (A.hsize+1, A.hsize)
+function LimitedMemoryUpperHessenberg(mat::AbstractMatrix{T}, memory) where T
+    data = similar(mat, memory, memory)
+    fill!(data, 0)
+    data[end-memory+1:end, end-size(mat, 2)+1:end] .= mat
+    return LimitedMemoryUpperHessenberg{T, typeof(data)}(data, memory, size(mat, 2))
+end
+function LimitedMemoryUpperHessenberg(vec::AbstractVector{T}, memory) where T
+    data = similar(vec, memory, memory)
+    fill!(data, 0)
+    data[:, end] .= vec[end-memory+1:end]
+    return LimitedMemoryUpperHessenberg{T, typeof(data)}(data, memory, 1)
+end
+function LimitedMemoryUpperHessenberg{T, V}(memory) where {T, V<:AbstractMatrix{T}}
+    data = V(undef, memory, memory)
+    fill!(data, 0)
+    return LimitedMemoryUpperHessenberg{T, typeof(data)}(data, memory, 0)
+end
+
+Base.@propagate_inbounds function Base.setindex!(A::LimitedMemoryUpperHessenberg, v, i::Integer, j::Integer)
+    @boundscheck checkbounds(A, i, j)
+    lowest_stored_index = size(A, 2) - A.memory + 1
+    if j < lowest_stored_index || (i > j+1) || (j+1-A.memory) >= i
+        throw(
+            ArgumentError(
+                "Cannot set index ($i, $j) out of memory range, memory kept from $(max(1, lowest_stored_index)) to $(size(A, 2)))"
+            )
+        )
+    else
+        jj = j - lowest_stored_index + 1
+        ii = i - lowest_stored_index - (size(A, 2) - j)
+        A.data[ii, jj] = v
+    end
+    return v
+end
+
+Base.@propagate_inbounds function Base.getindex(A::LimitedMemoryUpperHessenberg, i::Integer, j::Integer)
+    @boundscheck checkbounds(A, i, j)
+    lowest_stored_index_j = size(A, 2) - A.memory + 1
+    if j < lowest_stored_index_j || (i > j+1) || (j+1-A.memory) >= i
+        v = zero(eltype(A))
+    else
+        jj = j - lowest_stored_index_j + 1
+        ii = i - lowest_stored_index_j + (size(A, 2) - j)
+        @inbounds v = A.data[ii, jj]
+    end
+    return v
+end
+
+function _grow_hcat!(A, v)
+    if !isempty(A)
+        @assert size(v, 1) == size(A, 1)+1
+    end
+    @assert sum(abs, v[1:end-A.memory]) < eps(real(eltype(A)))
+    mem = min(A.memory, length(v))
+    A.hsize += 1
+    A.data[:, 1:end-1] .= A.data[:, 2:end]
+    A.data[end-mem+1:end, end] .= v[end-mem+1:end]
+    return A    
 end

test/limited_memory_matrices.jl

@@ -25,4 +25,48 @@ using Test
     @test A[3, 3] == 3
     @test A[3, 4] == 4
     @test A[3, 5] == 5
+end
+
+@testset "LimitedMemoryUpperTriangular" begin
+    A = IS.LimitedMemoryUpperTriangular{Float64, Matrix{Float64}}(3)
+    IS._grow_hcat!(A, fill(1.0, 1))
+    IS._grow_hcat!(A, fill(2.0, 2))
+    @test A ≈ [
+        1.0 2.0
+        0.0 2.0
+    ]
+    IS._grow_hcat!(A, fill(3.0, 3))
+    v = fill(4.0, 4)
+    v[1] = 0
+    IS._grow_hcat!(A, v)
+    @test A ≈ [
+        0.0 2.0 3.0 0.0
+        0.0 2.0 3.0 4.0
+        0.0 0.0 3.0 4.0
+        0.0 0.0 0.0 4.0
+    ]
+end
+
+@testset "LimitedMemoryUpperHessenberg" begin
+    A = IS.LimitedMemoryUpperHessenberg{Float64, Matrix{Float64}}(3)
+    IS._grow_hcat!(A, fill(1.0, 2))
+    IS._grow_hcat!(A, fill(2.0, 3))
+    @test A ≈ [
+        1.0 2.0    
+        1.0 2.0
+        0.0 2.0
+    ]
+    v = fill(3.0, 4)
+    v[1] = 0
+    IS._grow_hcat!(A, v)
+    v = fill(4.0, 5)
+    v[1:2] .= 0
+    IS._grow_hcat!(A, v)
+    @test A ≈ [
+        0.0 2.0 0.0 0.0
+        0.0 2.0 3.0 0.0
+        0.0 2.0 3.0 4.0
+        0.0 0.0 3.0 4.0
+        0.0 0.0 0.0 4.0
+    ]
 end