Merge branch 'master' of https://github.com/JuliaMatrices/InfiniteLinearAlgebra.jl

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
+                    @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}}
 
@@ -290,8 +290,13 @@ MemoryLayout(::Type{<:ConstRowMatrix}) = ConstRows()
 MemoryLayout(::Type{<:PertConstRowMatrix}) = PertConstRows()
 bandedcolumns(::ConstRows) = BandedToeplitzLayout()
 bandedcolumns(::PertConstRows) = PertToeplitzLayout()
-subarraylayout(::ConstRows, ::Type{<:Tuple{Any,AbstractInfUnitRange{Int}}}) = ConstRows() # no way to lose const rows
-subarraylayout(::PertConstRows, ::Type{<:Tuple{Any,AbstractInfUnitRange{Int}}}) = PertConstRows() # no way to lose const rows
+subarraylayout(::ConstRows, inds...) = subarraylayout(ApplyLayout{typeof(*)}(), inds...)
+subarraylayout(::PertConstRows, inds...) = subarraylayout(ApplyLayout{typeof(hcat)}(), inds...)
+for Typ in (:ConstRows, :PertConstRows)
+    @eval begin
+        subarraylayout(::$Typ, ::Type{<:Tuple{Any,AbstractInfUnitRange{Int}}}) = $Typ() # no way to lose const rows
+    end
+end
 
 const BandedToeplitzLayout = BandedColumns{ConstRows}
 const PertToeplitzLayout = BandedColumns{PertConstRows}

src/banded/infqltoeplitz.jl

@@ -103,6 +103,13 @@ end
 
 getindex(Q::ProductQ{<:Any,<:Tuple{Vararg{LowerHessenbergQ}}}, i::Integer, j::Integer) = (Q')[j,i]'
 
+
+function (*)(A::ProductQ{T}, x::AbstractVector{S}) where {T,S}
+    TS = promote_op(matprod, T, S)
+    lmul!(A, Base.copymutable(convert(AbstractVector{TS},x)))
+end
+
+
 # LQ where Q is a product of orthogonal operations
 struct QLProduct{T,QQ<:Tuple,LL} <: Factorization{T}
     Qs::QQ

src/infql.jl

@@ -106,7 +106,7 @@ getL(Q::QLHessenberg, ::Tuple{OneToInf{Int},OneToInf{Int}}) where T = LowerTrian
 nzzeros(A::AbstractArray, k) = size(A,k)
 nzzeros(::Zeros, k) = 0
 nzzeros(B::Vcat, k) = sum(size.(B.args[1:end-1],k))
-nzzeros(B::CachedArray, k) = max(size(B.data,k), nzzeros(B.array,k))
+nzzeros(B::CachedArray, k) = max(B.datasize[k], nzzeros(B.array,k))
 function nzzeros(B::AbstractMatrix, k) 
     l,u = bandwidths(B)
     k == 1 ? size(B,2) + l : size(B,1) + u
@@ -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

src/infqr.jl

@@ -91,7 +91,7 @@ end
 function lmul!(A::QRPackedQ{<:Any,<:AdaptiveQRFactors}, B::CachedVector{T,Vector{T},<:Zeros{T,1}}) where T
     require_one_based_indexing(B)
     mA, nA = size(A.factors)
-    sB = length(B.data)
+    sB = B.datasize[1]
     mB = length(B)
     if mA != mB
         throw(DimensionMismatch("matrix A has dimensions ($mA,$nA) but B has dimensions ($mB, $nB)"))
@@ -124,7 +124,7 @@ end
 
 function lmul!(adjA::Adjoint{<:Any,<:QRPackedQ{<:Any,<:AdaptiveQRFactors}}, B::CachedVector{T,Vector{T},<:Zeros{T,1}}) where T
     COLGROWTH = 1000 # rate to grow columns
-    tol = floatmin(T)
+    tol = floatmin(real(T))
 
     require_one_based_indexing(B)
     A = adjA.parent
@@ -164,7 +164,12 @@ function (*)(A::Adjoint{T,<:QRPackedQ{T,<:AdaptiveQRFactors}}, x::AbstractVector
     lmul!(A, Base.copymutable(convert(AbstractVector{TS},x)))
 end
 
-ldiv!(R::UpperTriangular{<:Any,<:AdaptiveQRFactors}, B::AbstractVector) = materialize!(Ldiv(R, B))
+function ldiv!(R::UpperTriangular{<:Any,<:AdaptiveQRFactors}, B::CachedVector{<:Any,<:Any,<:Zeros{<:Any,1}}) 
+    n = B.datasize[1]
+    partialqr!(parent(R).data, n)
+    materialize!(Ldiv(UpperTriangular(view(parent(R).data.data.data,OneTo(n),OneTo(n))), view(B.data,OneTo(n))))
+    B
+end
 
 
 

test/runtests.jl

@@ -7,7 +7,7 @@ import BlockArrays: _BlockArray
 import BlockBandedMatrices: isblockbanded, _BlockBandedMatrix
 import MatrixFactorizations: QLPackedQ
 import BandedMatrices: bandeddata, _BandedMatrix
-import LazyArrays: colsupport, ApplyStyle, MemoryLayout
+import LazyArrays: colsupport, ApplyStyle, MemoryLayout, ApplyLayout
 
 @testset "∞-Toeplitz and Pert-Toeplitz" begin
     A = BandedMatrix(1 => Fill(2im,∞), 2 => Fill(-1,∞), 3 => Fill(2,∞), -2 => Fill(-4,∞), -3 => Fill(-2im,∞))
@@ -47,6 +47,7 @@ end
         @test At[1:10,1:10] ≈ transpose(A)[1:10,1:10] ≈ transpose(A[1:10,1:10])
 
         A = BandedMatrix(-1 => Vcat(Float64[], Fill(1/4,∞)), 0 => Vcat([1.0+im],Fill(0,∞)), 1 => Vcat(Float64[], Fill(1,∞)))
+        @test MemoryLayout(typeof(view(A.data,:,1:10))) == ApplyLayout{typeof(hcat)}()
         Ac = BandedMatrix(A')
         At = BandedMatrix(transpose(A))
         @test Ac[1:10,1:10] ≈ (A')[1:10,1:10] ≈ A[1:10,1:10]'

test/test_infql.jl

@@ -13,6 +13,7 @@ import BandedMatrices: _BandedMatrix
             L = _BandedMatrix(Hcat([e; X[2,2]; X[2,1]], X[2,end:-1:1] * Ones{Float64}(1,∞)), ∞, 2, 0)
             Qn,Ln = ql(A[1:1000,1:1000])
             @test Q[1:10,1:10] ≈ Qn[1:10,1:10]
+            @test Array((Q'A)[1:10,1:10]) ≈ [(Q'A)[k,j] for k=1:10,j=1:10]
             @test (Q'A)[1:10,1:10] ≈ Ln[1:10,1:10] ≈ L[1:10,1:10]
 
             A = BandedMatrix(-1 => Fill(2,∞), 0 => Fill(5+im,∞), 1 => Fill(0.5,∞))
@@ -37,7 +38,7 @@ import BandedMatrices: _BandedMatrix
                 @test Q[1:10,1:11]*L[1:11,1:10] ≈ T[1:10,1:10]
             end
 
-                for (Z,A,B)  in ((2,5.0im,0.5),
+            for (Z,A,B)  in ((2,5.0im,0.5),
                         (2,2.1+0.1im,0.5),
                         (2,2.1-0.1im,0.5),
                         (2,1.5+0.1im,0.5),
@@ -158,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

test/test_infqr.jl

@@ -62,7 +62,8 @@ import InfiniteLinearAlgebra: partialqr!, AdaptiveQRData, AdaptiveLayout
         A = _BandedMatrix(Vcat(Ones(1,∞), (1:∞)', Ones(1,∞)), ∞, 1, 1) 
         Q,R = qr(A);
         b = Vcat([1.,2,3],Zeros(∞))
-        @test lmul!(Q, Base.copymutable(b)).data ≈ qr(A[1:4,1:3]).Q*[1,2,3]
+        @test lmul!(Q, Base.copymutable(b)).datasize[1] == 4
+        @test lmul!(Q, Base.copymutable(b)).data[1:4] ≈ qr(A[1:4,1:3]).Q*[1,2,3]
 
         @test Q[1,1] ≈ -1/sqrt(2)
         @test Q[200_000,200_000] ≈ -1.0
@@ -111,10 +112,12 @@ import InfiniteLinearAlgebra: partialqr!, AdaptiveQRData, AdaptiveLayout
         C = cache(AB);
         resizedata!(C,103,100); resizedata!(C,203,200);
         @test C[103,104] ≈ 1.0
+        C = qr(AB).factors.data.data;
+        resizedata!(C,102,104); resizedata!(C,202,204);
+        @test C[202,204] == AB[202,204]
         F = qr(AB);
-        partialqr!(F.factors.data, 100);
-        partialqr!(F.factors.data, 200);
-        @test norm(F.factors.data.data.data) ≤ 4000
+        partialqr!(F.factors.data, 100); partialqr!(F.factors.data, 200);
+        @test norm(F.factors.data.data.data[Base.OneTo.(F.factors.data.data.datasize)...]) ≤ 4000
         b = Vcat([3,4,5],Zeros(∞))
         @time x = qr(AB) \ b;
         @test x[1:300] ≈ AB[1:300,1:300] \ b[1:300]