Code clean-up

.github/workflows/ci.yml

@@ -14,6 +14,7 @@ jobs:
         version:
           - '1.6'
           - '1'
+          - 'nightly'
         os:
           - ubuntu-latest
           - macOS-latest

Project.toml

@@ -14,10 +14,9 @@ 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 = ["Base64", "Compat", "Random", "StableRNGs", "Test"]
+test = ["Base64", "Random", "StableRNGs", "Test"]

src/ArrayLayouts.jl

@@ -1,57 +1,39 @@
 module ArrayLayouts
 using Base: _typed_hcat
 using Base, Base.Broadcast, LinearAlgebra, FillArrays, SparseArrays
-import LinearAlgebra.BLAS
-
-import Base: AbstractArray, AbstractMatrix, AbstractVector,
-        ReinterpretArray, ReshapedArray, AbstractCartesianIndex, Slice,
-             RangeIndex, BroadcastStyle, copyto!, length, broadcastable, axes,
-             getindex, eltype, tail, IndexStyle, IndexLinear, getproperty,
-             *, +, -, /, \, ==, isinf, isfinite, sign, angle, show, isless,
-         fld, cld, div, min, max, minimum, maximum, mod,
-         <, ≤, >, ≥, promote_rule, convert, copy,
-         size, step, isempty, length, first, last, ndims,
-         getindex, setindex!, intersect, @_inline_meta, inv,
-         sort, sort!, issorted, sortperm, diff, cumsum, sum, in, broadcast,
-         eltype, parent, real, imag,
-         conj, transpose, adjoint, permutedims, vec,
-         exp, log, sqrt, cos, sin, tan, csc, sec, cot,
-                   cosh, sinh, tanh, csch, sech, coth,
-                   acos, asin, atan, acsc, asec, acot,
-                   acosh, asinh, atanh, acsch, asech, acoth, (:),
-         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, Ordered,
-      ArithmeticWraps, floatrange, reverse, unitrange_last,
-      AbstractArray, AbstractVector, axes, (:), _sub2ind_recurse, broadcast, promote_eltypeof,
-      similar, @_gc_preserve_end, @_gc_preserve_begin,
-      @nexprs, @ncall, @ntuple, tuple_type_tail,
-      all, any, isbitsunion, issubset, replace_in_print_matrix, replace_with_centered_mark,
-      strides, unsafe_convert, first_index, unalias, union
-
-import Base.Broadcast: BroadcastStyle, AbstractArrayStyle, Broadcasted, broadcasted,
-                        combine_eltypes, DefaultArrayStyle, instantiate, materialize,
-                        materialize!, eltypes
-
-import LinearAlgebra: AbstractTriangular, AbstractQ, QRCompactWYQ, QRPackedQ, checksquare, pinv,
-                        fill!, tilebufsize, factorize, qr, lu, cholesky,
-                        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, _copy_oftype
-
-import Base: require_one_based_indexing
+using LinearAlgebra.BLAS
+
+using Base: AbstractCartesianIndex, OneTo, RangeIndex, ReinterpretArray, ReshapedArray,
+            Slice, tuple_type_tail, unalias,
+            @propagate_inbounds, @_propagate_inbounds_meta
+
+import Base: axes, size, length, eltype, ndims, first, last, diff, isempty, union, sort!,
+                ==, *, +, -, /, \, copy, copyto!, similar, getproperty, getindex, strides,
+                reverse, unsafe_convert
+
+using Base.Broadcast: Broadcasted
+
+import Base.Broadcast: BroadcastStyle, broadcastable, instantiate, materialize, materialize!
+
+using LinearAlgebra: AbstractTriangular, AbstractQ, QRCompactWYQ, QRPackedQ, checksquare,
+                        pinv, tilebufsize, cholcopy,
+                        norm2, norm1, normInf, normMinusInf,
+                        AdjOrTrans, HermOrSym, RealHermSymComplexHerm, AdjointAbsVec, TransposeAbsVec,
+                        checknonsingular, _apply_ipiv_rows!, ipiv2perm, chkfullrank
+
+using LinearAlgebra.BLAS: BlasFloat, BlasReal, BlasComplex
+
+using FillArrays: AbstractFill, getindex_value, axes_print_matrix_row, _copy_oftype
+
+using Base: require_one_based_indexing
 
 export materialize, materialize!, MulAdd, muladd!, Ldiv, Rdiv, Lmul, Rmul, Dot,
-        lmul, rmul, mul, ldiv, rdiv, mul, MemoryLayout, AbstractStridedLayout,
+        lmul, mul, ldiv, rdiv, mul, MemoryLayout, AbstractStridedLayout,
         DenseColumnMajor, ColumnMajor, ZerosLayout, FillLayout, AbstractColumnMajor, RowMajor, AbstractRowMajor, UnitStride,
         DiagonalLayout, ScalarLayout, SymTridiagonalLayout, TridiagonalLayout, BidiagonalLayout,
         HermitianLayout, SymmetricLayout, TriangularLayout,
         UnknownLayout, AbstractBandedLayout, ApplyBroadcastStyle, ConjLayout, AbstractFillLayout, DualLayout,
-        colsupport, rowsupport, layout_getindex, QLayout, LayoutArray, LayoutMatrix, LayoutVector,
+        colsupport, rowsupport, layout_getindex, AbstractQLayout, LayoutArray, LayoutMatrix, LayoutVector,
         RangeCumsum
 
 if VERSION < v"1.7-"
@@ -81,7 +63,8 @@ abstract type LayoutArray{T,N} <: AbstractArray{T,N} end
 const LayoutMatrix{T} = LayoutArray{T,2}
 const LayoutVector{T} = LayoutArray{T,1}
 
-## TODO: Following are type piracy whch may be removed in Julia v1.5
+## TODO: Following are type piracy which may be removed in Julia v1.5
+## No, it can't, because strides(::AdjointAbsMat) is defined only for real eltype!
 _transpose_strides(a) = (a,1)
 _transpose_strides(a,b) = (b,a)
 strides(A::Adjoint) = _transpose_strides(strides(parent(A))...)
@@ -145,16 +128,16 @@ Base.@propagate_inbounds layout_getindex(A::AbstractArray, I::CartesianIndex) =
 
 macro _layoutgetindex(Typ)
     esc(quote
-        @inline Base.getindex(A::$Typ, kr::Colon, jr::Colon) = ArrayLayouts.layout_getindex(A, kr, jr)
-        @inline Base.getindex(A::$Typ, kr::Colon, jr::AbstractUnitRange) = ArrayLayouts.layout_getindex(A, kr, jr)
-        @inline Base.getindex(A::$Typ, kr::AbstractUnitRange, jr::Colon) = ArrayLayouts.layout_getindex(A, kr, jr)
-        @inline Base.getindex(A::$Typ, kr::AbstractUnitRange, jr::AbstractUnitRange) = ArrayLayouts.layout_getindex(A, kr, jr)
-        @inline Base.getindex(A::$Typ, kr::AbstractVector, jr::AbstractVector) = ArrayLayouts.layout_getindex(A, kr, jr)
-        @inline Base.getindex(A::$Typ, kr::Colon, jr::AbstractVector) = ArrayLayouts.layout_getindex(A, kr, jr)
-        @inline Base.getindex(A::$Typ, kr::Colon, jr::Integer) = ArrayLayouts.layout_getindex(A, kr, jr)
-        @inline Base.getindex(A::$Typ, kr::AbstractVector, jr::Colon) = ArrayLayouts.layout_getindex(A, kr, jr)
-        @inline Base.getindex(A::$Typ, kr::Integer, jr::Colon) = ArrayLayouts.layout_getindex(A, kr, jr)
-        @inline Base.getindex(A::$Typ, kr::Integer, jr::AbstractVector) = ArrayLayouts.layout_getindex(A, kr, jr)
+        @inline getindex(A::$Typ, kr::Colon, jr::Colon) = ArrayLayouts.layout_getindex(A, kr, jr)
+        @inline getindex(A::$Typ, kr::Colon, jr::AbstractUnitRange) = ArrayLayouts.layout_getindex(A, kr, jr)
+        @inline getindex(A::$Typ, kr::AbstractUnitRange, jr::Colon) = ArrayLayouts.layout_getindex(A, kr, jr)
+        @inline getindex(A::$Typ, kr::AbstractUnitRange, jr::AbstractUnitRange) = ArrayLayouts.layout_getindex(A, kr, jr)
+        @inline getindex(A::$Typ, kr::AbstractVector, jr::AbstractVector) = ArrayLayouts.layout_getindex(A, kr, jr)
+        @inline getindex(A::$Typ, kr::Colon, jr::AbstractVector) = ArrayLayouts.layout_getindex(A, kr, jr)
+        @inline getindex(A::$Typ, kr::Colon, jr::Integer) = ArrayLayouts.layout_getindex(A, kr, jr)
+        @inline getindex(A::$Typ, kr::AbstractVector, jr::Colon) = ArrayLayouts.layout_getindex(A, kr, jr)
+        @inline getindex(A::$Typ, kr::Integer, jr::Colon) = ArrayLayouts.layout_getindex(A, kr, jr)
+        @inline getindex(A::$Typ, kr::Integer, jr::AbstractVector) = ArrayLayouts.layout_getindex(A, kr, jr)
     end)
 end
 
@@ -243,7 +226,7 @@ end
 
 
 _copyto!(_, _, dest::AbstractArray{T,N}, src::AbstractArray{V,N}) where {T,V,N} =
-    Base.invoke(copyto!, Tuple{AbstractArray{T,N},AbstractArray{V,N}}, dest, src)
+    invoke(copyto!, Tuple{AbstractArray{T,N},AbstractArray{V,N}}, dest, src)
 
 
 _copyto!(dest, src) = _copyto!(MemoryLayout(dest), MemoryLayout(src), dest, src)
@@ -281,7 +264,7 @@ function zero!(_, A::AbstractArray{<:AbstractArray})
     A
 end
 
-_norm(_, A, p) = Base.invoke(norm, Tuple{Any,Real}, A, p)
+_norm(_, A, p) = invoke(norm, Tuple{Any,Real}, A, p)
 LinearAlgebra.norm(A::LayoutArray, p::Real=2) = _norm(MemoryLayout(A), A, p)
 LinearAlgebra.norm(A::SubArray{<:Any,N,<:LayoutArray}, p::Real=2) where N = _norm(MemoryLayout(A), A, p)
 
@@ -353,7 +336,7 @@ layout_replace_in_print_matrix(_, A, i, j, s) =
     i in colsupport(A,j) ? s : Base.replace_with_centered_mark(s)
 
 Base.replace_in_print_matrix(A::Union{LayoutVector,
-                                        LayoutMatrix,
+                                      LayoutMatrix,
                                       UpperTriangular{<:Any,<:LayoutMatrix},
                                       UnitUpperTriangular{<:Any,<:LayoutMatrix},
                                       LowerTriangular{<:Any,<:LayoutMatrix},

src/cumsum.jl

@@ -31,7 +31,7 @@ last(r::RangeCumsum) = sum(r.range)
 diff(r::RangeCumsum) = r.range[2:end]
 isempty(r::RangeCumsum) = isempty(r.range)
 
-union(a::RangeCumsum{<:Any,<:Base.OneTo}, b::RangeCumsum{<:Any,<:Base.OneTo}) = 
-    RangeCumsum(Base.OneTo(max(last(a.range),last(b.range))))
+union(a::RangeCumsum{<:Any,<:OneTo}, b::RangeCumsum{<:Any,<:OneTo}) = 
+    RangeCumsum(OneTo(max(last(a.range), last(b.range))))
 
 sort!(a::RangeCumsum{<:Any,<:AbstractUnitRange}) = a

src/factorizations.jl

@@ -123,7 +123,7 @@ adjointlayout(::Type, ::QRPackedQLayout{SLAY,TLAY}) where {SLAY,TLAY} = AdjQRPac
 adjointlayout(::Type, ::QRCompactWYQLayout{SLAY,TLAY}) where {SLAY,TLAY} = AdjQRCompactWYQLayout{SLAY,TLAY}()
 
 copy(M::Lmul{<:AbstractQLayout}) = copyto!(similar(M), M)
-mulreduce(M::Mul{<:AbstractQLayout,<:AbstractQLayout}) = MulAdd(M)
+mulreduce(M::Mul{<:AbstractQLayout,<:AbstractQLayout}) = Lmul(M)
 mulreduce(M::Mul{<:AbstractQLayout}) = Lmul(M)
 mulreduce(M::Mul{<:Any,<:AbstractQLayout}) = Rmul(M)
 mulreduce(M::Mul{<:TriangularLayout,<:AbstractQLayout}) = Rmul(M)
@@ -151,17 +151,17 @@ materialize!(M::Rmul{LAY}) where LAY<:AbstractQLayout = error("Overload material
 
 materialize!(M::Ldiv{<:AbstractQLayout}) = materialize!(Lmul(M.A',M.B))
 
-materialize!(M::Lmul{<:QRPackedQLayout{<:AbstractColumnMajor,<:AbstractColumnMajor},<:AbstractColumnMajor,<:AbstractMatrix{T},<:AbstractVecOrMat{T}}) where T<:BlasFloat =
+materialize!(M::Lmul{<:QRPackedQLayout{<:AbstractColumnMajor,<:AbstractColumnMajor},<:AbstractColumnMajor,<:AbstractQ{T},<:AbstractVecOrMat{T}}) where T<:BlasFloat =
     LAPACK.ormqr!('L','N',M.A.factors,M.A.τ,M.B)
 
-materialize!(M::Lmul{<:QRCompactWYQLayout{<:AbstractColumnMajor,<:AbstractColumnMajor},<:AbstractColumnMajor,<:AbstractMatrix{T},<:AbstractVecOrMat{T}}) where T<:BlasFloat =
+materialize!(M::Lmul{<:QRCompactWYQLayout{<:AbstractColumnMajor,<:AbstractColumnMajor},<:AbstractColumnMajor,<:AbstractQ{T},<:AbstractVecOrMat{T}}) where T<:BlasFloat =
     LAPACK.gemqrt!('L','N',M.A.factors,M.A.T,M.B)
 
-function materialize!(M::Lmul{<:QRCompactWYQLayout{<:AbstractColumnMajor,<:AbstractColumnMajor},<:AbstractRowMajor,<:AbstractMatrix{T},<:AbstractMatrix{T}}) where T<:BlasReal
+function materialize!(M::Lmul{<:QRCompactWYQLayout{<:AbstractColumnMajor,<:AbstractColumnMajor},<:AbstractRowMajor,<:AbstractQ{T},<:AbstractMatrix{T}}) where T<:BlasReal
     LAPACK.gemqrt!('R','T',M.A.factors,M.A.T,transpose(M.B))
     M.B
 end
-function materialize!(M::Lmul{<:QRCompactWYQLayout{<:AbstractColumnMajor,<:AbstractColumnMajor},<:ConjLayout{<:AbstractRowMajor},<:AbstractMatrix{T},<:AbstractMatrix{T}}) where T<:BlasComplex
+function materialize!(M::Lmul{<:QRCompactWYQLayout{<:AbstractColumnMajor,<:AbstractColumnMajor},<:ConjLayout{<:AbstractRowMajor},<:AbstractQ{T},<:AbstractMatrix{T}}) where T<:BlasComplex
     LAPACK.gemqrt!('R','C',M.A.factors,M.A.T,(M.B)')
     M.B
 end
@@ -207,7 +207,7 @@ materialize!(M::Lmul{<:AdjQRCompactWYQLayout{<:AbstractStridedLayout,<:AbstractS
 function materialize!(M::Lmul{<:AdjQRPackedQLayout})
     adjA,B = M.A, M.B
     require_one_based_indexing(B)
-    A = adjA.parent
+    A = parent(adjA)
     mA, nA = size(A.factors)
     mB, nB = size(B,1), size(B,2)
     if mA != mB
@@ -233,9 +233,9 @@ function materialize!(M::Lmul{<:AdjQRPackedQLayout})
 end
 
 ## AQ
-materialize!(M::Rmul{<:AbstractStridedLayout,<:QRPackedQLayout{<:AbstractStridedLayout,<:AbstractStridedLayout},<:AbstractVecOrMat{T},<:AbstractMatrix{T}}) where T<:BlasFloat =
+materialize!(M::Rmul{<:AbstractStridedLayout,<:QRPackedQLayout{<:AbstractStridedLayout,<:AbstractStridedLayout},<:AbstractVecOrMat{T},<:AbstractQ{T}}) where T<:BlasFloat =
     LAPACK.ormqr!('R', 'N', M.B.factors, M.B.τ, M.A)
-materialize!(M::Rmul{<:AbstractStridedLayout,<:QRCompactWYQLayout{<:AbstractStridedLayout,<:AbstractStridedLayout},<:AbstractVecOrMat{T},<:AbstractMatrix{T}}) where T<:BlasFloat =
+materialize!(M::Rmul{<:AbstractStridedLayout,<:QRCompactWYQLayout{<:AbstractStridedLayout,<:AbstractStridedLayout},<:AbstractVecOrMat{T},<:AbstractQ{T}}) where T<:BlasFloat =
     LAPACK.gemqrt!('R','N', M.B.factors, M.B.T, M.A)
 function materialize!(M::Rmul{<:Any,<:QRPackedQLayout})
     A,Q = M.A,M.B
@@ -274,7 +274,7 @@ materialize!(M::Rmul{<:AbstractStridedLayout,<:AdjQRCompactWYQLayout{<:AbstractS
     (B = M.B.parent; LAPACK.gemqrt!('R','C',B.factors,B.T,M.A))
 function materialize!(M::Rmul{<:Any,<:AdjQRPackedQLayout})
     A,adjQ = M.A,M.B
-    Q = adjQ.parent
+    Q = parent(adjQ)
     mQ, nQ = size(Q.factors)
     mA, nA = size(A,1), size(A,2)
     if nA != mQ
@@ -300,12 +300,12 @@ function materialize!(M::Rmul{<:Any,<:AdjQRPackedQLayout})
 end
 
 
-__qr(layout, lengths, A; kwds...) = Base.invoke(qr, Tuple{AbstractMatrix{eltype(A)}}, A; kwds...)
+__qr(layout, lengths, A; kwds...) = 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...)
+_qr(layout, axes, A, pivot::P; kwds...) where P = invoke(qr, Tuple{AbstractMatrix{eltype(A)},P}, A, pivot; kwds...)
 _qr!(layout, axes, A, args...; kwds...) = error("Overload _qr!(::$(typeof(layout)), axes, A)")
-_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; kwds...) = invoke(lu, Tuple{AbstractMatrix{eltype(A)}}, A; kwds...)
+_lu(layout, axes, A, pivot::P; kwds...) where P = 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, ::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)
@@ -384,7 +384,7 @@ function _cholesky!(::SymmetricLayout{<:AbstractColumnMajor}, axes, A::AbstractM
 end
 
 
-_inv_eye(_, ::Type{T}, axs::NTuple{2,Base.OneTo{Int}}) where T = Matrix{T}(I, map(length,axs)...)
+_inv_eye(_, ::Type{T}, axs::NTuple{2,OneTo{Int}}) where T = Matrix{T}(I, map(length,axs)...)
 function _inv_eye(A, ::Type{T}, (rows,cols)) where T
     dest = zero!(similar(A, T, (cols,rows)))
     dest[diagind(dest)] .= one(T)
@@ -414,7 +414,7 @@ macro _layoutfactorizations(Typ)
         LinearAlgebra.lu(A::$Typ{T}; kwds...) where T = ArrayLayouts._lu(ArrayLayouts.MemoryLayout(A), axes(A), A; kwds...)
         LinearAlgebra.lu!(A::$Typ, args...; kwds...) = ArrayLayouts._lu!(ArrayLayouts.MemoryLayout(A), axes(A), A, args...; kwds...)
         LinearAlgebra.factorize(A::$Typ) = ArrayLayouts._factorize(ArrayLayouts.MemoryLayout(A), axes(A), A)
-        LinearAlgebra.inv(A::$Typ) = ArrayLayouts._inv(ArrayLayouts.MemoryLayout(A), axes(A), A)
+        Base.inv(A::$Typ) = ArrayLayouts._inv(ArrayLayouts.MemoryLayout(A), axes(A), A)
         LinearAlgebra.ldiv!(L::LU{<:Any,<:$Typ}, B) = ArrayLayouts.ldiv!(L, B)
         LinearAlgebra.ldiv!(L::LU{<:Any,<:$Typ}, B::$Typ) = ArrayLayouts.ldiv!(L, B)
     end)