Support degenerate bands in QRPackedQ (#469)

* Support degenerate bands in QRPackedQ

* Update bandedqr.jl

Author: dlfivefifty

Project.toml

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

src/banded/bandedqr.jl

@@ -48,15 +48,15 @@ function banded_qr_lmul!(A, B)
     l,u = bandwidths(Afactors)
     D = bandeddata(Afactors)
     @inbounds begin
-        for k = min(mA,nA):-1:1
+        for k = min(mA,nA,last(colsupport(B)),length(A.τ)):-1:1
             for j = 1:nB
                 vBj = B[k,j]
-                for i = k+1:min(k+l,mB)
+                for i = k+max(1,-u):min(k+l,mB)
                     vBj += conj(D[i-k+u+1,k])*B[i,j]
                 end
                 vBj = A.τ[k]*vBj
                 B[k,j] -= vBj
-                for i = k+1:min(k+l,mB)
+                for i = k+max(1,-u):min(k+l,mB)
                     B[i,j] -= D[i-k+u+1,k]*vBj
                 end
             end
@@ -81,12 +81,12 @@ function banded_qr_lmul!(adjA::AdjointQType, B)
             cs = colsupport(B,j)
             for k = max(1,minimum(cs)):min(mA,nA,maximum(cs)+l,length(A.τ))
                 vBj = B[k,j]
-                for i = k+1:min(k+l,mB)
+                for i = k+max(1,-u):min(k+l,mB)
                     vBj += conj(D[i-k+u+1,k])*B[i,j]
                 end
                 vBj = conj(A.τ[k])*vBj
                 B[k,j] -= vBj
-                for i = k+1:min(k+l,mB)
+                for i = k+max(1,-u):min(k+l,mB)
                     B[i,j] -= D[i-k+u+1,k]*vBj
                 end
             end

test/test_bandedqr.jl

@@ -142,6 +142,14 @@ Random.seed!(0)
             @test ArrayLayouts.rmul!(copy(B), Q) == rmul!(copy(B), Q) ≈ B*Q
         end
     end
+
+    @testset "lowerbanded B" begin
+        n = 10
+        c = sqrt.(2*(5:2:2n+4) ./ ((3:n+2) .* (4:n+3) ))
+        s = sqrt.(((1:n) .* (2:n+1)) ./ ((3:n+2) .* (4:n+3) ))
+        Q = BandedMatrices.QRPackedQ(BandedMatrix(-2 => -(s ./ (c .+ 1))), c .+ 1)
+        @test Q'Q  ≈ Matrix(Q)'Matrix(Q) ≈ I
+    end
 end
 
 end # module