Allow more general *TridiagonalLayout, improve interface (#56)

* Allow more general *TridiagonalLayout, improve interface

* tests pass again

* Hopefully increase coverage

* restore sub/sup/diagonaldata

* special mul

* Require Julia v1.5

* Drop < v"1.5", fix empty range col/rowsupport for bidiagonal

* Update diagonal.jl

Author: dlfivefifty

.github/workflows/ci.yml

@@ -10,8 +10,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.0'
-          - '1'
+          - '1.5'
           - '^1.6.0-0'
         os:
           - ubuntu-latest

Project.toml

@@ -1,7 +1,7 @@
 name = "ArrayLayouts"
 uuid = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.5.4"
+version = "0.6.0"
 
 [deps]
 Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
@@ -12,7 +12,7 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 [compat]
 Compat = "3.16"
 FillArrays = "0.11"
-julia = "1"
+julia = "1.5"
 
 [extras]
 Base64 = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

src/ArrayLayouts.jl

@@ -41,12 +41,7 @@ import LinearAlgebra.BLAS: BlasFloat, BlasReal, BlasComplex
 
 import FillArrays: AbstractFill, getindex_value, axes_print_matrix_row
 
-if VERSION < v"1.2-"
-    import Base: has_offset_axes
-    require_one_based_indexing(A...) = !has_offset_axes(A...) || throw(ArgumentError("offset arrays are not supported but got an array with index other than 1"))
-else
-    import Base: require_one_based_indexing
-end
+import Base: require_one_based_indexing
 
 export materialize, materialize!, MulAdd, muladd!, Ldiv, Rdiv, Lmul, Rmul, Dot,
         lmul, rmul, mul, ldiv, rdiv, mul, MemoryLayout, AbstractStridedLayout,
@@ -208,10 +203,8 @@ copyto!(dest::SubArray{<:Any,2,<:LayoutArray}, src::AdjOrTrans{<:Any,<:LayoutArr
 copyto!(dest::SubArray{<:Any,2,<:LayoutMatrix}, src::SubArray{<:Any,2,<:AdjOrTrans{<:Any,<:LayoutArray}}) = _copyto!(dest, src)
 copyto!(dest::AbstractMatrix, src::SubArray{<:Any,2,<:AdjOrTrans{<:Any,<:LayoutArray}}) = _copyto!(dest, src)
 # ambiguity from sparsematrix.jl
-if VERSION ≥ v"1.5"
-    copyto!(dest::LayoutMatrix, src::SparseArrays.AbstractSparseMatrixCSC) = _copyto!(dest, src)
-    copyto!(dest::SubArray{<:Any,2,<:LayoutMatrix}, src::SparseArrays.AbstractSparseMatrixCSC) = _copyto!(dest, src)
-end
+copyto!(dest::LayoutMatrix, src::SparseArrays.AbstractSparseMatrixCSC) = _copyto!(dest, src)
+copyto!(dest::SubArray{<:Any,2,<:LayoutMatrix}, src::SparseArrays.AbstractSparseMatrixCSC) = _copyto!(dest, src)
 
 # avoid bad copy in Base
 Base.map(::typeof(copy), D::Diagonal{<:LayoutArray}) = Diagonal(map(copy, D.diag))

src/diagonal.jl

@@ -22,7 +22,6 @@ copy(M::Rmul{<:Any,<:DiagonalLayout}) = M.A .* permutedims(diagonaldata(M.B))
 
 
 
-
 # Diagonal multiplication never changes structure
 similar(M::Rmul{<:Any,<:DiagonalLayout}, ::Type{T}, axes) where T = similar(M.A, T, axes)
 # equivalent to rescaling
@@ -45,4 +44,21 @@ copy(M::Ldiv{<:DiagonalLayout{<:AbstractFillLayout},<:DiagonalLayout}) = Diagona
 copy(M::Rdiv{<:DiagonalLayout,<:DiagonalLayout}) = Diagonal(M.A.diag .* inv.(M.B.diag))
 copy(M::Rdiv{<:Any,<:DiagonalLayout}) = M.A .* inv.(permutedims(M.B.diag))
 copy(M::Rdiv{<:Any,<:DiagonalLayout{<:AbstractFillLayout}}) = M.A .* inv(getindex_value(M.B.diag)) 
-copy(M::Rdiv{<:DiagonalLayout,<:DiagonalLayout{<:AbstractFillLayout}}) = Diagonal(M.A.diag .* inv(getindex_value(M.B.diag)))
+copy(M::Rdiv{<:DiagonalLayout,<:DiagonalLayout{<:AbstractFillLayout}}) = Diagonal(M.A.diag .* inv(getindex_value(M.B.diag)))
+
+
+## bi/tridiagonal copy
+copy(M::Rmul{<:BidiagonalLayout,<:DiagonalLayout}) = convert(Bidiagonal, M.A) * M.B
+copy(M::Lmul{<:DiagonalLayout,<:BidiagonalLayout}) = M.A * convert(Bidiagonal, M.B)
+copy(M::Rmul{<:TridiagonalLayout,<:DiagonalLayout}) = convert(Tridiagonal, M.A) * M.B
+copy(M::Lmul{<:DiagonalLayout,<:TridiagonalLayout}) = M.A * convert(Tridiagonal, M.B)
+copy(M::Rmul{<:SymTridiagonalLayout,<:DiagonalLayout}) = convert(SymTridiagonal, M.A) * M.B
+copy(M::Lmul{<:DiagonalLayout,<:SymTridiagonalLayout}) = M.A * convert(SymTridiagonal, M.B)
+
+
+copy(M::Rmul{<:BidiagonalLayout,<:DiagonalLayout{<:AbstractFillLayout}}) = M.A * getindex_value(diagonaldata(M.B))
+copy(M::Lmul{<:DiagonalLayout{<:AbstractFillLayout},<:BidiagonalLayout}) = getindex_value(diagonaldata(M.A)) * M.B
+copy(M::Rmul{<:TridiagonalLayout,<:DiagonalLayout{<:AbstractFillLayout}}) = M.A * getindex_value(diagonaldata(M.B))
+copy(M::Lmul{<:DiagonalLayout{<:AbstractFillLayout},<:TridiagonalLayout}) = getindex_value(diagonaldata(M.A)) * M.B
+copy(M::Rmul{<:SymTridiagonalLayout,<:DiagonalLayout{<:AbstractFillLayout}}) = M.A * getindex_value(diagonaldata(M.B))
+copy(M::Lmul{<:DiagonalLayout{<:AbstractFillLayout},<:SymTridiagonalLayout}) = getindex_value(diagonaldata(M.A)) * M.B

src/ldiv.jl

@@ -95,12 +95,7 @@ __ldiv!(_, F, B) = LinearAlgebra.ldiv!(F, B)
 @inline materialize!(M::Ldiv) = _ldiv!(M.A, M.B)
 @inline materialize!(M::Rdiv) = ldiv!(M.B', M.A')'
 @inline copyto!(dest::AbstractArray, M::Rdiv) = copyto!(dest', Ldiv(M.B', M.A'))'
-
-if VERSION ≥ v"1.1-pre"
-    @inline copyto!(dest::AbstractArray, M::Ldiv) = _ldiv!(dest, M.A, M.B)
-else
-    @inline copyto!(dest::AbstractArray, M::Ldiv) = _ldiv!(dest, M.A, copy(M.B))
-end
+@inline copyto!(dest::AbstractArray, M::Ldiv) = _ldiv!(dest, M.A, copy(M.B))
 
 const MatLdivVec{styleA, styleB, T, V} = Ldiv{styleA, styleB, <:AbstractMatrix{T}, <:AbstractVector{V}}
 const MatLdivMat{styleA, styleB, T, V} = Ldiv{styleA, styleB, <:AbstractMatrix{T}, <:AbstractMatrix{V}}

src/memorylayout.jl

@@ -542,20 +542,24 @@ abstract type AbstractBandedLayout <: MemoryLayout end
 abstract type AbstractTridiagonalLayout <: AbstractBandedLayout end
 
 struct DiagonalLayout{ML} <: AbstractBandedLayout end
-struct BidiagonalLayout{ML} <: AbstractBandedLayout end
-struct SymTridiagonalLayout{ML} <: AbstractTridiagonalLayout end
-struct TridiagonalLayout{ML} <: AbstractTridiagonalLayout end
+struct BidiagonalLayout{DV,EV} <: AbstractBandedLayout end
+struct SymTridiagonalLayout{DV,EV} <: AbstractTridiagonalLayout end
+struct TridiagonalLayout{DL,D,DU} <: AbstractTridiagonalLayout end
 
-bidiagonallayout(_) = BidiagonalLayout{UnknownLayout}()
-tridiagonallayout(_) = TridiagonalLayout{UnknownLayout}()
-symtridiagonallayout(_) = SymTridiagonalLayout{UnknownLayout}()
+bidiagonallayout(dv, ev) = BidiagonalLayout{UnknownLayout,UnknownLayout}()
+tridiagonallayout(dl, d, du) = TridiagonalLayout{UnknownLayout,UnknownLayout,UnknownLayout}()
+
+symtridiagonallayout(d, ev) = SymTridiagonalLayout{UnknownLayout,UnknownLayout}()
+bidiagonallayout(d) = bidiagonallayout(d, d)
+tridiagonallayout(d) = tridiagonallayout(d,d,d)
+symtridiagonallayout(d) = symtridiagonallayout(d,d)
 diagonallayout(_) = DiagonalLayout{UnknownLayout}()
 
 
-diagonallayout(lay::Union{AbstractStridedLayout, AbstractFillLayout}) = DiagonalLayout{typeof(lay)}()
-bidiagonallayout(lay::Union{AbstractStridedLayout, AbstractFillLayout}) = BidiagonalLayout{typeof(lay)}()
-tridiagonallayout(lay::Union{AbstractStridedLayout, AbstractFillLayout}) = TridiagonalLayout{typeof(lay)}()
-symtridiagonallayout(lay::Union{AbstractStridedLayout, AbstractFillLayout}) = SymTridiagonalLayout{typeof(lay)}()
+diagonallayout(::Lay) where Lay<:Union{AbstractStridedLayout, AbstractFillLayout} = DiagonalLayout{Lay}()
+bidiagonallayout(::Lay, ::Lay) where Lay<:Union{AbstractStridedLayout, AbstractFillLayout} = BidiagonalLayout{Lay,Lay}()
+tridiagonallayout(::Lay, ::Lay, ::Lay) where Lay<:Union{AbstractStridedLayout, AbstractFillLayout} = TridiagonalLayout{Lay,Lay,Lay}()
+symtridiagonallayout(::Lay, ::Lay) where Lay<:Union{AbstractStridedLayout, AbstractFillLayout} = SymTridiagonalLayout{Lay,Lay}()
 
 
 MemoryLayout(D::Type{Diagonal{T,P}}) where {T,P} = diagonallayout(MemoryLayout(P))
@@ -586,9 +590,11 @@ transposelayout(ml::SymTridiagonalLayout) = ml
 transposelayout(ml::TridiagonalLayout) = ml
 transposelayout(ml::ConjLayout{DiagonalLayout}) = ml
 
-triangularlayout(::Type{<:TriangularLayout{UPLO,'N'}}, ::TridiagonalLayout{ML}) where {UPLO,ML} = BidiagonalLayout{ML}()
-triangularlayout(::Type{<:TriangularLayout{UPLO,'N'}}, ::TridiagonalLayout{FillLayout}) where UPLO = BidiagonalLayout{FillLayout}()
-triangularlayout(::Type{<:TriangularLayout{UPLO,'U'}}, ::TridiagonalLayout{FillLayout}) where UPLO = BidiagonalLayout{FillLayout}()
+triangularlayout(::Type{<:TriangularLayout{'L','N'}}, ::TridiagonalLayout{DL,D,DU}) where {DL,D,DU} = BidiagonalLayout{D,DL}()
+triangularlayout(::Type{<:TriangularLayout{'U','N'}}, ::TridiagonalLayout{DL,D,DU}) where {UPLO,DL,D,DU} = BidiagonalLayout{D,DU}()
+triangularlayout(::Type{<:TriangularLayout{'L','N'}}, ::TridiagonalLayout{FillLayout,FillLayout,FillLayout}) = BidiagonalLayout{FillLayout,FillLayout}()
+triangularlayout(::Type{<:TriangularLayout{'U','N'}}, ::TridiagonalLayout{FillLayout,FillLayout,FillLayout}) = BidiagonalLayout{FillLayout,FillLayout}()
+triangularlayout(::Type{<:TriangularLayout{UPLO,'U'}}, ::TridiagonalLayout{FillLayout,FillLayout,FillLayout}) where UPLO = BidiagonalLayout{FillLayout,FillLayout}()
 
 bidiagonaluplo(::Union{UpperTriangular,UnitUpperTriangular}) = 'U'
 bidiagonaluplo(::Union{LowerTriangular,UnitLowerTriangular}) = 'L'
@@ -598,20 +604,20 @@ supdiagonaldata(U::Union{UnitUpperTriangular,UpperTriangular}) = supdiagonaldata
 subdiagonaldata(U::Union{UnitLowerTriangular,LowerTriangular}) = subdiagonaldata(triangulardata(U))
 
 adjointlayout(::Type{<:Real}, ml::SymTridiagonalLayout) = ml
-adjointlayout(::Type{<:Real}, ml::TridiagonalLayout) = ml
+adjointlayout(::Type{<:Real}, ::TridiagonalLayout{DL,D,DU}) where {DL,D,DU} = TridiagonalLayout{DU,D,DL}()
 adjointlayout(::Type{<:Real}, ml::BidiagonalLayout) = ml
 
-symmetriclayout(B::BidiagonalLayout{ML}) where ML = SymTridiagonalLayout{ML}()
-hermitianlayout(::Type{<:Real}, B::BidiagonalLayout{ML}) where ML = SymTridiagonalLayout{ML}()
+symmetriclayout(B::BidiagonalLayout{DV,EV}) where {DV,EV} = SymTridiagonalLayout{DV,EV}()
+hermitianlayout(::Type{<:Real}, B::BidiagonalLayout{DV,EV}) where {DV,EV} = SymTridiagonalLayout{DV,EV}()
 hermitianlayout(_, B::BidiagonalLayout) = HermitianLayout{typeof(B)}()
 
+diagonaldata(D::Transpose) = diagonaldata(parent(D))
 subdiagonaldata(D::Transpose) = supdiagonaldata(parent(D))
 supdiagonaldata(D::Transpose) = subdiagonaldata(parent(D))
-diagonaldata(D::Transpose) = diagonaldata(parent(D))
 
+diagonaldata(D::Adjoint{<:Real}) = diagonaldata(parent(D))
 subdiagonaldata(D::Adjoint{<:Real}) = supdiagonaldata(parent(D))
 supdiagonaldata(D::Adjoint{<:Real}) = subdiagonaldata(parent(D))
-diagonaldata(D::Adjoint{<:Real}) = diagonaldata(parent(D))
 
 diagonaldata(S::HermOrSym) = diagonaldata(parent(S))
 subdiagonaldata(S::HermOrSym) = symmetricuplo(S) == 'L' ? subdiagonaldata(parent(S)) : supdiagonaldata(parent(S))
@@ -653,10 +659,14 @@ colsupport(::ZerosLayout, A, _) = 1:0
 rowsupport(::DiagonalLayout, _, k) = isempty(k) ? (1:0) : minimum(k):maximum(k)
 colsupport(::DiagonalLayout, _, j) = isempty(j) ? (1:0) : minimum(j):maximum(j)
 
-colsupport(::BidiagonalLayout, A, j) = 
+function colsupport(::BidiagonalLayout, A, j)
+    isempty(j) && return 1:0
     bidiagonaluplo(A) == 'L' ? (minimum(j):min(size(A,1),maximum(j)+1)) : (max(minimum(j)-1,1):maximum(j))
-rowsupport(::BidiagonalLayout, A, j) = 
+end
+function rowsupport(::BidiagonalLayout, A, j)
+    isempty(j) && return 1:0
     bidiagonaluplo(A) == 'U' ? (minimum(j):min(size(A,2),maximum(j)+1)) : (max(minimum(j)-1,1):maximum(j))
+end
 
 colsupport(::AbstractTridiagonalLayout, A, j) = max(minimum(j)-1,1):min(size(A,1),maximum(j)+1)
 rowsupport(::AbstractTridiagonalLayout, A, j) = max(minimum(j)-1,1):min(size(A,2),maximum(j)+1)

test/test_layoutarray.jl

@@ -171,15 +171,14 @@ MemoryLayout(::Type{MyVector}) = DenseColumnMajor()
         C = randn(ComplexF64,5,5)
         @test ArrayLayouts.lmul!(2, Hermitian(copy(C))) == ArrayLayouts.rmul!(Hermitian(copy(C)), 2) == 2Hermitian(C)
 
-        if VERSION ≥ v"1.5"
-            @test ldiv!(2, deepcopy(b)) == rdiv!(deepcopy(b), 2) == 2\b
-            @test ldiv!(2, deepcopy(A)) == rdiv!(deepcopy(A), 2) == 2\A
-            @test ldiv!(2, deepcopy(A)') == rdiv!(deepcopy(A)', 2) == 2\A'
-            @test ldiv!(2, transpose(deepcopy(A))) == rdiv!(transpose(deepcopy(A)), 2) == 2\transpose(A)
-            @test ldiv!(2, Symmetric(deepcopy(A))) == rdiv!(Symmetric(deepcopy(A)), 2) == 2\Symmetric(A)
-            @test ldiv!(2, Hermitian(deepcopy(A))) == rdiv!(Hermitian(deepcopy(A)), 2) == 2\Hermitian(A)
-            @test ArrayLayouts.ldiv!(2, Hermitian(copy(C))) == ArrayLayouts.rdiv!(Hermitian(copy(C)), 2) == 2\Hermitian(C)
-        end
+        
+        @test ldiv!(2, deepcopy(b)) == rdiv!(deepcopy(b), 2) == 2\b
+        @test ldiv!(2, deepcopy(A)) == rdiv!(deepcopy(A), 2) == 2\A
+        @test ldiv!(2, deepcopy(A)') == rdiv!(deepcopy(A)', 2) == 2\A'
+        @test ldiv!(2, transpose(deepcopy(A))) == rdiv!(transpose(deepcopy(A)), 2) == 2\transpose(A)
+        @test ldiv!(2, Symmetric(deepcopy(A))) == rdiv!(Symmetric(deepcopy(A)), 2) == 2\Symmetric(A)
+        @test ldiv!(2, Hermitian(deepcopy(A))) == rdiv!(Hermitian(deepcopy(A)), 2) == 2\Hermitian(A)
+        @test ArrayLayouts.ldiv!(2, Hermitian(copy(C))) == ArrayLayouts.rdiv!(Hermitian(copy(C)), 2) == 2\Hermitian(C)
     end
 
     @testset "pow/I" begin

test/test_layouts.jl

@@ -139,8 +139,8 @@ struct FooNumber <: Number end
         @testset "Triangular of Tridiagonal" begin
             U = UpperTriangular(T)
             L = LowerTriangular(T)
-            @test MemoryLayout(U) isa BidiagonalLayout{DenseColumnMajor}
-            @test MemoryLayout(L) isa BidiagonalLayout{DenseColumnMajor}
+            @test MemoryLayout(U) isa BidiagonalLayout{DenseColumnMajor,DenseColumnMajor}
+            @test MemoryLayout(L) isa BidiagonalLayout{DenseColumnMajor,DenseColumnMajor}
             @test diagonaldata(U) == diagonaldata(L) == diagonaldata(T)
             @test subdiagonaldata(L) == subdiagonaldata(T)
             @test supdiagonaldata(U) == supdiagonaldata(T)
@@ -304,6 +304,12 @@ struct FooNumber <: Number end
         v = SubArray(Fill(1,10),(1:3,))
         @test ArrayLayouts.sub_materialize(v) ≡ Fill(1,3)
         @test ArrayLayouts._copyto!(Vector{Float64}(undef,3), v) == ones(3)
+
+        T = Tridiagonal(Fill(1,10), Fill(2,11), Fill(3,10))
+        @test MemoryLayout(UpperTriangular(T)) isa BidiagonalLayout{FillLayout,FillLayout}
+        @test MemoryLayout(LowerTriangular(T)) isa BidiagonalLayout{FillLayout,FillLayout}
+        @test MemoryLayout(UnitUpperTriangular(T)) isa BidiagonalLayout{FillLayout,FillLayout}
+        @test MemoryLayout(UnitLowerTriangular(T)) isa BidiagonalLayout{FillLayout,FillLayout}
     end
 
     @testset "Triangular col/rowsupport" begin

test/test_muladd.jl

@@ -1,4 +1,4 @@
-using ArrayLayouts, FillArrays, Random, Test
+using ArrayLayouts, FillArrays, Random, LinearAlgebra, Test
 import ArrayLayouts: DenseColumnMajor, AbstractStridedLayout, AbstractColumnMajor, DiagonalLayout, mul, Mul, zero!
 
 Random.seed!(0)
@@ -669,4 +669,21 @@ Random.seed!(0)
         @test copy(MulAdd(A,B̃)) == A*B̃
         @test eltype(MulAdd(A,B̃)) == eltype(B̃)
     end
+
+    @testset "Bidiagonal" begin
+        BidiagU = Bidiagonal(randn(5), randn(4), :U)
+        BidiagL = Bidiagonal(randn(5), randn(4), :L)
+        Tridiag = Tridiagonal(randn(4), randn(5), randn(4))
+        SymTri = SymTridiagonal(randn(5), randn(4))
+        Diag = Diagonal(randn(5))
+        @test typeof(mul(BidiagU,Diag)) <: Bidiagonal
+        @test typeof(mul(BidiagL,Diag)) <: Bidiagonal
+        @test typeof(mul(Tridiag,Diag)) <: Tridiagonal
+        @test typeof(mul(SymTri,Diag))  <: Tridiagonal
+
+        @test typeof(mul(BidiagU,Diag)) <: Bidiagonal
+        @test typeof(mul(Diag,BidiagL)) <: Bidiagonal
+        @test typeof(mul(Diag,Tridiag)) <: Tridiagonal
+        @test typeof(mul(Diag,SymTri))  <: Tridiagonal
+    end
 end