extend custom documentation

Author: dkarrasch

docs/src/custom.jl

@@ -229,3 +229,48 @@ 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`.

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)