Support FillArrays.OffDiagonal and remove sublayout piracy for DiagonalLayout (require ArrayLayouts v1.15)

.github/workflows/downstream.yml

@@ -40,6 +40,7 @@ jobs:
         os: [ubuntu-latest]
         package:
           - {repo: BlockBandedMatrices.jl, group: JuliaLinearAlgebra}
+          - {repo: LazyArrays.jl, group: JuliaArrays}
           - {repo: LazyBandedMatrices.jl, group: JuliaLinearAlgebra}
           - {repo: ApproxFunBase.jl, group: JuliaApproximation}
           - {repo: ApproxFun.jl, group: JuliaApproximation}

Project.toml

@@ -1,6 +1,6 @@
 name = "BandedMatrices"
 uuid = "aae01518-5342-5314-be14-df237901396f"
-version = "1.13.1"
+version = "1.14.0"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -18,11 +18,11 @@ CliqueTreesExt = "CliqueTrees"
 
 [compat]
 Aqua = "0.8"
-ArrayLayouts = "1.6.0"
+ArrayLayouts = "1.15"
 BlockArrays = "1"
 CliqueTrees = "0.5, 1"
 Documenter = "1"
-FillArrays = "1.3"
+FillArrays = "1.18"
 GenericLinearAlgebra = "0.3, 0.4"
 InfiniteArrays = "0.12, 0.13, 0.14, 0.15"
 LinearAlgebra = "1"

src/BandedMatrices.jl

@@ -32,9 +32,10 @@ import ArrayLayouts: AbstractTridiagonalLayout, BidiagonalLayout, BlasMatLdivVec
                      hermitianlayout, materialize, materialize!, reflector!, reflectorApply!, rowsupport,
                      sub_materialize, subdiagonaldata, sublayout, supdiagonaldata, symmetricdata, symmetriclayout,
                      symmetricuplo, transposelayout, triangulardata, triangularlayout, zero!,
-                     QRPackedQLayout, AdjQRPackedQLayout, LayoutVecOrMats
+                     QRPackedQLayout, AdjQRPackedQLayout, LayoutVecOrMats, OffDiagonalLayout, offdiagonaldata,
+                     offdiagonalindex
 
-import FillArrays: AbstractFill, getindex_value, _broadcasted_zeros, unique_value, OneElement, RectDiagonal, OneElementMatrix, OneElementVector, ZerosMatrix, ZerosVector
+import FillArrays: AbstractFill, getindex_value, _broadcasted_zeros, unique_value, OneElement, RectDiagonal, OffDiagonal, OneElementMatrix, OneElementVector, ZerosMatrix, ZerosVector
 
 const libblas = LinearAlgebra.BLAS.libblas
 const liblapack = LinearAlgebra.BLAS.liblapack

src/banded/BandedMatrix.jl

@@ -290,33 +290,40 @@ function _BandedMatrix(::BandedRows, A::AbstractMatrix)
     return _BandedMatrix(data_new, axes(A, 1), u, ℓ)
 end
 function _BandedMatrix(::DiagonalLayout, A::AbstractMatrix{T}) where T
-    m,n = size(A)
-    dat = Matrix{T}(undef, 1, n)
-    copyto!(vec(dat), diagonaldata(A))
-    _BandedMatrix(dat, m, 0, 0)
+    ret = BandedMatrix{T}(undef, axes(A), (0,0))
+    copyto!(ret.data, diagonaldata(A))
+    ret
+end
+function _BandedMatrix(::OffDiagonalLayout, A::AbstractMatrix{T}) where T
+    k = offdiagonalindex(A)
+    band = offdiagonaldata(A)
+    ret = BandedMatrix{T}(undef, axes(A), (-k,k))
+    copyto!(view(ret.data, range(max(1,k+1); length=length(band))), band)
+    ret
 end
 function _BandedMatrix(::BidiagonalLayout, A::AbstractMatrix{T}) where T
     m,n = size(A)
-    dat = Matrix{T}(undef, 2, n)
     if bidiagonaluplo(A) == 'L'
-        copyto!(view(dat, 1, :), diagonaldata(A))
-        copyto!(view(dat, 2, 1:n-1), subdiagonaldata(A))
-        _BandedMatrix(dat, m, 1, 0)
+        ret = BandedMatrix{T}(undef, axes(A), (1,0))
+        copyto!(view(ret.data, 1, :), diagonaldata(A))
+        copyto!(view(ret.data, 2, 1:n-1), subdiagonaldata(A))
+        ret
     else # u
-        copyto!(view(dat, 1, 2:n), supdiagonaldata(A))
-        copyto!(view(dat, 2, :), diagonaldata(A))
-        _BandedMatrix(dat, m, 0, 1)
+        ret = BandedMatrix{T}(undef, axes(A), (0,1))
+        copyto!(view(ret.data, 1, 2:n), supdiagonaldata(A))
+        copyto!(view(ret.data, 2, :), diagonaldata(A))
+        ret
     end
 end
 
 function _BandedMatrix(::AbstractTridiagonalLayout, A::AbstractMatrix{T}) where T
-    m,n = size(A)
-    dat = Matrix{T}(undef, 3, n)
-
-    copyto!(view(dat, 1, 2:n), supdiagonaldata(A))
-    copyto!(view(dat, 2, :), diagonaldata(A))
-    copyto!(view(dat, 3, 1:n-1), subdiagonaldata(A))
-    _BandedMatrix(dat, m, 1, 1)
+    ret = BandedMatrix{T}(undef, axes(A), (1,1))
+    d = diagonaldata(A)
+    n = length(d)
+    copyto!(view(ret.data, 1, 2:n), supdiagonaldata(A))
+    copyto!(view(ret.data, 2, :), d)
+    copyto!(view(ret.data, 3, 1:n-1), subdiagonaldata(A))
+    ret
 end
 
 

src/interfaceimpl.jl

@@ -4,7 +4,7 @@
 # Sparse BroadcastStyle
 ##
 
-for Typ in (:Diagonal, :SymTridiagonal, :Tridiagonal, :Bidiagonal)
+for Typ in (:Diagonal, :SymTridiagonal, :Tridiagonal, :Bidiagonal, :OffDiagonal)
     @eval begin
         BroadcastStyle(::StructuredMatrixStyle{<:$Typ}, ::BandedStyle) =
             BandedStyle()
@@ -35,11 +35,18 @@ bandedrowsdata(Ac::Transpose{<:Number}) = Transpose(bandeddata(parent(Ac)))
 bandedrowsdata(Ac::Adjoint{<:Real}) = Transpose(bandeddata(parent(Ac)))
 bandedrowsdata(V::SubArray{<:Number}) = Transpose(bandeddata(view(parent(parent(V)), reverse(parentindices(V))...)))
 
-# treat subinds as banded
-sublayout(::DiagonalLayout{L}, inds::Type) where L = sublayout(bandedcolumns(L()), inds)
-sublayout(::DiagonalLayout{L}, inds::Type{<:NTuple{2,AbstractUnitRange{Int}}}) where L = sublayout(bandedcolumns(L()), inds)
-
-# bandeddata(V::SubArray{<:Any,2,<:Diagonal}) = view(bandeddata(parent(V)), :, parentindices(V)[2])
+# OffDiagonal, and views of Diagonal/OffDiagonal by unit ranges, are supported on a single band
+function bandwidths(::OffDiagonalLayout, A)
+    b = offdiagonalindex(A)
+    m,n = size(A)
+    if -b ≥ m || b ≥ n
+        (-720, -720) # empty
+    else
+        (-b, b)
+    end
+end
+inbands_getindex(A::OffDiagonal, k::Integer, j::Integer) = A.band[min(k,j)]
+inbands_setindex!(A::OffDiagonal, v, k::Integer, j::Integer) = (A.band[min(k,j)] = v)
 
 isbanded(::Bidiagonal) = true
 bandwidths(A::Bidiagonal) = A.uplo == 'U' ? (0,1) : (1,0)

test/test_interface.jl

@@ -53,7 +53,7 @@ LinearAlgebra.fill!(A::PseudoBandedMatrix, v) = fill!(A.data,v)
 
         V = view(Eye(5),2:3,:)
         @test bandwidths(V) == (-1,1)
-        @test MemoryLayout(typeof(V)) isa BandedColumns{OnesLayout}
+        @test MemoryLayout(typeof(V)) isa OffDiagonalLayout{OnesLayout}
         @test BandedMatrix(V) == [0 1 0 0 0; 0 0 1 0 0]
 
         B = brand(5,5,1,2)
@@ -79,7 +79,9 @@ LinearAlgebra.fill!(A::PseudoBandedMatrix, v) = fill!(A.data,v)
 
         @test bandeddata(A) == bandeddata(Adjoint(A)) == bandeddata(Transpose(A)) == diag(A)'
 
-        @test MemoryLayout(view(A, 1:3,2:4)) isa BandedColumns{DenseColumnMajor}
+        @test MemoryLayout(view(A, 1:3,2:4)) isa OffDiagonalLayout{DenseColumnMajor}
+        @test bandwidths(view(A, 1:3,2:4)) == (1,-1)
+        @test BandedMatrix(view(A, 1:3,2:4)) == A[1:3,2:4]
         @test MemoryLayout(view(A, [1,2,3],2:4)) isa UnknownLayout
 
         A[band(0)][1] = 3
@@ -99,6 +101,44 @@ LinearAlgebra.fill!(A::PseudoBandedMatrix, v) = fill!(A.data,v)
         @test @inferred(B[Band(0)]) == 1:∞
     end
 
+    @testset "OffDiagonal" begin
+        for A in (FillArrays.OffDiagonal([1.,2,3], 1), FillArrays.OffDiagonal([1.,2,3], -2), FillArrays.OffDiagonal(Fill(2.0,3), 1, 3, 5),
+                  FillArrays.OffDiagonal([1.,2,3], 1)', view(Diagonal([1.,2,3,4,5]), 2:5, 1:3))
+            k = ArrayLayouts.offdiagonalindex(A)
+            @test isbanded(A)
+            @test bandwidths(A) == (-k,k)
+            @test BandedMatrix(A) isa BandedMatrix
+            @test BandedMatrix(A) == A
+            @test bandwidths(BandedMatrix(A)) == (-k,k)
+            B = brand(size(A,2), 4, 1, 1)
+            C = brand(4, size(A,1), 1, 1)
+            @test A*B ≈ Matrix(A)*B
+            @test C*A ≈ C*Matrix(A)
+        end
+
+        A = FillArrays.OffDiagonal([1.,2,3], 1)
+        @test MemoryLayout(A) isa OffDiagonalLayout{DenseColumnMajor}
+        @test BandedMatrices.inbands_getindex(A, 2, 3) == 2
+        BandedMatrices.inbands_setindex!(A, 4, 2, 3)
+        @test A[2,3] == 4
+        @test A[band(1)] == [1,4,3]
+        @test A[band(0)] == zeros(4)
+        @test MemoryLayout(view(A, 1:3, 2:4)) isa OffDiagonalLayout{DenseColumnMajor}
+        @test bandwidths(view(A, 1:3, 2:4)) == (0,0)
+
+        @testset "broadcast with BandedMatrix" begin
+            A = FillArrays.OffDiagonal([1.,2,3], 1)
+            B = brand(4, 4, 1, 1)
+            for C in (A .+ B, B .+ A, A .* B, sin.(A) .+ 2 .* B)
+                @test C isa BandedMatrix
+            end
+            @test A .+ B == Matrix(A) .+ Matrix(B)
+            @test bandwidths(A .+ B) == (1,1)
+            @test A .* B == Matrix(A) .* Matrix(B)
+            @test sin.(A) .+ 2 .* B ≈ sin.(Matrix(A)) .+ 2 .* Matrix(B)
+        end
+    end
+
     @testset "SymTridiagonal" begin
         A = SymTridiagonal([1,2,3],[4,5])
         @test @inferred(A[Band(0)]) == [1,2,3]
@@ -303,7 +343,8 @@ end
     @test bandwidths(A) == bandwidths(Ã) == (0,2)
     @test à == A == Matrix(B) + Matrix(Bu)
 
-    @test layout_getindex(D,1:10,1:10) isa BandedMatrix
+    @test layout_getindex(D,1:10,1:10) isa FillArrays.OffDiagonal
+    @test layout_getindex(D,1:10,1:10) == D
     @test layout_getindex(B,1:10,1:10) isa BandedMatrix
     @test layout_getindex(Bu,1:10,1:10) isa BandedMatrix
     @test layout_getindex(J,1:10,1:10) isa BandedMatrix