solves with QL

Author: dlfivefifty

src/InfiniteLinearAlgebra.jl

@@ -13,7 +13,7 @@ import LinearAlgebra: lmul!,  ldiv!, matprod, qr, AbstractTriangular, AbstractQ,
 import LazyArrays: CachedArray, CachedMatrix, CachedVector, DenseColumnMajor, FillLayout, ApplyMatrix, check_mul_axes, ApplyStyle, LazyArrayApplyStyle, LazyArrayStyle,
                     CachedMatrix, CachedArray, resizedata!, MemoryLayout, mulapplystyle, LmulStyle, RmulStyle,
                     colsupport, rowsupport, triangularlayout, factorize, subarraylayout, sub_materialize,
-                    @lazymul, ApplyLayout
+                    @lazymul, ApplyLayout, TriangularLayout, PaddedLayout, materialize!, MatLdivVec, triangulardata
 import MatrixFactorizations: ql, ql!, QLPackedQ, getL, getR, reflector!, reflectorApply!, QL, QR, QRPackedQ
 
 import BlockArrays: BlockSizes, cumulsizes, _find_block, AbstractBlockVecOrMat, sizes_from_blocks

src/banded/infbanded.jl

@@ -1,6 +1,6 @@
 const TriToeplitz{T} = Tridiagonal{T,Fill{T,1,Tuple{OneToInf{Int}}}}
 const ConstRowMatrix{T} = ApplyMatrix{T,typeof(*),<:Tuple{<:AbstractVector,<:AbstractFill{<:Any,2,Tuple{OneTo{Int},OneToInf{Int}}}}}
-const PertConstRowMatrix{T} = Hcat{T,<:Tuple{Matrix{T},<:ConstRowMatrix{T}}}
+const PertConstRowMatrix{T} = Hcat{T,<:Tuple{Array{T},<:ConstRowMatrix{T}}}
 const InfToeplitz{T} = BandedMatrix{T,<:ConstRowMatrix{T},OneToInf{Int}}
 const PertToeplitz{T} = BandedMatrix{T,<:PertConstRowMatrix{T},OneToInf{Int}}
 

src/infql.jl

@@ -276,4 +276,28 @@ end
 function (*)(A::Adjoint{T,<:QLPackedQ{T,<:InfBlockBandedMatrix}}, x::AbstractVector{S}) where {T,S}
     TS = promote_op(matprod, T, S)
     lmul!(A, cache(convert(AbstractVector{TS},x)))
-end
+end
+
+ldiv!(F::QLProduct, b::AbstractVector) = ldiv!(F.L, lmul!(F.Q',b))
+
+function materialize!(M::MatLdivVec{<:TriangularLayout{'L','N',BandedColumns{PertConstRows}},<:PaddedLayout})
+    A,b = M.A,M.B
+    require_one_based_indexing(A, b)
+    n = size(A, 2)
+    if !(n == length(b))
+        throw(DimensionMismatch("second dimension of left hand side A, $n, and length of right hand side b, $(length(b)), must be equal"))
+    end
+    data = triangulardata(A)
+    nz = nzzeros(b,1)
+    @inbounds for j in 1:n
+        iszero(data[j,j]) && throw(SingularException(j))
+        bj = b[j] = data[j,j] \ b[j]
+        allzero = j > nz && iszero(bj)
+        for i in (j+1:n) ∩ colsupport(data,j)
+            b[i] -= data[i,j] * bj
+            allzero = allzero && iszero(b[i])
+        end
+        allzero && break
+    end
+    b
+end

test/test_infql.jl

@@ -159,4 +159,9 @@ import BandedMatrices: _BandedMatrix
         @test Qn[1:10,1:10] * diagm(0 => [Ones(5); -(-1).^(1:5)]) ≈ Q[1:10,1:10]
         @test diagm(0 => [Ones(5); -(-1).^(1:5)]) * Ln[1:10,1:10] ≈ L[1:10,1:10]
     end
+
+    @testset "solve with QL" begin
+        A = BandedMatrix(-1 => Fill(2,∞), 0 => Fill(5,∞), 1 => Fill(0.5,∞))
+        @test (qr(A)\Vcat(1.0,Zeros(∞)))[1:1000] ≈ (ql(A)\Vcat(1.0,Zeros(∞)))[1:1000]
+    end
 end