Q'b for adaptive QR

Author: dlfivefifty

src/infql.jl

@@ -128,9 +128,10 @@ function lmul!(A::QLPackedQ{<:Any,<:InfBandedMatrix}, B::AbstractVecOrMat)
         for j = 1:nB
             vBj = B[k,j]
             for i = max(1,ν-u):k-1
-                if !iszero(B[i,j])
+                Bij = B[i,j]
+                if !iszero(Bij)
                     allzero = false
-                    vBj += conj(D[i-ν+u+1,ν])*B[i,j]
+                    vBj += conj(D[i-ν+u+1,ν])*Bij
                 end
             end
             vBj = A.τ[k]*vBj

src/infqr.jl

@@ -83,6 +83,10 @@ function qr(A::BandedMatrix{<:Any,<:Any,<:OneToInf})
     QR(AdaptiveQRFactors(data), AdaptiveQRTau(data))
 end
 
+#########
+# lmul!
+#########
+
 
 function lmul!(A::QRPackedQ{<:Any,<:AdaptiveQRFactors}, B::CachedVector{T,Vector{T},<:Zeros{T,1}}) where T
     require_one_based_indexing(B)
@@ -114,4 +118,38 @@ function lmul!(A::QRPackedQ{<:Any,<:AdaptiveQRFactors}, B::CachedVector{T,Vector
         end
     end
     B
+end
+
+
+function lmul!(adjA::Adjoint{<:Any,<:QRPackedQ{<:Any,<:AdaptiveQRFactors}}, B::CachedVector{T,Vector{T},<:Zeros{T,1}}) where T
+    require_one_based_indexing(B)
+    A = adjA.parent
+    mA, nA = size(A.factors)
+    mB = length(B)
+    if mA != mB
+        throw(DimensionMismatch("matrix A has dimensions ($mA,$nA) but B has dimensions ($mB, $nB)"))
+    end
+    Afactors = A.factors
+    sB = length(B.data)
+    @inbounds begin
+        for k = 1:min(mA,nA)
+            vBj = B[k]
+            allzero = k > sB ? iszero(vBj) : false
+            cs = colsupport(Afactors,k)
+            for i = (k+1:mB) ∩ cs
+                Bi = B[i]
+                if !iszero(Bi)
+                    allzero = false
+                    vBj += conj(Afactors[i,k])*Bi
+                end
+            end
+            vBj = conj(A.τ[k])*vBj
+            B[k] -= vBj
+            for i = (k+1:mB) ∩ cs
+                B[i] -= Afactors[i,k]*vBj
+            end
+            allzero && break
+        end
+    end
+    B
 end

test/test_infqr.jl

@@ -67,7 +67,11 @@ import InfiniteLinearAlgebra: partialqr!, AdaptiveQRData, AdaptiveLayout
         @test Q[200_000,200_000] ≈ -1.0
         @test Q[1:101,1:100] ≈ qr(A[1:101,1:100]).Q[:,1:100]
 
-        # lmul!(Q', Base.copymutable(b))
+        r = lmul!(Q', Base.copymutable(b))
+        nr = length(r.data)
+        @test qr(A[1:nr+1,1:nr]).Q'b[1:nr+1] ≈ r[1:nr+1]
+
+        materialize!(Ldiv(R, r))
     end
 end