Minor polish before v3.7 release (#182)

Author: dkarrasch

README.md

@@ -39,5 +39,5 @@ julia> import Pkg; Pkg.add("LinearMaps")
 [license-img]: http://img.shields.io/badge/license-MIT-brightgreen.svg?style=flat
 [license-url]: LICENSE.md
 
-[aqua-img]: https://img.shields.io/badge/Aqua.jl-%F0%9F%8C%A2-aqua.svg
+[aqua-img]: https://raw.githubusercontent.com/JuliaTesting/Aqua.jl/master/badge.svg
 [aqua-url]: https://github.com/JuliaTesting/Aqua.jl

docs/src/custom.jl

@@ -73,7 +73,7 @@ mul!(ones(3,3), A, reshape(collect(1:9), 3, 3), 2, 2)
 using BenchmarkTools
 
 @benchmark mul!($(zeros(3)), $A, $x)
- 
+
 #-
 
 @benchmark mul!($(zeros(3)), $A, $x, $(rand()), $(rand()))
@@ -150,12 +150,13 @@ transpose(A)*x
 
 try MyFillMap(5.0, (3, 4))' * ones(3) catch e println(e) end
 
-# which require explicit adjoint/transpose handling, for which there exist two *distinct* paths.
+# which require explicit adjoint/transpose handling, for which there exist two *distinct*
+# paths.
 
 # ### Path 1: Generic, non-invariant `LinearMap` subtypes
 
-# The first option is to write `LinearAlgebra.mul!` methods for the corresponding wrapped
-# map types; for instance,
+# The first option is to write `LinearMaps._unsafe_mul!` methods for the corresponding
+# wrapped map types; for instance,
 
 function LinearMaps._unsafe_mul!(
     y::AbstractVecOrMat,
@@ -228,3 +229,56 @@ mul!(similar(x)', x', A)
 # corresponding 5-arg `mul!` methods. This may seem like a lot of methods to
 # be implemented, but note that adding such methods is only necessary/recommended
 # for performance.
+
+# ## Computing a matrix representation
+
+# In some cases, it might be necessary to compute a matrix representation of a `LinearMap`.
+# This is essentially done via the
+# `[LinearMaps._unsafe_mul!(::Matrix,::LinearMap,::Number)]`(@ref) method, for which a
+# generic fallback exists: it applies the `LinearMap` successively to the standard unit
+# vectors.
+
+F = MyFillMap(5, (100,100))
+M = Matrix{eltype(F)}(undef, size(F))
+@benchmark Matrix($F)
+
+#-
+
+@benchmark LinearMaps._unsafe_mul!($(Matrix{Int}(undef, (100,100))), $(MyFillMap(5, (100,100))), true)
+
+# If a more performant implementation exists, it is recommended to overwrite this method,
+# for instance (as before, size checks need not be included here since they are handled by
+# the corresponding `LinearAlgebra.mul!` method):
+
+LinearMaps._unsafe_mul!(M::AbstractMatrix, A::MyFillMap, s::Number) = fill!(M, A.λ*s)
+@benchmark Matrix($F)
+
+#-
+
+@benchmark LinearMaps._unsafe_mul!($(Matrix{Int}(undef, (100,100))), $(MyFillMap(5, (100,100))), true)
+
+# As one can see, the above runtimes are dominated by the allocation of the output matrix,
+# but still overwriting the multiplication kernel yields a speed-up of about factor 3 for
+# the matrix filling part.
+
+# ## Slicing
+
+# As usual, generic fallbacks for `LinearMap` slicing exist and are handled by the following
+# method hierarchy, where at least one of `I` and `J` has to be a `Colon`:
+#
+#     Base.getindex(::LinearMap, I, J)
+#     -> LinearMaps._getindex(::LinearMap, I, J)
+#
+# The method `Base.getindex` checks the validity of the the requested indices and calls
+# `LinearMaps._getindex`, which should be overloaded for custom `LinearMap`s subtypes.
+# For instance:
+
+@benchmark F[1,:]
+
+#-
+
+LinearMaps._getindex(A::MyFillMap, ::Integer, J::Base.Slice) = fill(A.λ, axes(J))
+@benchmark F[1,:]
+
+# Note that in `Base.getindex` `Colon`s are converted to `Base.Slice` via
+# `Base.to_indices`, thus the dispatch must be on `Base.Slice` rather than on `Colon`.

docs/src/index.md

@@ -50,17 +50,17 @@ import Arpack, IterativeSolvers, KrylovKit, TSVD, ArnoldiMethod
 # Example 1, 1-dimensional Laplacian with periodic boundary conditions
 function leftdiff!(y::AbstractVector, x::AbstractVector) # left difference assuming periodic boundary conditions
     N = length(x)
-    length(y) == N || throw(DimensionMismatch())
-    @inbounds for i=1:N
+    axes(y) == axes(x) || throw(DimensionMismatch())
+    @inbounds for i in eachindex(x, y)
         y[i] = x[i] - x[mod1(i-1, N)]
     end
     return y
 end
 
 function mrightdiff!(y::AbstractVector, x::AbstractVector) # minus right difference
     N = length(x)
-    length(y) == N || throw(DimensionMismatch())
-    @inbounds for i=1:N
+    axes(y) == axes(x) || throw(DimensionMismatch())
+    @inbounds for i in eachindex(x, y)
         y[i] = x[i] - x[mod1(i+1, N)]
     end
     return y
@@ -107,7 +107,7 @@ In Julia v1.3 and above, the last line can be simplified to
 
 ```julia
 KrylovKit.eigsolve(Δ, size(Δ, 1), 3, :LR)
-````
+```
 
 leveraging the fact that objects of type `L <: LinearMap` are callable.
 

docs/src/related.md

@@ -8,6 +8,7 @@ The following open-source packages provide similar or even extended functionalit
   [`LinearOperators.jl`](https://github.com/JuliaSmoothOptimizers/LinearOperators.jl)
   and the Python package [`PyLops`](https://github.com/equinor/pylops)
 * [`fastmat`: fast linear transforms in Python](https://pypi.org/project/fastmat/)
+* [`JOLI.jl`: Julia Operators LIbrary](https://github.com/slimgroup/JOLI.jl)
 * [`FunctionOperators.jl`](https://github.com/hakkelt/FunctionOperators.jl)
   and [`LinearMapsAA.jl`](https://github.com/JeffFessler/LinearMapsAA.jl)
   also support mappings between `Array`s, inspired by the `fatrix` object type in the
@@ -19,8 +20,10 @@ there exist further related packages in the Julia ecosystem:
 * [`LazyArrays.jl`](https://github.com/JuliaArrays/LazyArrays.jl)
 * [`BlockArrays.jl`](https://github.com/JuliaArrays/BlockArrays.jl)
 * [`BlockDiagonals.jl`](https://github.com/invenia/BlockDiagonals.jl)
+* [`BlockFactorizations.jl`](https://github.com/SebastianAment/BlockFactorizations.jl)
 * [`Kronecker.jl`](https://github.com/MichielStock/Kronecker.jl)
 * [`FillArrays.jl`](https://github.com/JuliaArrays/FillArrays.jl)
+* [`LiftedMaps.jl](https://github.com/krcools/LiftedMaps.jl)
 
 Since these packages provide types that are subtypes of Julia `Base`'s `AbstractMatrix` type,
 objects of those types can be wrapped by a `LinearMap` and freely mixed with, for instance,

docs/src/types.md

@@ -92,6 +92,10 @@ Type for lazily representing constantly filled matrices.
 LinearMaps.FillMap
 ```
 
+### `EmbeddedMap`
+
+Type for representing linear maps that are embedded in larger zero maps.
+
 ## Methods
 
 ### Multiplication methods
@@ -104,7 +108,7 @@ LinearAlgebra.mul!(::AbstractVecOrMat,::LinearMap,::AbstractVector)
 LinearAlgebra.mul!(::AbstractVecOrMat,::LinearMap,::AbstractVector,::Number,::Number)
 LinearAlgebra.mul!(::AbstractMatrix,::AbstractMatrix,::LinearMap)
 LinearAlgebra.mul!(::AbstractVecOrMat,::LinearMap,::Number)
-LinearAlgebra.mul!(::AbstractVecOrMat,::LinearMap,::Number,::Number,::Number)
+LinearAlgebra.mul!(::AbstractMatrix,::LinearMap,::Number,::Number,::Number)
 *(::LinearAlgebra.AdjointAbsVec,::LinearMap)
 *(::LinearAlgebra.TransposeAbsVec,::LinearMap)
 ```
@@ -139,3 +143,13 @@ as in the usual matrix case: `transpose(A) * x` and `mul!(y, A', x)`, for instan
   purposes or if you want to have the explicit sparse matrix representation of
   a linear map for which you only have a function definition (e.g. to be able
   to use its `transpose` or `adjoint`).
+
+### Slicing methods
+
+Complete slicing, i.e., `A[:,j]`, `A[:,J]`, `A[i,:]`, `A[I,:]` and `A[:,:]` for `i`, `j`
+`Integer` and `I`, `J` `AbstractVector{<:Integer}` is generically available for any
+`A::LinearMap` subtype via application of `A` (or `A'` for (predominantly) horizontal
+slicing) to standard unit vectors of appropriate length. By complete slicing we refer
+two-dimensional Cartesian indexing where at least one of the "indices" is a colon. This is
+facilitated by overloads of `Base.getindex`. Partial slicing à la `A[I,J]` and scalar or
+linear indexing are _not_ supported.

src/LinearMaps.jl

@@ -258,16 +258,16 @@ end
 _unsafe_mul!(y, A::MapOrVecOrMat, x) = mul!(y, A, x)
 _unsafe_mul!(y, A::AbstractVecOrMat, x, α, β) = mul!(y, A, x, α, β)
 _unsafe_mul!(y::AbstractVecOrMat, A::LinearMap, x::AbstractVector, α, β) =
-    _generic_mapvec_mul!(y, A, x, α, β)
+    _generic_map_mul!(y, A, x, α, β)
 _unsafe_mul!(y::AbstractMatrix, A::LinearMap, x::AbstractMatrix) =
-    _generic_mapmat_mul!(y, A, x)
+    _generic_map_mul!(y, A, x)
 _unsafe_mul!(y::AbstractMatrix, A::LinearMap, x::AbstractMatrix, α::Number, β::Number) =
-    _generic_mapmat_mul!(y, A, x, α, β)
-_unsafe_mul!(Y::AbstractMatrix, A::LinearMap, s::Number) = _generic_mapnum_mul!(Y, A, s)
+    _generic_map_mul!(y, A, x, α, β)
+_unsafe_mul!(Y::AbstractMatrix, A::LinearMap, s::Number) = _generic_map_mul!(Y, A, s)
 _unsafe_mul!(Y::AbstractMatrix, A::LinearMap, s::Number, α::Number, β::Number) =
-    _generic_mapnum_mul!(Y, A, s, α, β)
+    _generic_map_mul!(Y, A, s, α, β)
 
-function _generic_mapvec_mul!(y, A, x, α, β)
+function _generic_map_mul!(y, A, x::AbstractVector, α, β)
     # this function needs to call mul! for, e.g.,  AdjointMap{...,<:CustomMap}
     if isone(α)
         iszero(β) && return mul!(y, A, x)
@@ -294,21 +294,19 @@ function _generic_mapvec_mul!(y, A, x, α, β)
         return y
     end
 end
-
-function _generic_mapmat_mul!(Y, A, X)
+function _generic_map_mul!(Y, A, X::AbstractMatrix)
     for (Xi, Yi) in zip(eachcol(X), eachcol(Y))
         mul!(Yi, A, Xi)
     end
     return Y
 end
-function _generic_mapmat_mul!(Y, A, X, α, β)
+function _generic_map_mul!(Y, A, X::AbstractMatrix, α, β)
     for (Xi, Yi) in zip(eachcol(X), eachcol(Y))
         mul!(Yi, A, Xi, α, β)
     end
     return Y
 end
-
-function _generic_mapnum_mul!(Y, A, s)
+function _generic_map_mul!(Y, A, s::Number)
     T = promote_type(eltype(A), typeof(s))
     ax2 = axes(A)[2]
     xi = zeros(T, ax2)
@@ -319,7 +317,7 @@ function _generic_mapnum_mul!(Y, A, s)
     end
     return Y
 end
-function _generic_mapnum_mul!(Y, A, s, α, β)
+function _generic_map_mul!(Y, A, s::Number, α, β)
     T = promote_type(eltype(A), typeof(s))
     ax2 = axes(A)[2]
     xi = zeros(T, ax2)

src/blockmap.jl

@@ -26,7 +26,7 @@ MulStyle(A::BlockMap) = MulStyle(A.maps...)
 function _getranges(maps, dim, inds=1:length(maps))
     ends = map(i -> size(maps[i], dim)::Int, inds)
     cumsum!(ends, ends)
-    starts = vcat(1, 1 .+ @views ends[1:end-1])
+    starts = vcat(1, 1 .+ @views ends[firstindex(ends):lastindex(ends)-1])
     return UnitRange.(starts, ends)
 end
 
@@ -172,7 +172,7 @@ function Base.hvcat(rows::Tuple{Vararg{Int}},
         throw(ArgumentError("mismatch between row sizes and number of arguments"))
     n = fill(-1, length(As))
     j = 0
-    for i in 1:nr # infer UniformScaling sizes from row counts, if possible:
+    for i in eachindex(rows) # infer UniformScaling sizes from row counts, if possible:
         ni = -1 # number of rows in this block-row, -1 indicates unknown
         for k in 1:rows[i]
             if !isa(As[j+k], UniformScaling)
@@ -190,10 +190,10 @@ function Base.hvcat(rows::Tuple{Vararg{Int}},
     # check for consistent total column number
     nc = -1
     j = 0
-    for i in 1:nr
+    for i in eachindex(rows)
         nci = 0
         rows[i] > 0 && n[j+1] == -1 && (j += rows[i]; continue)
-        for k = 1:rows[i]
+        for k in 1:rows[i]
             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"))
@@ -202,11 +202,11 @@ function Base.hvcat(rows::Tuple{Vararg{Int}},
     end
     nc == -1 && throw(ArgumentError("sizes of UniformScalings could not be inferred"))
     j = 0
-    for i in 1:nr
+    for i in eachindex(rows)
         if rows[i] > 0 && n[j+1] == -1 # this row consists entirely of UniformScalings
             nci, r = divrem(nc, rows[i])
             r != 0 && throw(DimensionMismatch("indivisible UniformScaling sizes"))
-            for k = 1:rows[i]
+            for k in 1:rows[i]
                 n[j+k] = nci
             end
         end

src/conversion.jl

@@ -1,12 +1,12 @@
 # Matrix: create matrix representation of LinearMap
 function Base.Matrix{T}(A::LinearMap) where {T}
     mat = Matrix{T}(undef, size(A))
-    mul!(mat, A, one(T))
+    _unsafe_mul!(mat, A, true)
 end
 
 function Base.AbstractMatrix{T}(A::LinearMap) where {T}
     mat = similar(Array{T}, axes(A))
-    mul!(mat, A, one(T))
+    _unsafe_mul!(mat, A, true)
 end
 
 Base.Matrix(A::LinearMap{T}) where {T} = Matrix{T}(A)
@@ -28,7 +28,7 @@ function SparseArrays.sparse(A::LinearMap{T}) where {T}
     v = fill(zero(T), N)
     Av = Vector{T}(undef, M)
 
-    @inbounds for i in 1:N
+    @inbounds for i in eachindex(v)
         v[i] = one(T)
         _unsafe_mul!(Av, A, v)
         js = findall(!iszero, Av)

src/transpose.jl

@@ -50,44 +50,24 @@ Base.:(==)(A::LinearMap, B::TransposeMap)    = issymmetric(A) && B.lmap == A
 Base.:(==)(A::LinearMap, B::AdjointMap)      = ishermitian(A) && B.lmap == A
 
 # multiplication with vector/matrices
-# # TransposeMap
-_unsafe_mul!(y::AbstractVecOrMat, A::TransposeMap, x::AbstractVector) =
-    issymmetric(A.lmap) ?
-        _unsafe_mul!(y, A.lmap, x) : error("transpose not implemented for $(A.lmap)")
-_unsafe_mul!(y::AbstractMatrix, A::TransposeMap, x::AbstractMatrix) =
-    issymmetric(A.lmap) ?
-        _unsafe_mul!(y, A.lmap, x) : _generic_mapmat_mul!(y, A, x)
-_unsafe_mul!(y::AbstractMatrix, A::TransposeMap, x::Number) =
-    issymmetric(A.lmap) ?
-        _unsafe_mul!(y, A.lmap, x) : _generic_mapnum_mul!(y, A, x)
-_unsafe_mul!(y::AbstractVecOrMat, A::TransposeMap, x::AbstractVector, α::Number, β::Number)=
-    issymmetric(A.lmap) ?
-        _unsafe_mul!(y, A.lmap, x, α, β) : _generic_mapvec_mul!(y, A, x, α, β)
-_unsafe_mul!(y::AbstractMatrix, A::TransposeMap, x::AbstractMatrix, α::Number, β::Number) =
-    issymmetric(A.lmap) ?
-        _unsafe_mul!(y, A.lmap, x, α, β) : _generic_mapmat_mul!(y, A, x, α, β)
-_unsafe_mul!(y::AbstractMatrix, A::TransposeMap, x::Number, α::Number, β::Number) =
-    issymmetric(A.lmap) ?
-        _unsafe_mul!(y, A.lmap, x, α, β) : _generic_mapnum_mul!(y, A, x, α, β)
-# # AdjointMap
-_unsafe_mul!(y::AbstractVecOrMat, A::AdjointMap, x::AbstractVector) =
-    ishermitian(A.lmap) ?
-        _unsafe_mul!(y, A.lmap, x) : error("adjoint not implemented for $(A.lmap)")
-_unsafe_mul!(y::AbstractMatrix, A::AdjointMap, x::AbstractMatrix) =
-    ishermitian(A.lmap) ?
-        _unsafe_mul!(y, A.lmap, x) : _generic_mapmat_mul!(y, A, x)
-_unsafe_mul!(y::AbstractMatrix, A::AdjointMap, x::Number) =
-    ishermitian(A.lmap) ?
-        _unsafe_mul!(y, A.lmap, x) : _generic_mapnum_mul!(y, A, x)
-_unsafe_mul!(y::AbstractVecOrMat, A::AdjointMap, x::AbstractVector, α::Number, β::Number) =
-    ishermitian(A.lmap) ?
-        _unsafe_mul!(y, A.lmap, x, α, β) : _generic_mapvec_mul!(y, A, x, α, β)
-_unsafe_mul!(y::AbstractMatrix, A::AdjointMap, x::AbstractMatrix, α::Number, β::Number) =
-    ishermitian(A.lmap) ?
-        _unsafe_mul!(y, A.lmap, x, α, β) : _generic_mapmat_mul!(y, A, x, α, β)
-_unsafe_mul!(y::AbstractMatrix, A::AdjointMap, x::Number, α::Number, β::Number) =
-    ishermitian(A.lmap) ?
-        _unsafe_mul!(y, A.lmap, x, α, β) : _generic_mapnum_mul!(y, A, x, α, β)
+for (Typ, prop, text) in ((AdjointMap, ishermitian, "adjoint"), (TransposeMap, issymmetric, "transpose"))
+    @eval _unsafe_mul!(y::AbstractVecOrMat, A::$Typ, x::AbstractVector) =
+        $prop(A.lmap) ?
+            _unsafe_mul!(y, A.lmap, x) : error($text * " not implemented for $(A.lmap)")
+    @eval _unsafe_mul!(y::AbstractVecOrMat, A::$Typ, x::AbstractVector, α::Number, β::Number) =
+        $prop(A.lmap) ?
+            _unsafe_mul!(y, A.lmap, x, α, β) : _generic_map_mul!(y, A, x, α, β)
+
+    for In in (Number, AbstractMatrix)
+        @eval _unsafe_mul!(y::AbstractMatrix, A::$Typ, x::$In) =
+            $prop(A.lmap) ?
+                _unsafe_mul!(y, A.lmap, x) : _generic_map_mul!(y, A, x)
+
+        @eval _unsafe_mul!(y::AbstractMatrix, A::$Typ, x::$In, α::Number, β::Number) =
+            ishermitian(A.lmap) ?
+                _unsafe_mul!(y, A.lmap, x, α, β) : _generic_map_mul!(y, A, x, α, β)
+    end
+end
 
 # # ConjugateMap
 const ConjugateMap = AdjointMap{<:Any, <:TransposeMap}
@@ -104,9 +84,10 @@ for (In, Out) in ((AbstractVector, AbstractVecOrMat), (AbstractMatrix, AbstractM
         end
     end
 end
+_unsafe_mul!(y::AbstractMatrix, Ac::ConjugateMap, x::Number) = _conjmul!(y, Ac.lmap.lmap, x)
 
 # multiplication helper function
-_conjmul!(y, A, x) = conj!(mul!(y, A, conj(x)))
+_conjmul!(y, A, x) = conj!(_unsafe_mul!(y, A, conj(x)))
 function _conjmul!(y, A, x::AbstractVector, α, β)
     xca = conj!(x * α)
     z = A * xca