Julia v1.6+ overhaul (#163)

Author: dkarrasch

.github/workflows/ci.yml

@@ -14,7 +14,7 @@ jobs:
     strategy:
       matrix:
         version:
-          - '1.0'
+          - '1.6'
           - '1'
           - 'nightly'
         os:

Project.toml

@@ -7,7 +7,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [compat]
-julia = "1"
+julia = "1.6"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"

src/LinearMaps.jl

@@ -8,12 +8,7 @@ using LinearAlgebra
 import LinearAlgebra: mul!
 using SparseArrays
 
-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
+using Base: require_one_based_indexing
 
 abstract type LinearMap{T} end
 
@@ -33,11 +28,7 @@ MulStyle(::ThreeArg, ::FiveArg) = ThreeArg()
 MulStyle(::FiveArg, ::ThreeArg) = ThreeArg()
 MulStyle(::ThreeArg, ::ThreeArg) = ThreeArg()
 MulStyle(::LinearMap) = ThreeArg() # default
-@static if VERSION ≥ v"1.3.0-alpha.115"
-    MulStyle(::AbstractVecOrMat) = FiveArg()
-else
-    MulStyle(::AbstractVecOrMat) = ThreeArg()
-end
+MulStyle(::AbstractVecOrMat) = FiveArg()
 MulStyle(A::LinearMap, As::LinearMap...) = MulStyle(MulStyle(A), MulStyle(As...))
 
 Base.isreal(A::LinearMap) = eltype(A) <: Real
@@ -117,9 +108,8 @@ function Base.:(*)(A::LinearMap, x::AbstractVector)
     y = similar(x, T, axes(A)[1])
     return mul!(y, A, x)
 end
-if VERSION ≥ v"1.3"
-    (L::LinearMap)(x::AbstractVector) = L*x
-end
+
+(L::LinearMap)(x::AbstractVector) = L*x
 
 """
     mul!(Y::AbstractVecOrMat, A::LinearMap, B::AbstractVector) -> Y
@@ -225,13 +215,9 @@ function mul!(Y::AbstractMatrix, A::LinearMap, X::AbstractMatrix, α::Number, β
 end
 
 function _generic_mapmat_mul!(Y, A, X, α=true, β=false)
-    @views for i in 1:size(X, 2)
-        _unsafe_mul!(Y[:, i], A, X[:, i], α, β)
+    for (Xi, Yi) in zip(eachcol(X), eachcol(Y))
+        mul!(Yi, A, Xi, α, β)
     end
-    # starting from Julia v1.1, we could use the `eachcol` iterator
-    # for (Xi, Yi) in zip(eachcol(X), eachcol(Y))
-    #     mul!(Yi, A, Xi, α, β)
-    # end
     return Y
 end
 

src/blockmap.jl

@@ -24,6 +24,13 @@ BlockMap{T}(maps::As, rows::Rs) where {T, As<:LinearMapTuple, Rs} =
 
 MulStyle(A::BlockMap) = MulStyle(A.maps...)
 
+function _getranges(maps, dim, inds=ntuple(identity, Val(length(maps))))
+    sizes = map(i -> size(maps[i], dim)::Int, inds)
+    ends = cumsum(sizes)
+    starts = (1, (1 .+ Base.front(ends))...)
+    return UnitRange.(starts, ends)
+end
+
 """
     rowcolranges(maps, rows)
 
@@ -32,24 +39,20 @@ map in `maps`, according to its position in a virtual matrix representation of t
 block linear map obtained from `hvcat(rows, maps...)`.
 """
 function rowcolranges(maps, rows)
-    rowranges = ntuple(n->1:0, Val(length(rows)))
-    colranges = ntuple(n->1:0, Val(length(maps)))
-    mapind = 0
-    rowstart = 1
-    for (i, row) in enumerate(rows)
-        mapind += 1
-        rowend = rowstart + Int(size(maps[mapind], 1))::Int - 1
-        rowranges = Base.setindex(rowranges, rowstart:rowend, i)
-        colstart = 1
-        colend = Int(size(maps[mapind], 2))::Int
-        colranges = Base.setindex(colranges, colstart:colend, mapind)
-        for colind in 2:row
-            mapind += 1
-            colstart = colend + 1
-            colend += Int(size(maps[mapind], 2))::Int
-            colranges = Base.setindex(colranges, colstart:colend, mapind)
-        end
-        rowstart = rowend + 1
+    # find indices of the row-wise first maps
+    firstmapinds = cumsum((1, Base.front(rows)...))
+    # compute rowranges from first dimension of the row-wise first maps
+    rowranges = _getranges(maps, 1, firstmapinds)
+
+    # compute ranges from second dimension as if all in one row
+    temp = _getranges(maps, 2)
+    # introduce "line breaks"
+    colranges = ntuple(Val(length(maps))) do i
+        # for each map find the index of the respective row-wise first map
+        # something-trick just to assure the compiler that the index is an Int
+        @inbounds firstmapind = firstmapinds[something(findlast(<=(i), firstmapinds), 1)]
+        # shift ranges by the first col-index of the row-wise first map
+        return @inbounds temp[i] .- first(temp[firstmapind]) .+ 1
     end
     return rowranges, colranges
 end
@@ -82,17 +85,13 @@ function Base.hcat(As::Union{LinearMap, UniformScaling, AbstractVecOrMat}...)
     T = promote_type(map(eltype, As)...)
     nbc = length(As)
 
-    nrows = -1
     # find first non-UniformScaling to detect number of rows
-    for A in As
-        if !(A isa UniformScaling)
-            nrows = size(A, 1)
-            break
-        end
-    end
-    @assert nrows != -1
+    j = findfirst(A -> !isa(A, UniformScaling), As)
     # this should not happen, function should only be called with at least one LinearMap
-    return BlockMap{T}(promote_to_lmaps(ntuple(i->nrows, nbc), 1, 1, As...), (nbc,))
+    @assert !isnothing(j)
+    @inbounds nrows = size(As[j], 1)::Int
+    
+    return BlockMap{T}(promote_to_lmaps(ntuple(_ -> nrows, Val(nbc)), 1, 1, As...), (nbc,))
 end
 
 ############
@@ -124,18 +123,14 @@ function Base.vcat(As::Union{LinearMap,UniformScaling,AbstractVecOrMat}...)
     T = promote_type(map(eltype, As)...)
     nbr = length(As)
 
-    ncols = -1
-    # find first non-UniformScaling to detect number of columns
-    for A in As
-        if !(A isa UniformScaling)
-            ncols = size(A, 2)
-            break
-        end
-    end
-    @assert ncols != -1
+    # find first non-UniformScaling to detect number of rows
+    j = findfirst(A -> !isa(A, UniformScaling), As)
     # this should not happen, function should only be called with at least one LinearMap
-    rows = ntuple(i->1, nbr)
-    return BlockMap{T}(promote_to_lmaps(ntuple(i->ncols, nbr), 1, 2, As...), rows)
+    @assert !isnothing(j)
+    @inbounds ncols = size(As[j], 2)::Int
+
+    rows = ntuple(_ -> 1, Val(nbr))
+    return BlockMap{T}(promote_to_lmaps(ntuple(_ -> ncols, Val(nbr)), 1, 2, As...), rows)
 end
 
 ############
@@ -181,7 +176,7 @@ function Base.hvcat(rows::Tuple{Vararg{Int}},
         ni = -1 # number of rows in this block-row, -1 indicates unknown
         for k in 1:rows[i]
             if !isa(As[j+k], UniformScaling)
-                na = size(As[j+k], 1)
+                na = size(As[j+k], 1)::Int
                 ni >= 0 && ni != na &&
                     throw(DimensionMismatch("mismatch in number of rows"))
                 ni = na
@@ -199,7 +194,7 @@ function Base.hvcat(rows::Tuple{Vararg{Int}},
         nci = 0
         rows[i] > 0 && n[j+1] == -1 && (j += rows[i]; continue)
         for k = 1:rows[i]
-            nci += isa(As[j+k], UniformScaling) ? n[j+k] : size(As[j+k], 2)
+            nci += isa(As[j+k], UniformScaling) ? n[j+k] : size(As[j+k], 2)::Int
         end
         nc >= 0 && nc != nci && throw(DimensionMismatch("mismatch in number of columns"))
         nc = nci
@@ -412,14 +407,8 @@ struct BlockDiagonalMap{T,
             promote_type(T, TA) == T ||
                 error("eltype $TA cannot be promoted to $T in BlockDiagonalMap constructor")
         end
-        # row ranges
-        inds = vcat(1, size.(maps, 1)...)
-        cumsum!(inds, inds)
-        rowranges = ntuple(i -> inds[i]:inds[i+1]-1, Val(length(maps)))
-        # column ranges
-        inds[2:end] .= size.(maps, 2)
-        cumsum!(inds, inds)
-        colranges = ntuple(i -> inds[i]:inds[i+1]-1, Val(length(maps)))
+        rowranges = _getranges(maps, 1)
+        colranges = _getranges(maps, 2)
         return new{T, As, typeof(rowranges)}(maps, rowranges, colranges)
     end
 end
@@ -476,7 +465,7 @@ object among the first 8 arguments.
 """
 Base.cat
 
-Base.size(A::BlockDiagonalMap) = (last(A.rowranges[end]), last(A.colranges[end]))
+Base.size(A::BlockDiagonalMap) = (last(last(A.rowranges)), last(last(A.colranges)))
 
 MulStyle(A::BlockDiagonalMap) = MulStyle(A.maps...)
 

src/composition.jl

@@ -5,8 +5,7 @@ struct CompositeMap{T, As<:LinearMapTuple} <: LinearMap{T}
         for n in 2:N
             check_dim_mul(maps[n], maps[n-1])
         end
-        for TA in Base.Generator(eltype, maps)
-            # like lazy map; could use Base.Iterators.map in Julia >= 1.6
+        for TA in Base.Iterators.map(eltype, maps)
             promote_type(T, TA) == T ||
                 error("eltype $TA cannot be promoted to $T in CompositeMap constructor")
         end

src/kronecker.jl

@@ -1,7 +1,7 @@
 struct KroneckerMap{T, As<:LinearMapTuple} <: LinearMap{T}
     maps::As
     function KroneckerMap{T}(maps::LinearMapTuple) where {T}
-        for TA in Base.Generator(eltype, maps)
+        for TA in Base.Iterators.map(eltype, maps)
             promote_type(T, TA) == T ||
                 error("eltype $TA cannot be promoted to $T in KroneckerMap constructor")
         end
@@ -221,7 +221,7 @@ struct KroneckerSumMap{T, As<:Tuple{LinearMap, LinearMap}} <: LinearMap{T}
         A1, A2 = maps
         (size(A1, 1) == size(A1, 2) && size(A2, 1) == size(A2, 2)) ||
             throw(ArgumentError("operators need to be square in Kronecker sums"))
-        for TA in Base.Generator(eltype, maps)
+        for TA in Base.Iterators.map(eltype, maps)
             promote_type(T, TA) == T ||
                 error("eltype $TA cannot be promoted to $T in KroneckerSumMap constructor")
         end

src/linearcombination.jl

@@ -18,7 +18,7 @@ MulStyle(A::LinearCombination) = MulStyle(A.maps...)
 
 # basic methods
 Base.size(A::LinearCombination) = size(A.maps[1])
-Base.axes(A::LinearMaps.LinearCombination) = axes(A.maps[1])
+Base.axes(A::LinearCombination) = axes(A.maps[1])
 # following conditions are sufficient but not necessary
 LinearAlgebra.issymmetric(A::LinearCombination) = all(issymmetric, A.maps)
 LinearAlgebra.ishermitian(A::LinearCombination) = all(ishermitian, A.maps)

src/wrappedmap.jl

@@ -87,40 +87,38 @@ mul!(Y::AbstractMatrix, X::AbstractMatrix, A::VecOrMatMap) = mul!(Y, X, A.lmap)
 # the following method is needed for disambiguation with left-multiplication
 mul!(Y::AbstractMatrix, X::TransposeAbsVecOrMat, A::VecOrMatMap) = mul!(Y, X, A.lmap)
 
-if VERSION ≥ v"1.3.0-alpha.115"
-    for (In, Out) in ((AbstractVector, AbstractVecOrMat), (AbstractMatrix, AbstractMatrix))
-        @eval begin
-            function _unsafe_mul!(y::$Out, A::WrappedMap, x::$In, α::Number, β::Number)
-                return _unsafe_mul!(y, A.lmap, x, α, β)
-            end
-            function _unsafe_mul!(y::$Out, At::TransposeMap{<:Any,<:WrappedMap}, x::$In,
-                                    α::Number, β::Number)
-                A = At.lmap
-                return (issymmetric(A) || (isreal(A) && ishermitian(A))) ?
-                    _unsafe_mul!(y, A.lmap, x, α, β) :
-                    _unsafe_mul!(y, transpose(A.lmap), x, α, β)
-            end
-            function _unsafe_mul!(y::$Out, Ac::AdjointMap{<:Any,<:WrappedMap}, x::$In, α::Number, β::Number)
-                A = Ac.lmap
-                return ishermitian(A) ?
-                    _unsafe_mul!(y, A.lmap, x, α, β) :
-                    _unsafe_mul!(y, adjoint(A.lmap), x, α, β)
-            end
+for (In, Out) in ((AbstractVector, AbstractVecOrMat), (AbstractMatrix, AbstractMatrix))
+    @eval begin
+        function _unsafe_mul!(y::$Out, A::WrappedMap, x::$In, α::Number, β::Number)
+            return _unsafe_mul!(y, A.lmap, x, α, β)
+        end
+        function _unsafe_mul!(y::$Out, At::TransposeMap{<:Any,<:WrappedMap}, x::$In,
+                                α::Number, β::Number)
+            A = At.lmap
+            return (issymmetric(A) || (isreal(A) && ishermitian(A))) ?
+                _unsafe_mul!(y, A.lmap, x, α, β) :
+                _unsafe_mul!(y, transpose(A.lmap), x, α, β)
+        end
+        function _unsafe_mul!(y::$Out, Ac::AdjointMap{<:Any,<:WrappedMap}, x::$In, α::Number, β::Number)
+            A = Ac.lmap
+            return ishermitian(A) ?
+                _unsafe_mul!(y, A.lmap, x, α, β) :
+                _unsafe_mul!(y, adjoint(A.lmap), x, α, β)
         end
     end
+end
 
-    mul!(X::AbstractMatrix, Y::AbstractMatrix, A::VecOrMatMap, α::Number, β::Number) =
-        mul!(X, Y, A.lmap, α, β)
-    # the following method is needed for disambiguation with left-multiplication
-    function mul!(Y::AbstractMatrix{<:RealOrComplex}, X::AbstractMatrix{<:RealOrComplex}, A::VecOrMatMap{<:RealOrComplex},
-                    α::RealOrComplex, β::RealOrComplex)
-        return mul!(Y, X, A.lmap, α, β)
-    end
-    function mul!(Y::AbstractMatrix{<:RealOrComplex}, X::TransposeAbsVecOrMat{<:RealOrComplex}, A::VecOrMatMap{<:RealOrComplex},
-                    α::RealOrComplex, β::RealOrComplex)
-        return mul!(Y, X, A.lmap, α, β)
-    end
-end # VERSION
+mul!(X::AbstractMatrix, Y::AbstractMatrix, A::VecOrMatMap, α::Number, β::Number) =
+    mul!(X, Y, A.lmap, α, β)
+# the following 2 methods are needed for disambiguation with left-multiplication
+function mul!(Y::AbstractMatrix{<:RealOrComplex}, X::AbstractMatrix{<:RealOrComplex}, A::VecOrMatMap{<:RealOrComplex},
+                α::RealOrComplex, β::RealOrComplex)
+    return mul!(Y, X, A.lmap, α, β)
+end
+function mul!(Y::AbstractMatrix{<:RealOrComplex}, X::TransposeAbsVecOrMat{<:RealOrComplex}, A::VecOrMatMap{<:RealOrComplex},
+                α::RealOrComplex, β::RealOrComplex)
+    return mul!(Y, X, A.lmap, α, β)
+end
 
 # combine LinearMap and Matrix objects: linear combinations and map composition
 Base.:(+)(A₁::LinearMap, A₂::AbstractMatrix) = +(A₁, WrappedMap(A₂))

test/blockmap.jl

@@ -1,4 +1,5 @@
 using Test, LinearMaps, LinearAlgebra, SparseArrays, BenchmarkTools, InteractiveUtils
+using LinearMaps: FiveArg, ThreeArg
 
 @testset "block maps" begin
     @testset "hcat" begin
@@ -8,7 +9,7 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays, BenchmarkTools, Interactive
             v = rand(elty, 10)
             L = @inferred hcat(LinearMap(A11), LinearMap(A12))
             @test occursin("10×$(10+n2) LinearMaps.BlockMap{$elty}", sprint((t, s) -> show(t, "text/plain", s), L))
-            @test @inferred(LinearMaps.MulStyle(L)) === matrixstyle
+            @test @inferred(LinearMaps.MulStyle(L)) === FiveArg()
             @test L isa LinearMaps.BlockMap{elty}
             if elty <: Complex
                 @test_throws ErrorException LinearMaps.BlockMap{Float64}((LinearMap(A11), LinearMap(A12)), (2,))
@@ -54,7 +55,7 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays, BenchmarkTools, Interactive
             L = @inferred vcat(LinearMap(A11), LinearMap(A21))
             @test occursin("30×10 LinearMaps.BlockMap{$elty}", sprint((t, s) -> show(t, "text/plain", s), L))
             @test L isa LinearMaps.BlockMap{elty}
-            @test @inferred(LinearMaps.MulStyle(L)) === matrixstyle
+            @test @inferred(LinearMaps.MulStyle(L)) === FiveArg()
             @test (@which [A11; A21]).module != LinearMaps
             A = [A11; A21]
             x = rand(10)
@@ -84,7 +85,7 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays, BenchmarkTools, Interactive
             @test (@which [A11 A12; A21 A22]).module != LinearMaps
             @inferred hvcat((2,2), LinearMap(A11), LinearMap(A12), LinearMap(A21), LinearMap(A22))
             L = [LinearMap(A11) LinearMap(A12); LinearMap(A21) LinearMap(A22)]
-            @test @inferred(LinearMaps.MulStyle(L)) === matrixstyle
+            @test @inferred(LinearMaps.MulStyle(L)) === FiveArg()
             @test @inferred !issymmetric(L)
             @test @inferred !ishermitian(L)
             x = rand(30)
@@ -142,7 +143,7 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays, BenchmarkTools, Interactive
             A12 = rand(elty, 10, 10)
             A = [I A12; transform(A12) I]
             L = [I LinearMap(A12); transform(LinearMap(A12)) I]
-            @test @inferred(LinearMaps.MulStyle(L)) === matrixstyle
+            @test @inferred(LinearMaps.MulStyle(L)) === FiveArg()
             if elty <: Complex
                 if transform == transpose
                     @test @inferred issymmetric(L)
@@ -201,7 +202,7 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays, BenchmarkTools, Interactive
             @test (@which cat(M1, M2, M3, M2, M1; dims=(1,2))).module != LinearMaps
             x = randn(elty, size(Md, 2))
             Bd = @inferred blockdiag(L1, L2, L3, L2, L1)
-            @test @inferred(LinearMaps.MulStyle(Bd)) === matrixstyle
+            @test @inferred(LinearMaps.MulStyle(Bd)) === FiveArg()
             @test occursin("25×39 LinearMaps.BlockDiagonalMap{$elty}", sprint((t, s) -> show(t, "text/plain", s), Bd))
             @test Matrix(Bd) == Md
             @test convert(AbstractMatrix, Bd) isa SparseMatrixCSC