Use a copyto! approach in _BandedMatrix(::BandedRows, A), and switch to PermutedDimsArray in bandedrowsdata (#488)

* Use a copyto! approach, and switch to PermutedDimsArray

* Fix for rectangular matrices

* Remove Module prefix

* Use size instead of sum

* Use axes

* Use Transpose for now

* Change test

---------

Co-authored-by: Sheehan Olver <[email protected]>

Author: DanielVandH

Project.toml

@@ -1,6 +1,6 @@
 name = "BandedMatrices"
 uuid = "aae01518-5342-5314-be14-df237901396f"
-version = "1.10.1"
+version = "1.10.2"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"

src/banded/BandedMatrix.jl

@@ -266,14 +266,21 @@ BandedMatrix(A::AbstractMatrix) = _BandedMatrix(MemoryLayout(A), A)
 _BandedMatrix(::BandedColumns, A::AbstractMatrix) = _BandedMatrix(copy(bandeddata(A)), axes(A,1), bandwidths(A)...)
 function _BandedMatrix(::BandedRows, A::AbstractMatrix)
     bdata = bandedrowsdata(A)
-    m,ν = size(bdata)
-    ℓ,u = bandwidths(A)
-    data = similar(bdata, eltype(bdata), ν, size(A,2))
-    n = size(A, 2)
-    for j in axes(A, 1), i in max(1, j - ℓ):min(n, j + u)
-        data[u + 1 + j - i, i] = bdata[j, ℓ+1+i-j]
+    u, ℓ = bandwidths(A)
+    m, n = size(A)
+    data_new = similar(bdata, eltype(bdata), size(bdata, 2), n)
+    for i in axes(bdata, 2)
+        d = u + 1 - i
+        row_new = ℓ + 1 + d
+        k_start = max(1, 1 - d)
+        k_end = min(n, m - d)
+        if k_start ≤ k_end
+            col_start_A = k_start + d
+            col_end_A = k_end + d
+            copyto!(view(data_new, row_new, k_start:k_end), view(bdata, col_start_A:col_end_A, i))
+        end
     end
-    return _BandedMatrix(data, axes(A, 1), ℓ, u)
+    return _BandedMatrix(data_new, axes(A, 1), u, ℓ)
 end
 function _BandedMatrix(::DiagonalLayout, A::AbstractMatrix{T}) where T
     m,n = size(A)

src/interfaceimpl.jl

@@ -30,9 +30,10 @@ function bandeddata(D::Union{Diagonal, Transpose{<:Any, <:Diagonal}, Adjoint{<:A
     permutedims(diagonaldata(D))
 end
 
-bandedrowsdata(Ac::Transpose) = permutedims(bandeddata(parent(Ac)))
-bandedrowsdata(Ac::Adjoint{<:Real}) = permutedims(bandeddata(parent(Ac)))
-bandedrowsdata(V::SubArray) = permutedims(bandeddata(view(parent(parent(V)), reverse(parentindices(V))...)))
+#TODO: PermutedDimsArray(bandeddata(parent(Ac)), (2, 1))
+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)

test/test_banded.jl

@@ -620,6 +620,9 @@ include("mymatrix.jl")
             @test BandedMatrix(V) == V
         end
 
+        @test BandedMatrices.bandedrowsdata(A') isa Transpose
+        @test BandedMatrices.bandedrowsdata(transpose(A)) isa Transpose
+        @test BandedMatrices.bandedrowsdata(view(A', 1:2, 1:2)) isa Transpose
     end
 end