Switch to _copy_oftype (#88)

* Change to _copy_oftype

With this, tests still pass on Julia 1.7.2. It doesn't fix 1.8-beta1
yet because there is some other problem too.

* Add StableRNGs

* 1.5 -> 1.6

* Require Julia v1.6

* Update test_muladd.jl

* Fix typo (fix broken test)

* Fix another test by allowing minor differences

* tests pass

Co-authored-by: Sheehan Olver <solver@mac.com>

Author: Rik Huijzer

.github/workflows/ci.yml

@@ -1,7 +1,9 @@
 name: CI
 on:
-  - push
-  - pull_request
+  push:
+    branches:
+      - master
+  pull_request:
 jobs:
   test:
     name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }}
@@ -10,9 +12,9 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.5'
+          - '1.6'
           - '1'
-          - '^1.7.0-0'
+          - '^1.8.0-0'
         os:
           - ubuntu-latest
           - macOS-latest
@@ -62,4 +64,4 @@ jobs:
       - run: julia --project=docs docs/make.jl
         env:
           GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
-          DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }}
+          DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }}

Project.toml

@@ -1,22 +1,23 @@
 name = "ArrayLayouts"
 uuid = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
 authors = ["Sheehan Olver <solver@mac.com>"]
-version = "0.7.10"
+version = "0.8"
 
 [deps]
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [compat]
-FillArrays = "0.12.6, 0.13"
-julia = "1.5"
+FillArrays = "0.13.1"
+julia = "1.6"
 
 [extras]
 Base64 = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
 Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "Random", "Base64", "Compat"]
+test = ["Base64", "Compat", "Random", "StableRNGs", "Test"]

src/ArrayLayouts.jl

@@ -22,7 +22,7 @@ import Base: AbstractArray, AbstractMatrix, AbstractVector,
          AbstractMatrix, AbstractArray, checkindex, unsafe_length, OneTo, one, zero,
         to_shape, _sub2ind, print_matrix, print_matrix_row, print_matrix_vdots,
       checkindex, Slice, @propagate_inbounds, @_propagate_inbounds_meta,
-      _in_range, _range, _rangestyle, Ordered,
+      _in_range, _range, Ordered,
       ArithmeticWraps, floatrange, reverse, unitrange_last,
       AbstractArray, AbstractVector, axes, (:), _sub2ind_recurse, broadcast, promote_eltypeof,
       similar, @_gc_preserve_end, @_gc_preserve_begin,
@@ -35,12 +35,12 @@ import Base.Broadcast: BroadcastStyle, AbstractArrayStyle, Broadcasted, broadcas
                         materialize!, eltypes
 
 import LinearAlgebra: AbstractTriangular, AbstractQ, checksquare, pinv, fill!, tilebufsize, factorize, qr, lu, cholesky,
-                        norm2, norm1, normInf, normMinusInf, qr, lu, qr!, lu!, AdjOrTrans, HermOrSym, copy_oftype,
-                        AdjointAbsVec, TransposeAbsVec, cholcopy, checknonsingular, _apply_ipiv_rows!, ipiv2perm, RealHermSymComplexHerm, chkfullrank
+                        norm2, norm1, normInf, normMinusInf, qr, lu, qr!, lu!, AdjOrTrans, HermOrSym, AdjointAbsVec,
+                        TransposeAbsVec, cholcopy, checknonsingular, _apply_ipiv_rows!, ipiv2perm, RealHermSymComplexHerm, chkfullrank
 
 import LinearAlgebra.BLAS: BlasFloat, BlasReal, BlasComplex
 
-import FillArrays: AbstractFill, getindex_value, axes_print_matrix_row
+import FillArrays: AbstractFill, getindex_value, axes_print_matrix_row, _copy_oftype
 
 import Base: require_one_based_indexing
 
@@ -58,6 +58,14 @@ if VERSION < v"1.7-"
     const RowMaximum = Val{true}
     const NoPivot = Val{false}
 end
+
+if VERSION < v"1.8-"
+    const CRowMaximum = Val{true}
+    const CNoPivot = Val{false}
+else
+    const CRowMaximum = RowMaximum
+    const CNoPivot = NoPivot
+end
 
 
 struct ApplyBroadcastStyle <: BroadcastStyle end
@@ -81,7 +89,9 @@ strides(A::Transpose) = _transpose_strides(strides(parent(A))...)
 """
     ConjPtr{T}
 
-represents that the entry is the complex-conjugate of the pointed to entry.
+A memory address referring to complex conjugated data of type T. However, there is no guarantee
+that the memory is actually valid, or that it actually represents the complex conjugate of data of
+the specified type.
 """
 struct ConjPtr{T}
     ptr::Ptr{T}
@@ -96,6 +106,14 @@ function unsafe_convert(::Type{ConjPtr{T}}, V::SubArray{T,2}) where {T,N,P}
     unsafe_convert(Ptr{T}, view(parent(V)', jr, kr))
 end
 
+Base.elsize(::Type{<:Adjoint{<:Complex,P}}) where P<:AbstractVecOrMat = conjelsize(P)
+conjelsize(::Type{<:Adjoint{<:Complex,P}}) where P<:AbstractVecOrMat = Base.elsize(P)
+conjelsize(::Type{<:Transpose{<:Any, P}}) where {P<:AbstractVecOrMat} = conjelsize(P)
+conjelsize(::Type{<:PermutedDimsArray{<:Any, <:Any, <:Any, <:Any, P}}) where {P} = conjelsize(P)
+conjelsize(::Type{<:ReshapedArray{<:Any,<:Any,P}}) where {P} = conjelsize(P)
+conjelsize(::Type{<:SubArray{<:Any,<:Any,P}}) where {P} = conjelsize(P)
+conjelsize(A::AbstractArray) = conjelsize(typeof(A))
+
 include("memorylayout.jl")
 include("mul.jl")
 include("muladd.jl")
@@ -367,4 +385,4 @@ Base.typed_vcat(::Type{T}, A::LayoutVecOrMats, B::LayoutVecOrMats, C::AbstractVe
 Base.typed_hcat(::Type{T}, A::LayoutVecOrMats, B::LayoutVecOrMats, C::AbstractVecOrMat...) where T = typed_hcat(T, A, B, C...)
 Base.typed_vcat(::Type{T}, A::AbstractVecOrMat, B::LayoutVecOrMats, C::AbstractVecOrMat...) where T = typed_vcat(T, A, B, C...)
 Base.typed_hcat(::Type{T}, A::AbstractVecOrMat, B::LayoutVecOrMats, C::AbstractVecOrMat...) where T = typed_hcat(T, A, B, C...)
-end # module
+end # module

src/diagonal.jl

@@ -56,11 +56,11 @@ copy(M::Lmul{<:DiagonalLayout,<:TridiagonalLayout}) = M.A * convert(Tridiagonal,
 copy(M::Rmul{<:SymTridiagonalLayout,<:DiagonalLayout}) = convert(SymTridiagonal, M.A) * M.B
 copy(M::Lmul{<:DiagonalLayout,<:SymTridiagonalLayout}) = M.A * convert(SymTridiagonal, M.B)
 
-copy(M::Lmul{DiagonalLayout{OnesLayout}}) = copy_oftype(M.B, eltype(M))
-copy(M::Lmul{DiagonalLayout{OnesLayout},<:DiagonalLayout}) = Diagonal(copy_oftype(diagonaldata(M.B), eltype(M)))
-copy(M::Lmul{<:DiagonalLayout,DiagonalLayout{OnesLayout}}) = Diagonal(copy_oftype(diagonaldata(M.A), eltype(M)))
-copy(M::Lmul{DiagonalLayout{OnesLayout},DiagonalLayout{OnesLayout}}) = copy_oftype(M.B, eltype(M))
-copy(M::Rmul{<:Any,DiagonalLayout{OnesLayout}}) = copy_oftype(M.A, eltype(M))
+copy(M::Lmul{DiagonalLayout{OnesLayout}}) = _copy_oftype(M.B, eltype(M))
+copy(M::Lmul{DiagonalLayout{OnesLayout},<:DiagonalLayout}) = Diagonal(_copy_oftype(diagonaldata(M.B), eltype(M)))
+copy(M::Lmul{<:DiagonalLayout,DiagonalLayout{OnesLayout}}) = Diagonal(_copy_oftype(diagonaldata(M.A), eltype(M)))
+copy(M::Lmul{DiagonalLayout{OnesLayout},DiagonalLayout{OnesLayout}}) = _copy_oftype(M.B, eltype(M))
+copy(M::Rmul{<:Any,DiagonalLayout{OnesLayout}}) = _copy_oftype(M.A, eltype(M))
 
 copy(M::Lmul{<:DiagonalLayout{<:AbstractFillLayout}}) = getindex_value(diagonaldata(M.A)) * M.B
 copy(M::Lmul{<:DiagonalLayout{<:AbstractFillLayout},<:DiagonalLayout}) = getindex_value(diagonaldata(M.A)) * M.B
@@ -74,9 +74,9 @@ copy(M::Rmul{<:SymTridiagonalLayout,<:DiagonalLayout{<:AbstractFillLayout}}) = M
 copy(M::Lmul{<:DiagonalLayout{<:AbstractFillLayout},<:SymTridiagonalLayout}) = getindex_value(diagonaldata(M.A)) * M.B
 
 
-copy(M::Rmul{<:BidiagonalLayout,DiagonalLayout{OnesLayout}}) = copy_oftype(M.A, eltype(M))
-copy(M::Lmul{DiagonalLayout{OnesLayout},<:BidiagonalLayout}) =  copy_oftype(M.B, eltype(M))
-copy(M::Rmul{<:TridiagonalLayout,DiagonalLayout{OnesLayout}}) = copy_oftype(M.A, eltype(M))
-copy(M::Lmul{DiagonalLayout{OnesLayout},<:TridiagonalLayout}) =  copy_oftype(M.B, eltype(M))
-copy(M::Rmul{<:SymTridiagonalLayout,DiagonalLayout{OnesLayout}}) = copy_oftype(M.A, eltype(M))
-copy(M::Lmul{DiagonalLayout{OnesLayout},<:SymTridiagonalLayout}) =  copy_oftype(M.B, eltype(M))
+copy(M::Rmul{<:BidiagonalLayout,DiagonalLayout{OnesLayout}}) = _copy_oftype(M.A, eltype(M))
+copy(M::Lmul{DiagonalLayout{OnesLayout},<:BidiagonalLayout}) =  _copy_oftype(M.B, eltype(M))
+copy(M::Rmul{<:TridiagonalLayout,DiagonalLayout{OnesLayout}}) = _copy_oftype(M.A, eltype(M))
+copy(M::Lmul{DiagonalLayout{OnesLayout},<:TridiagonalLayout}) =  _copy_oftype(M.B, eltype(M))
+copy(M::Rmul{<:SymTridiagonalLayout,DiagonalLayout{OnesLayout}}) = _copy_oftype(M.A, eltype(M))
+copy(M::Lmul{DiagonalLayout{OnesLayout},<:SymTridiagonalLayout}) =  _copy_oftype(M.B, eltype(M))

src/factorizations.jl

@@ -305,20 +305,32 @@ _qr!(layout, axes, A, args...; kwds...) = error("Overload _qr!(::$(typeof(layout
 _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...)
 _lu!(layout, axes, A, args...; kwds...) = error("Overload _lu!(::$(typeof(layout)), axes, A)")
-_cholesky(layout, axes, A, ::Val{false}=Val(false); check::Bool = true) = cholesky!(cholcopy(A); check = check)
-_cholesky(layout, axes, A, ::Val{true}; tol = 0.0, check::Bool = true) = cholesky!(cholcopy(A), Val(true); tol = tol, check = check)
+_cholesky(layout, axes, A, ::CNoPivot=CNoPivot(); check::Bool = true) = cholesky!(cholcopy(A); check = check)
+_cholesky(layout, axes, A, ::CRowMaximum; tol = 0.0, check::Bool = true) = cholesky!(cholcopy(A), CRowMaximum(); tol = tol, check = check)
 _factorize(layout, axes, A) = qr(A) # Default to QR
 
 
 _factorize(::AbstractStridedLayout, axes, A) = lu(A)
-function _lu!(::AbstractColumnMajor, axes, A::AbstractMatrix{T}, pivot::Union{NoPivot, RowMaximum} = RowMaximum();
-            check::Bool = true) where T<:BlasFloat
-    if pivot === NoPivot()
-        return generic_lufact!(A, pivot; check = check)
+if VERSION < v"1.8-"
+    function _lu!(::AbstractColumnMajor, axes, A::AbstractMatrix{T}, pivot::Union{NoPivot, RowMaximum} = RowMaximum();
+                check::Bool = true) where T<:BlasFloat
+        if pivot === NoPivot()
+            return generic_lufact!(A, pivot; check = check)
+        end
+        lpt = LAPACK.getrf!(A)
+        check && checknonsingular(lpt[3])
+        return LU{T,typeof(A)}(lpt[1], lpt[2], lpt[3])
+    end
+else
+    function _lu!(::AbstractColumnMajor, axes, A::AbstractMatrix{T}, pivot::Union{NoPivot, RowMaximum} = RowMaximum();
+                check::Bool = true) where T<:BlasFloat
+        if pivot === NoPivot()
+            return generic_lufact!(A, pivot; check = check)
+        end
+        lpt = LAPACK.getrf!(A)
+        check && checknonsingular(lpt[3])
+        return LU{T,typeof(A),typeof(lpt[2])}(lpt[1], lpt[2], lpt[3])
     end
-    lpt = LAPACK.getrf!(A)
-    check && checknonsingular(lpt[3])
-    return LU{T,typeof(A)}(lpt[1], lpt[2], lpt[3])
 end
 
 # for some reason only defined for StridedMatrix in LinearAlgebra
@@ -355,14 +367,14 @@ end
 
 _chol!(_, A, UL) = LinearAlgebra._chol!(A, UL)
 
-function _cholesky!(layout, axes, A::RealHermSymComplexHerm, ::Val{false}; check::Bool = true)
+function _cholesky!(layout, axes, A::RealHermSymComplexHerm, ::CNoPivot; check::Bool = true)
     C, info = _chol!(layout, A.data, A.uplo == 'U' ? UpperTriangular : LowerTriangular)
     check && LinearAlgebra.checkpositivedefinite(info)
     return Cholesky(C.data, A.uplo, info)
 end
 
 function _cholesky!(::SymmetricLayout{<:AbstractColumnMajor}, axes, A::AbstractMatrix{<:BlasReal},
-    ::Val{true}; tol = 0.0, check::Bool = true)
+    ::CRowMaximum; tol = 0.0, check::Bool = true)
     AA, piv, rank, info = LAPACK.pstrf!(A.uplo, A.data, tol)
     C = CholeskyPivoted{eltype(AA),typeof(AA)}(AA, A.uplo, piv, rank, tol, info)
     check && chkfullrank(C)
@@ -390,9 +402,10 @@ end
 
 macro _layoutfactorizations(Typ)
     esc(quote
-        LinearAlgebra.cholesky(A::$Typ, args...; kwds...) = ArrayLayouts._cholesky(ArrayLayouts.MemoryLayout(A), axes(A), A, args...; kwds...)
-        LinearAlgebra.cholesky!(A::LinearAlgebra.RealHermSymComplexHerm{<:Real,<:$Typ}, v::Val{false}=Val(false); check::Bool = true) = ArrayLayouts._cholesky!(ArrayLayouts.MemoryLayout(A), axes(A), A, v; check=check)
-        LinearAlgebra.cholesky!(A::LinearAlgebra.RealHermSymComplexHerm{<:Real,<:$Typ}, v::Val{true}; check::Bool = true, tol = 0.0) = ArrayLayouts._cholesky!(ArrayLayouts.MemoryLayout(A), axes(A), A, v; check=check, tol=tol)
+        LinearAlgebra.cholesky(A::$Typ, v::CNoPivot = CNoPivot(); kwds...) = ArrayLayouts._cholesky(ArrayLayouts.MemoryLayout(A), axes(A), A, v; kwds...)
+        LinearAlgebra.cholesky(A::$Typ, v::CRowMaximum; kwds...) = ArrayLayouts._cholesky(ArrayLayouts.MemoryLayout(A), axes(A), A, v; kwds...)
+        LinearAlgebra.cholesky!(A::LinearAlgebra.RealHermSymComplexHerm{<:Real,<:$Typ}, v::CNoPivot = CNoPivot(); check::Bool = true) = ArrayLayouts._cholesky!(ArrayLayouts.MemoryLayout(A), axes(A), A, v; check=check)
+        LinearAlgebra.cholesky!(A::LinearAlgebra.RealHermSymComplexHerm{<:Real,<:$Typ}, v::CRowMaximum; check::Bool = true, tol = 0.0) = ArrayLayouts._cholesky!(ArrayLayouts.MemoryLayout(A), axes(A), A, v; check=check, tol=tol)
         LinearAlgebra.qr(A::$Typ, args...; kwds...) = ArrayLayouts._qr(ArrayLayouts.MemoryLayout(A), axes(A), A, args...; kwds...)
         LinearAlgebra.qr!(A::$Typ, args...; kwds...) = ArrayLayouts._qr!(ArrayLayouts.MemoryLayout(A), axes(A), A, args...; kwds...)
         LinearAlgebra.lu(A::$Typ, pivot::Union{ArrayLayouts.NoPivot,ArrayLayouts.RowMaximum}; kwds...) = ArrayLayouts._lu(ArrayLayouts.MemoryLayout(A), axes(A), A, pivot; kwds...)
@@ -424,6 +437,3 @@ end
 
 LinearAlgebra.ldiv!(L::LU{<:Any,<:LayoutMatrix}, B::LayoutVector) = ArrayLayouts.ldiv!(L, B)
 
-if VERSION ≥ v"1.7-"
-    @deprecate qr(A::LayoutMatrix, ::Val{true}) qr(A, ColumnNorm())
-end

src/mul.jl

@@ -309,4 +309,4 @@ for Typ in (:LayoutMatrix, :(Symmetric{<:Any,<:LayoutMatrix}), :(Hermitian{<:Any
         @inline -(A::$Typ, Λ::UniformScaling) = _apply(MemoryLayout(A), size(A), -, A, Λ)
         @inline -(Λ::UniformScaling, A::$Typ) = _apply(MemoryLayout(A), size(A), -, Λ, A)
     end
-end
+end

src/muladd.jl

@@ -74,7 +74,7 @@ copy(M::MulAdd) = copyto!(similar(M), M)
 _fill_copyto!(dest, C) = copyto!(dest, C)
 _fill_copyto!(dest, C::Zeros) = zero!(dest) # exploit special fill! overload
 
-@inline copyto!(dest::AbstractArray{T}, M::MulAdd) where T = 
+@inline copyto!(dest::AbstractArray{T}, M::MulAdd) where T =
     muladd!(M.α, unalias(dest,M.A), unalias(dest,M.B), M.β, _fill_copyto!(dest, M.C))
 
 # Modified from LinearAlgebra._generic_matmatmul!
@@ -420,4 +420,4 @@ end
 const ZerosLayouts = Union{ZerosLayout,DualLayout{ZerosLayout}}
 copy(M::MulAdd{<:ZerosLayouts, <:ZerosLayouts, <:ZerosLayouts}) = M.C
 copy(M::MulAdd{<:ZerosLayouts, <:Any, <:ZerosLayouts}) = M.C
-copy(M::MulAdd{<:Any, <:ZerosLayouts, <:ZerosLayouts}) = M.C
+copy(M::MulAdd{<:Any, <:ZerosLayouts, <:ZerosLayouts}) = M.C

test/test_layoutarray.jl

@@ -1,5 +1,5 @@
 using ArrayLayouts, LinearAlgebra, FillArrays, Base64, Test
-import ArrayLayouts: sub_materialize, MemoryLayout, ColumnNorm, RowMaximum
+import ArrayLayouts: sub_materialize, MemoryLayout, ColumnNorm, RowMaximum, CRowMaximum
 
 
 struct MyMatrix <: LayoutMatrix{Float64}
@@ -10,6 +10,7 @@ Base.getindex(A::MyMatrix, k::Int, j::Int) = A.A[k,j]
 Base.setindex!(A::MyMatrix, v, k::Int, j::Int) = setindex!(A.A, v, k, j)
 Base.size(A::MyMatrix) = size(A.A)
 Base.strides(A::MyMatrix) = strides(A.A)
+Base.elsize(::Type{MyMatrix}) = sizeof(Float64)
 Base.unsafe_convert(::Type{Ptr{T}}, A::MyMatrix) where T = Base.unsafe_convert(Ptr{T}, A.A)
 MemoryLayout(::Type{MyMatrix}) = DenseColumnMajor()
 Base.copy(A::MyMatrix) = MyMatrix(copy(A.A))
@@ -22,6 +23,7 @@ Base.getindex(A::MyVector, k::Int) = A.A[k]
 Base.setindex!(A::MyVector, v, k::Int) = setindex!(A.A, v, k)
 Base.size(A::MyVector) = size(A.A)
 Base.strides(A::MyVector) = strides(A.A)
+Base.elsize(::Type{MyVector}) = sizeof(Float64)
 Base.unsafe_convert(::Type{Ptr{T}}, A::MyVector) where T = Base.unsafe_convert(Ptr{T}, A.A)
 MemoryLayout(::Type{MyVector}) = DenseColumnMajor()
 
@@ -100,17 +102,14 @@ MemoryLayout(::Type{MyVector}) = DenseColumnMajor()
 
             @test qr(A) isa LinearAlgebra.QRCompactWY
             @test inv(A) ≈ inv(A.A)
-            if VERSION ≥ v"1.7-"
-                @test qr(A, Val(true)) == qr(A, ColumnNorm())
-            end
 
             S = Symmetric(MyMatrix(reshape(inv.(1:25),5,5) + 10I))
             @test cholesky(S).U ≈ @inferred(cholesky!(deepcopy(S))).U
-            @test cholesky(S,Val(true)).U ≈ cholesky(Matrix(S),Val(true)).U
+            @test cholesky(S,CRowMaximum()).U ≈ cholesky(Matrix(S),CRowMaximum()).U
 
             S = Symmetric(MyMatrix(reshape(inv.(1:25),5,5) + 10I),:L)
             @test cholesky(S).U ≈ @inferred(cholesky!(deepcopy(S))).U
-            @test cholesky(S,Val(true)).U ≈ cholesky(Matrix(S),Val(true)).U
+            @test cholesky(S,CRowMaximum()).U ≈ cholesky(Matrix(S),CRowMaximum()).U
         end
         Bin = randn(5,5)
         B = MyMatrix(copy(Bin))