Add extra step in _qr to allow dispatch on length, Bidiagonal solve (#26)

* Add extra step in _qr to allow dispatch on length (for InfiniteLinearAlgebra.jl)

* Add _mul(layoutA, layoutB, A, B), Bidiagonal substitution

* v0.3.6

* Update ArrayLayouts.jl

* Add Bidiagonal tests

* Update runtests.jl

Author: dlfivefifty

Project.toml

@@ -1,7 +1,7 @@
 name = "ArrayLayouts"
 uuid = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.3.5"
+version = "0.3.6"
 
 [deps]
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"

src/ArrayLayouts.jl

@@ -249,6 +249,7 @@ Base.print_matrix_row(io::IO,
         LayoutVector,
         AbstractTriangular{<:Any,<:LayoutMatrix},
         AdjOrTrans{<:Any,<:LayoutMatrix},
+        AdjOrTrans{<:Any,<:LayoutVector},
         HermOrSym{<:Any,<:LayoutMatrix},
         SubArray{<:Any,2,<:LayoutMatrix},
         Diagonal{<:Any,<:LayoutVector}}, A::Vector,

src/factorizations.jl

@@ -274,7 +274,8 @@ function materialize!(M::Rmul{<:Any,<:AdjQRPackedQLayout})
 end
 
 
-_qr(layout, axes, A; kwds...) = Base.invoke(qr, Tuple{AbstractMatrix{eltype(A)}}, A; kwds...)
+__qr(layout, lengths, A; kwds...) = Base.invoke(qr, Tuple{AbstractMatrix{eltype(A)}}, A; kwds...)
+_qr(layout, axes, A; kwds...) = __qr(layout, map(length, axes), A; kwds...)
 _qr(layout, axes, A, pivot::P; kwds...) where P = Base.invoke(qr, Tuple{AbstractMatrix{eltype(A)},P}, A, pivot; kwds...)
 _lu(layout, axes, A; kwds...) = Base.invoke(lu, Tuple{AbstractMatrix{eltype(A)}}, A; kwds...)
 _lu(layout, axes, A, pivot::P; kwds...) where P = Base.invoke(lu, Tuple{AbstractMatrix{eltype(A)},P}, A, pivot; kwds...)

src/ldiv.jl

@@ -109,28 +109,35 @@ const BlasMatRdivMat{styleA, styleB, T<:BlasFloat} = MatRdivMat{styleA, styleB,
 
 macro _layoutldiv(Typ)
     esc(quote
-        LinearAlgebra.ldiv!(A::$Typ, x::AbstractVector) = ArrayLayouts.materialize!(ArrayLayouts.Ldiv(A,x))
-        LinearAlgebra.ldiv!(A::$Typ, x::AbstractMatrix) = ArrayLayouts.materialize!(ArrayLayouts.Ldiv(A,x))
-        LinearAlgebra.ldiv!(A::$Typ, x::StridedVector) = ArrayLayouts.materialize!(ArrayLayouts.Ldiv(A,x))
-        LinearAlgebra.ldiv!(A::$Typ, x::StridedMatrix) = ArrayLayouts.materialize!(ArrayLayouts.Ldiv(A,x))
+        LinearAlgebra.ldiv!(A::$Typ, x::AbstractVector) = ArrayLayouts.ldiv!(A,x)
+        LinearAlgebra.ldiv!(A::$Typ, x::AbstractMatrix) = ArrayLayouts.ldiv!(A,x)
+        LinearAlgebra.ldiv!(A::$Typ, x::StridedVector) = ArrayLayouts.ldiv!(A,x)
+        LinearAlgebra.ldiv!(A::$Typ, x::StridedMatrix) = ArrayLayouts.ldiv!(A,x)
 
-        LinearAlgebra.ldiv!(A::Factorization, x::$Typ) = ArrayLayouts.materialize!(ArrayLayouts.Ldiv(A,x))
+        LinearAlgebra.ldiv!(A::Factorization, x::$Typ) = ArrayLayouts.ldiv!(A,x)
 
-        Base.:\(A::$Typ, x::AbstractVector) = ArrayLayouts.materialize(ArrayLayouts.Ldiv(A,x))
-        Base.:\(A::$Typ, x::AbstractMatrix) = ArrayLayouts.materialize(ArrayLayouts.Ldiv(A,x))
+        Base.:\(A::$Typ, x::AbstractVector) = ArrayLayouts.ldiv(A,x)
+        Base.:\(A::$Typ, x::AbstractMatrix) = ArrayLayouts.ldiv(A,x)
 
-        Base.:\(x::AbstractMatrix, A::$Typ) = ArrayLayouts.materialize(ArrayLayouts.Ldiv(x,A))
-        Base.:\(x::Diagonal, A::$Typ) = ArrayLayouts.materialize(ArrayLayouts.Ldiv(x,A))
+        Base.:\(x::AbstractMatrix, A::$Typ) = ArrayLayouts.ldiv(x,A)
+        Base.:\(x::Diagonal, A::$Typ) = ArrayLayouts.ldiv(x,A)
 
-        Base.:\(x::$Typ, A::$Typ) = ArrayLayouts.materialize(ArrayLayouts.Ldiv(x,A))
+        Base.:\(A::Bidiagonal{<:Number}, B::$Typ{<:Number}) = ArrayLayouts.ldiv(A,B)
+        Base.:\(A::Bidiagonal, B::$Typ) = ArrayLayouts.ldiv(A,B)
+        Base.:\(transA::Transpose{<:Number,<:Bidiagonal{<:Number}}, B::$Typ{<:Number}) = ArrayLayouts.ldiv(transA,B)
+        Base.:\(transA::Transpose{<:Any,<:Bidiagonal}, B::$Typ) = ArrayLayouts.ldiv(transA,B)
+        Base.:\(adjA::Adjoint{<:Number,<:Bidiagonal{<:Number}}, B::$Typ{<:Number}) = ArrayLayouts.ldiv(adjA,B)
+        Base.:\(adjA::Adjoint{<:Any,<:Bidiagonal}, B::$Typ) = ArrayLayouts.ldiv(adjA,B)
 
-        Base.:/(A::$Typ, x::AbstractVector) = ArrayLayouts.materialize(ArrayLayouts.Rdiv(A,x))
-        Base.:/(A::$Typ, x::AbstractMatrix) = ArrayLayouts.materialize(ArrayLayouts.Rdiv(A,x))
+        Base.:\(x::$Typ, A::$Typ) = ArrayLayouts.ldiv(x,A)
 
-        Base.:/(x::AbstractMatrix, A::$Typ) = ArrayLayouts.materialize(ArrayLayouts.Rdiv(x,A))
-        Base.:/(x::Diagonal, A::$Typ) = ArrayLayouts.materialize(ArrayLayouts.Rdiv(x,A))
+        Base.:/(A::$Typ, x::AbstractVector) = ArrayLayouts.rdiv(A,x)
+        Base.:/(A::$Typ, x::AbstractMatrix) = ArrayLayouts.rdiv(A,x)
 
-        Base.:/(x::$Typ, A::$Typ) = ArrayLayouts.materialize(ArrayLayouts.Rdiv(x,A))
+        Base.:/(x::AbstractMatrix, A::$Typ) = ArrayLayouts.rdiv(x,A)
+        Base.:/(x::Diagonal, A::$Typ) = ArrayLayouts.rdiv(x,A)
+
+        Base.:/(x::$Typ, A::$Typ) = ArrayLayouts.rdiv(x,A)
     end)
 end
 
@@ -148,3 +155,5 @@ macro layoutldiv(Typ)
         ArrayLayouts.@_layoutldiv UnitLowerTriangular{T, <:SubArray{T,2,<:$Typ{T}}} where T
     end)
 end
+
+@_layoutldiv LayoutVector

src/muladd.jl

@@ -444,7 +444,8 @@ function MulAdd(A::AbstractArray{T}, B::AbstractMatrix{V}) where {T,V}
     MulAdd(scalarone(TV), A, B, scalarzero(TV), fillzeros(TV,(axes(A,1),axes(B,2))))
 end
 
-mul(A::AbstractArray, B::AbstractArray) = materialize(MulAdd(A,B))
+@inline _mul(layoutA, layoutB, A, B) = materialize(MulAdd(A,B))
+@inline mul(A::AbstractArray, B::AbstractArray) = _mul(MemoryLayout(A), MemoryLayout(B), A, B)
 
 macro layoutmul(Typ)
     ret = quote

src/triangular.jl

@@ -116,7 +116,7 @@ materialize!(M::MatRmulMat{<:AbstractStridedLayout,<:TriangularLayout}) = Linear
 ########
 
 
-@inline function copyto!(dest::AbstractArray, M::Ldiv{<:TriangularLayout})
+@inline function copyto!(dest::AbstractArray, M::Ldiv{<:Union{TriangularLayout,BidiagonalLayout}})
     A, B = M.A, M.B
     dest ≡ B || copyto!(dest, B)
     ldiv!(A, dest)
@@ -149,7 +149,7 @@ for UNIT in ('U', 'N')
     end
 end
 
-function materialize!(M::MatLdivMat{<:TriangularLayout})
+function materialize!(M::MatLdivMat{<:Union{TriangularLayout,BidiagonalLayout}})
     A,X = M.A,M.B
     size(A,2) == size(X,1) || thow(DimensionMismatch("Dimensions must match"))
     @views for j in axes(X,2)
@@ -229,3 +229,63 @@ function materialize!(M::MatLdivVec{<:TriangularLayout{'L','U'}})
     end
     b
 end
+
+function _bidiag_backsub!(M)
+    A,b = M.A, M.B
+    N = last(colsupport(b,1))
+    dv = diagonaldata(A)
+    ev = supdiagonaldata(A)
+    b[N] = bj1 = dv[N]\b[N]
+    
+    @inbounds for j = (N - 1):-1:1
+        bj  = b[j]
+        bj -= ev[j] * bj1
+        dvj = dv[j]
+        if iszero(dvj)
+            throw(SingularEbception(j))
+        end
+        bj   = dvj\bj
+        b[j] = bj1 = bj
+    end
+    
+    b
+end
+
+function _bidiag_forwardsub!(M)
+    A, b = M.A, M.B
+    dv = diagonaldata(A)
+    ev = subdiagonaldata(A)
+    N = length(b)
+    b[1] = bj1 = dv[1]\b[1]
+    @inbounds for j = 2:N
+        bj  = b[j]
+        bj -= ev[j - 1] * bj1
+        dvj = dv[j]
+        if iszero(dvj)
+            throw(SingularEbception(j))
+        end
+        bj   = dvj\bj
+        b[j] = bj1 = bj
+    end
+    b
+end
+
+#Generic solver using naive substitution, based on LinearAlgebra/src/bidiag.jl
+function materialize!(M::MatLdivVec{<:BidiagonalLayout}) 
+    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 A, $N, does not match one of the length of b, $(length(b))"))
+    end
+
+    if N == 0
+        return b
+    end
+
+    if bidiagonaluplo(A) == 'L' #do forward substitution
+        _bidiag_forwardsub!(M)
+    else #do backward substitution
+        _bidiag_backsub!(M)
+    end
+end

test/runtests.jl

@@ -72,19 +72,45 @@ MemoryLayout(::Type{MyVector}) = DenseColumnMajor()
         muladd!(1.0, A, A, 2.0, B)
         @test all(B .=== A.A^2 + 2Bin)
 
-        # tiled_blasmul!
-        B = MyMatrix(copy(Bin))
-        muladd!(1.0, Ones(5,5), A, 2.0, B)
+        @testset "tiled_blasmul!" begin
+            B = MyMatrix(copy(Bin))
+            muladd!(1.0, Ones(5,5), A, 2.0, B)
+        end
 
-        #generic_blasmul!
-        A = BigFloat.(randn(5,5))
-        Bin = BigFloat.(randn(5,5))
-        B = copy(Bin)
-        muladd!(1.0, Ones(5,5), A, 2.0, B)
-        @test B == Ones(5,5)*A + 2.0Bin
+        @testset "generic_blasmul!" begin
+            A = BigFloat.(randn(5,5))
+            Bin = BigFloat.(randn(5,5))
+            B = copy(Bin)
+            muladd!(1.0, Ones(5,5), A, 2.0, B)
+            @test B == Ones(5,5)*A + 2.0Bin
+        end
 
         C = MyMatrix([1 2; 3 4])
         @test stringmime("text/plain", C) == "2×2 MyMatrix:\n 1.0  2.0\n 3.0  4.0"
+
+        @testset "layoutldiv" begin
+            A = MyMatrix(randn(5,5))
+            x = randn(5)
+            X = randn(5,5)
+            t = view(randn(10),[1,3,4,6,7])
+            T = view(randn(10,5),[1,3,4,6,7],:)
+            t̃ = copy(t)
+            T̃ = copy(T)
+            B = Bidiagonal(randn(5),randn(4),:U)
+            @test ldiv!(A, copy(x)) ≈ A\x
+            @test A\t ≈ A\t̃
+            # QR is not general enough
+            @test_broken ldiv!(A, t) ≈ A\t
+            @test ldiv!(A, copy(X)) ≈ A\X
+            @test A\T ≈ A\T̃
+            @test_broken A/T ≈ A/T̃
+            @test_broken ldiv!(A, T) ≈ A\T
+            @test B\A ≈ B\Matrix(A)
+            @test transpose(B)\A ≈ transpose(B)\Matrix(A)
+            @test B'\A ≈ B'\Matrix(A)
+            @test A\A ≈ I
+            @test_broken A/A ≈ I
+        end
     end
 end
 

test/test_ldiv.jl

@@ -233,5 +233,21 @@ import ArrayLayouts: ApplyBroadcastStyle, QRCompactWYQLayout, QRCompactWYLayout,
             end
         end
     end
+
+    @testset "Bidiagonal" begin
+        n = 10
+        L = Bidiagonal(randn(n), randn(n-1), :L)
+        U = Bidiagonal(randn(n), randn(n-1), :U)
+        b = randn(n)
+        B = randn(n,2)
+
+        @test ArrayLayouts.ldiv(L,b) == L \ b
+        @test ArrayLayouts.ldiv(U,b) == U \ b
+        @test ArrayLayouts.ldiv(L,B) == L \ B
+        @test ArrayLayouts.ldiv(U,B) == U \ B
+
+        @test_throws DimensionMismatch ArrayLayouts.ldiv(L,randn(3))
+        @test_throws DimensionMismatch materialize!(Ldiv(L,randn(3)))
+    end
 end