Add BandedPlusSemiseparable and QR

.github/workflows/ci.yml

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

Project.toml

@@ -1,7 +1,7 @@
 name = "SemiseparableMatrices"
 uuid = "f8ebbe35-cbfb-4060-bf7f-b10e4670cf57"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.4.0"
+version = "0.5.0"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"

src/BandedPlusSemiseparableMatrix.jl

@@ -0,0 +1,46 @@
+struct BandedPlusSemiseparableMatrix{T} <: LayoutMatrix{T}
+    # representing B + tril(UV', -1) + triu(WS', 1)
+    B::BandedMatrix{T, Matrix{T}, Base.OneTo{Int}}
+    U::Matrix{T}
+    V::Matrix{T}
+    W::Matrix{T}
+    S::Matrix{T}
+end
+
+function BandedPlusSemiseparableMatrix(B, (U,V), (W,S))
+    if size(U,1) == size(V,1) == size(W,1) == size(S,1) == size(B,1) == size(B,2) && size(U,2) == size(V,2) && size(W,2) == size(S,2)
+        BandedPlusSemiseparableMatrix(B, U, V, W, S)
+    else
+throw(DimensionMismatch("Dimensions are not compatible."))
+    end
+end
+
+size(A::BandedPlusSemiseparableMatrix) = size(A.B)
+copy(A::BandedPlusSemiseparableMatrix) = A # not mutable
+
+function getindex(A::BandedPlusSemiseparableMatrix, k::Integer, j::Integer)
+    if j > k
+        view(A.W, k, :)' * view(A.S, j, :) + A.B[k,j]
+    elseif k > j
+        view(A.U, k, :)' * view(A.V, j, :) + A.B[k,j]
+    else
+        A.B[k,j]
+    end
+end
+
+function ldiv!(R::UpperTriangular{<:Any,<:BandedPlusSemiseparableMatrix}, b::StridedVector)
+    F = parent(R)
+    n, p = size(F.W)
+    l, m = bandwidths(F.B)
+    T = eltype(F.S)
+    sx = zeros(T, p)
+    for j = n : -1 : 1
+        residual = view(F.W, j, :)' * sx 
+        for k = j+1 : min(j+m, n)
+            residual += F.B[j, k] * b[k]
+        end
+        b[j] = (b[j] - residual) / F.B[j, j]
+        #sx += F.S[j, :] * b[j]
+        mul!(sx, I, view(F.S, j, :), b[j], one(T))
+    end
+end

src/SemiseparableMatrices.jl

@@ -2,16 +2,16 @@ module SemiseparableMatrices
 using LinearAlgebra: BlasFloat
 using ArrayLayouts, BandedMatrices, LazyArrays, LinearAlgebra, MatrixFactorizations, Base
 
-import Base: size, getindex, setindex!, convert, copyto!
+import Base: size, getindex, setindex!, convert, copyto!, copy, axes, getproperty
 import MatrixFactorizations: QR, QRPackedQ, getQ, getR, QRPackedQLayout, AdjQRPackedQLayout
-import LinearAlgebra: qr, qr!, lmul!, ldiv!, rmul!, triu!, factorize, rank
+import LinearAlgebra: qr, qr!, lmul!, ldiv!, rmul!, triu!, factorize, rank, AdjointQ
 import BandedMatrices: _banded_qr!, bandeddata, resize
 import LazyArrays: arguments, applylayout, _cache, CachedArray, CachedMatrix, ApplyLayout, resizedata!, PaddedRows
 import ArrayLayouts: MemoryLayout, sublayout, sub_materialize, MatLdivVec, materialize!, triangularlayout, 
                         triangulardata, zero!, _copyto!, colsupport, rowsupport,
                         _qr, _qr!, _factorize
 
-export SemiseparableMatrix, AlmostBandedMatrix, LowRankMatrix, ApplyMatrix, ApplyArray, almostbandwidths, almostbandedrank
+export SemiseparableMatrix, AlmostBandedMatrix, LowRankMatrix, ApplyMatrix, ApplyArray, almostbandwidths, almostbandedrank, BandedPlusSemiseparableMatrix
 
 LazyArraysBandedMatricesExt = Base.get_extension(LazyArrays, :LazyArraysBandedMatricesExt)
 BandedLayouts = LazyArraysBandedMatricesExt.BandedLayouts
@@ -30,5 +30,7 @@ separablerank(A) = size(arguments(ApplyLayout{typeof(*)}(),A)[1],2)
 include("SemiseparableMatrix.jl")
 include("AlmostBandedMatrix.jl")
 include("invbanded.jl")
+include("BandedPlusSemiseparableMatrix.jl")
+include("semiseparableqr.jl")
 
 end # module