From 1a52dea017c0d316025c78987fc969510bc51484 Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <Daniel.Karrasch@gmx.de>
Date: Thu, 2 Dec 2021 14:41:50 +0100
Subject: [PATCH 01/23] Add `getindex` functionality

---
 src/LinearMaps.jl  |   1 +
 src/getindex.jl    | 132 +++++++++++++++++++++++++++++++++++++++++++++
 test/getindex.jl   |  67 +++++++++++++++++++++++
 test/linearmaps.jl |  64 +++++++++++-----------
 test/runtests.jl   |   2 +
 5 files changed, 234 insertions(+), 32 deletions(-)
 create mode 100644 src/getindex.jl
 create mode 100644 test/getindex.jl

diff --git a/src/LinearMaps.jl b/src/LinearMaps.jl
index 52504975..80762c73 100644
--- a/src/LinearMaps.jl
+++ b/src/LinearMaps.jl
@@ -345,6 +345,7 @@ include("fillmap.jl") # linear maps representing constantly filled matrices
 include("embeddedmap.jl") # embedded linear maps
 include("conversion.jl") # conversion of linear maps to matrices
 include("show.jl") # show methods for LinearMap objects
+include("getindex.jl") # getindex functionality
 
 """
     LinearMap(A::LinearMap; kwargs...)::WrappedMap
diff --git a/src/getindex.jl b/src/getindex.jl
new file mode 100644
index 00000000..5c46c636
--- /dev/null
+++ b/src/getindex.jl
@@ -0,0 +1,132 @@
+# module GetIndex
+
+# using ..LinearMaps: LinearMap, AdjointMap, TransposeMap, FillMap, LinearCombination,
+#     ScaledMap, UniformScalingMap, WrappedMap
+
+# required in Base.to_indices for [:]-indexing
+Base.eachindex(::IndexLinear, A::LinearMap) = (Base.@_inline_meta; Base.OneTo(length(A)))
+Base.lastindex(A::LinearMap) = (Base.@_inline_meta; last(eachindex(IndexLinear(), A)))
+Base.firstindex(A::LinearMap) = (Base.@_inline_meta; first(eachindex(IndexLinear(), A)))
+
+function Base.checkbounds(A::LinearMap, i, j)
+    Base.@_inline_meta
+    Base.checkbounds_indices(Bool, axes(A), (i, j)) || throw(BoundsError(A, (i, j)))
+    nothing
+end
+# Linear indexing is explicitly allowed when there is only one (non-cartesian) index
+function Base.checkbounds(A::LinearMap, i)
+    Base.@_inline_meta
+    Base.checkindex(Bool, Base.OneTo(length(A)), i) || throw(BoundsError(A, i))
+    nothing
+end
+
+# main entry point
+Base.@propagate_inbounds function Base.getindex(A::LinearMap, I...)
+    # TODO: introduce some sort of switch?
+    Base.@_inline_meta
+    @boundscheck checkbounds(A, I...)
+    @inbounds _getindex(A, Base.to_indices(A, I)...)
+end
+# quick pass forward
+Base.@propagate_inbounds Base.getindex(A::ScaledMap, I...) = A.λ .* getindex(A.lmap, I...)
+Base.@propagate_inbounds Base.getindex(A::WrappedMap, I...) = A.lmap[I...]
+Base.@propagate_inbounds Base.getindex(A::WrappedMap, i::Integer) = A.lmap[i]
+Base.@propagate_inbounds Base.getindex(A::WrappedMap, i::Integer, j::Integer) = A.lmap[i,j]
+
+########################
+# linear indexing
+########################
+Base.@propagate_inbounds function _getindex(A::LinearMap, i::Integer)
+    Base.@_inline_meta
+    i1, i2 = Base._ind2sub(axes(A), i)
+    return _getindex(A, i1, i2)
+end
+Base.@propagate_inbounds _getindex(A::LinearMap, I::AbstractVector{<:Integer}) =
+    [_getindex(A, i) for i in I]
+_getindex(A::LinearMap, ::Base.Slice) = vec(Matrix(A))
+
+########################
+# Cartesian indexing
+########################
+Base.@propagate_inbounds _getindex(A::LinearMap, i::Integer, j::Integer) =
+    _getindex(A, Base.Slice(axes(A)[1]), j)[i]
+Base.@propagate_inbounds function _getindex(A::LinearMap, i::Integer, J::AbstractVector{<:Integer})
+    try
+        return (basevec(A, i)'A)[J]
+    catch
+        x = zeros(eltype(A), size(A, 2))
+        y = similar(x, eltype(A), size(A, 1))
+        r = similar(x, eltype(A), length(J))
+        for (ind, j) in enumerate(J)
+            x[j] = one(eltype(A))
+            _unsafe_mul!(y, A, x)
+            r[ind] = y[i]
+            x[j] = zero(eltype(A))
+        end
+        return r
+    end
+end
+function _getindex(A::LinearMap, i::Integer, J::Base.Slice)
+    try
+        return vec(basevec(A, i)'A)
+    catch
+        return vec(_getindex(A, i:i, J))
+    end
+end
+Base.@propagate_inbounds _getindex(A::LinearMap, I::AbstractVector{<:Integer}, j::Integer) =
+    _getindex(A, Base.Slice(axes(A)[1]), j)[I] # = A[:,j][I] w/o bounds check
+_getindex(A::LinearMap, ::Base.Slice, j::Integer) = A*basevec(A, j)
+Base.@propagate_inbounds function _getindex(A::LinearMap, Is::Vararg{AbstractVector{<:Integer},2})
+    shape = Base.index_shape(Is...)
+    dest = zeros(eltype(A), shape)
+    I, J = Is
+    for (ind, ij) in zip(eachindex(dest), Iterators.product(I, J))
+        i, j = ij
+        dest[ind] = _getindex(A, i, j)
+    end
+    return dest
+end
+Base.@propagate_inbounds function _getindex(A::LinearMap, I::AbstractVector{<:Integer}, ::Base.Slice)
+    x = zeros(eltype(A), size(A, 2))
+    y = similar(x, eltype(A), size(A, 1))
+    r = similar(x, eltype(A), (length(I), size(A, 2)))
+    @views for j in axes(A)[2]
+        x[j] = one(eltype(A))
+        _unsafe_mul!(y, A, x)
+        r[:,j] .= y[I]
+        x[j] = zero(eltype(A))
+    end
+    return r
+end
+Base.@propagate_inbounds function _getindex(A::LinearMap, ::Base.Slice, J::AbstractVector{<:Integer})
+    x = zeros(eltype(A), size(A, 2))
+    y = similar(x, eltype(A), (size(A, 1), length(J)))
+    for (i, j) in enumerate(J)
+        x[j] = one(eltype(A))
+        _unsafe_mul!(selectdim(y, 2, i), A, x)
+        x[j] = zero(eltype(A))
+    end
+    return y
+end
+_getindex(A::LinearMap, ::Base.Slice, ::Base.Slice) = Matrix(A)
+
+# specialized methods
+_getindex(A::FillMap, ::Integer, ::Integer) = A.λ
+Base.@propagate_inbounds _getindex(A::LinearCombination, i::Integer, j::Integer) =
+    sum(a -> A.maps[a][i, j], eachindex(A.maps))
+Base.@propagate_inbounds _getindex(A::AdjointMap, i::Integer, j::Integer) =
+    adjoint(A.lmap[j, i])
+Base.@propagate_inbounds _getindex(A::TransposeMap, i::Integer, j::Integer) =
+    transpose(A.lmap[j, i])
+_getindex(A::UniformScalingMap, i::Integer, j::Integer) = ifelse(i == j, A.λ, zero(eltype(A)))
+
+# helpers
+function basevec(A, i::Integer)
+    x = zeros(eltype(A), size(A, 2))
+    @inbounds x[i] = one(eltype(A))
+    return x
+end
+
+nogetindex_error() = error("indexing not allowed for LinearMaps; consider setting `LinearMaps.allowgetindex = true`")
+
+# end # module
diff --git a/test/getindex.jl b/test/getindex.jl
new file mode 100644
index 00000000..71cd4849
--- /dev/null
+++ b/test/getindex.jl
@@ -0,0 +1,67 @@
+using BenchmarkTools, LinearAlgebra, LinearMaps, Test
+# using LinearMaps.GetIndex
+
+function test_getindex(A::LinearMap, M::AbstractMatrix)
+    @assert size(A) == size(M)
+    @test all((A[i,j] == M[i,j] for i in axes(A, 1), j in axes(A, 2)))
+    @test all((A[i] == M[i] for i in 1:length(A)))
+    @test A[1,1] == M[1,1]
+    @test A[:] == M[:]
+    @test A[1,:] == M[1,:]
+    @test A[:,1] == M[:,1]
+    @test A[1:4,:] == M[1:4,:]
+    @test A[:,1:4] == M[:,1:4]
+    @test A[1,1:3] == M[1,1:3]
+    @test A[1:3,1] == M[1:3,1]
+    @test A[2:end,1] == M[2:end,1]
+    @test A[1:2,1:3] == M[1:2,1:3]
+    @test A[[2,1],1:3] == M[[2,1],1:3]
+    @test A[:,:] == M
+    @test A[7] == M[7]
+    @test_throws BoundsError A[firstindex(A)-1]
+    @test_throws BoundsError A[lastindex(A)+1]
+    @test_throws BoundsError A[6,1]
+    @test_throws BoundsError A[1,6]
+    @test_throws BoundsError A[2,1:6]
+    @test_throws BoundsError A[1:6,2]
+    return true
+end
+
+@testset "getindex" begin
+    A = rand(5,5)
+    L = LinearMap(A)
+    # @btime getindex($A, i) setup=(i = rand(1:9));
+    # @btime getindex($L, i) setup=(i = rand(1:9));
+    # @btime (getindex($A, i, j)) setup=(i = rand(1:3); j = rand(1:3));
+    # @btime (getindex($L, i, j)) setup=(i = rand(1:3); j = rand(1:3));
+
+    struct TwoMap <: LinearMaps.LinearMap{Float64} end
+    Base.size(::TwoMap) = (5,5)
+    LinearMaps._getindex(::TwoMap, i::Integer, j::Integer) = 2.0
+    LinearMaps._unsafe_mul!(y::AbstractVector, ::TwoMap, x::AbstractVector) = fill!(y, 2.0*sum(x))
+
+    @test test_getindex(TwoMap(), fill(2.0, 5, 5))
+    Base.adjoint(A::TwoMap) = A
+    @test test_getindex(TwoMap(), fill(2.0, 5, 5))
+
+    MA = rand(ComplexF64, 5, 5)
+    for FA in (
+        LinearMap{ComplexF64}((y, x) -> mul!(y, MA, x), (y, x) -> mul!(y, MA', x), 5, 5),
+        LinearMap{ComplexF64}((y, x) -> mul!(y, MA, x), 5, 5),
+    )
+        @test test_getindex(FA, MA)
+        @test test_getindex(3FA, 3MA)
+        @test test_getindex(FA + FA, 2MA)
+        if !isnothing(FA.fc)
+            @test test_getindex(transpose(FA), transpose(MA))
+            @test test_getindex(transpose(3FA), transpose(3MA))
+            @test test_getindex(3transpose(FA), transpose(3MA))
+            @test test_getindex(adjoint(FA), adjoint(MA))
+            @test test_getindex(adjoint(3FA), adjoint(3MA))
+            @test test_getindex(3adjoint(FA), adjoint(3MA))
+        end
+    end
+
+    @test test_getindex(FillMap(0.5, (5, 5)), fill(0.5, (5, 5)))
+    @test test_getindex(LinearMap(0.5I, 5), Matrix(0.5I, 5, 5))
+end
diff --git a/test/linearmaps.jl b/test/linearmaps.jl
index 2783cfe6..22c4efbf 100644
--- a/test/linearmaps.jl
+++ b/test/linearmaps.jl
@@ -33,20 +33,20 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays
     end
 end
 
-# new type
-struct SimpleFunctionMap <: LinearMap{Float64}
-    f::Function
-    N::Int
-end
-struct SimpleComplexFunctionMap <: LinearMap{Complex{Float64}}
-    f::Function
-    N::Int
-end
-Base.size(A::Union{SimpleFunctionMap,SimpleComplexFunctionMap}) = (A.N, A.N)
-Base.:(*)(A::Union{SimpleFunctionMap,SimpleComplexFunctionMap}, v::AbstractVector) = A.f(v)
-LinearAlgebra.mul!(y::AbstractVector, A::Union{SimpleFunctionMap,SimpleComplexFunctionMap}, x::AbstractVector) = copyto!(y, *(A, x))
-
 @testset "new LinearMap type" begin
+    # new type
+    struct SimpleFunctionMap <: LinearMap{Float64}
+        f::Function
+        N::Int
+    end
+    struct SimpleComplexFunctionMap <: LinearMap{Complex{Float64}}
+        f::Function
+        N::Int
+    end
+    Base.size(A::Union{SimpleFunctionMap,SimpleComplexFunctionMap}) = (A.N, A.N)
+    Base.:(*)(A::Union{SimpleFunctionMap,SimpleComplexFunctionMap}, v::AbstractVector) = A.f(v)
+    LinearAlgebra.mul!(y::AbstractVector, A::Union{SimpleFunctionMap,SimpleComplexFunctionMap}, x::AbstractVector) = copyto!(y, *(A, x))
+
     F = SimpleFunctionMap(cumsum, 10)
     FC = SimpleComplexFunctionMap(cumsum, 10)
     @test @inferred ndims(F) == 2
@@ -84,27 +84,27 @@ LinearAlgebra.mul!(y::AbstractVector, A::Union{SimpleFunctionMap,SimpleComplexFu
     @test Fs isa SparseMatrixCSC
 end
 
-struct MyFillMap{T} <: LinearMaps.LinearMap{T}
-    λ::T
-    size::Dims{2}
-    function MyFillMap(λ::T, dims::Dims{2}) where {T}
-        all(d -> d >= 0, dims) || throw(ArgumentError("dims of MyFillMap must be non-negative"))
-        promote_type(T, typeof(λ)) == T || throw(InexactError())
-        return new{T}(λ, dims)
+@testset "transpose of new LinearMap type" begin
+    struct MyFillMap{T} <: LinearMaps.LinearMap{T}
+        λ::T
+        size::Dims{2}
+        function MyFillMap(λ::T, dims::Dims{2}) where {T}
+            all(d -> d >= 0, dims) || throw(ArgumentError("dims of MyFillMap must be non-negative"))
+            promote_type(T, typeof(λ)) == T || throw(InexactError())
+            return new{T}(λ, dims)
+        end
+    end
+    Base.size(A::MyFillMap) = A.size
+    function LinearAlgebra.mul!(y::AbstractVecOrMat, A::MyFillMap, x::AbstractVector)
+        LinearMaps.check_dim_mul(y, A, x)
+        return fill!(y, iszero(A.λ) ? zero(eltype(y)) : A.λ*sum(x))
+    end
+    function LinearAlgebra.mul!(y::AbstractVecOrMat, transA::LinearMaps.TransposeMap{<:Any,<:MyFillMap}, x::AbstractVector)
+        LinearMaps.check_dim_mul(y, transA, x)
+        λ = transA.lmap.λ
+        return fill!(y, iszero(λ) ? zero(eltype(y)) : transpose(λ)*sum(x))
     end
-end
-Base.size(A::MyFillMap) = A.size
-function LinearAlgebra.mul!(y::AbstractVecOrMat, A::MyFillMap, x::AbstractVector)
-    LinearMaps.check_dim_mul(y, A, x)
-    return fill!(y, iszero(A.λ) ? zero(eltype(y)) : A.λ*sum(x))
-end
-function LinearAlgebra.mul!(y::AbstractVecOrMat, transA::LinearMaps.TransposeMap{<:Any,<:MyFillMap}, x::AbstractVector)
-    LinearMaps.check_dim_mul(y, transA, x)
-    λ = transA.lmap.λ
-    return fill!(y, iszero(λ) ? zero(eltype(y)) : transpose(λ)*sum(x))
-end
 
-@testset "transpose of new LinearMap type" begin
     A = MyFillMap(5.0, (3, 3))
     x = ones(3)
     @test A * x == fill(15.0, 3)
diff --git a/test/runtests.jl b/test/runtests.jl
index cd7065c3..a152bc0f 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -35,3 +35,5 @@ include("fillmap.jl")
 include("nontradaxes.jl")
 
 include("embeddedmap.jl")
+
+include("getindex.jl")

From 97c1fb3f2b05ee0dc4e90311327a590f596d443c Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <Daniel.Karrasch@gmx.de>
Date: Fri, 3 Dec 2021 18:28:10 +0100
Subject: [PATCH 02/23] increase code coverage

---
 src/getindex.jl  | 2 +-
 test/getindex.jl | 2 ++
 2 files changed, 3 insertions(+), 1 deletion(-)

diff --git a/src/getindex.jl b/src/getindex.jl
index 5c46c636..74934b9d 100644
--- a/src/getindex.jl
+++ b/src/getindex.jl
@@ -127,6 +127,6 @@ function basevec(A, i::Integer)
     return x
 end
 
-nogetindex_error() = error("indexing not allowed for LinearMaps; consider setting `LinearMaps.allowgetindex = true`")
+# nogetindex_error() = error("indexing not allowed for LinearMaps; consider setting `LinearMaps.allowgetindex = true`")
 
 # end # module
diff --git a/test/getindex.jl b/test/getindex.jl
index 71cd4849..eec4954f 100644
--- a/test/getindex.jl
+++ b/test/getindex.jl
@@ -18,6 +18,7 @@ function test_getindex(A::LinearMap, M::AbstractMatrix)
     @test A[[2,1],1:3] == M[[2,1],1:3]
     @test A[:,:] == M
     @test A[7] == M[7]
+    @test A[3:7] == M[3:7]
     @test_throws BoundsError A[firstindex(A)-1]
     @test_throws BoundsError A[lastindex(A)+1]
     @test_throws BoundsError A[6,1]
@@ -30,6 +31,7 @@ end
 @testset "getindex" begin
     A = rand(5,5)
     L = LinearMap(A)
+    @test test_getindex(L, A)
     # @btime getindex($A, i) setup=(i = rand(1:9));
     # @btime getindex($L, i) setup=(i = rand(1:9));
     # @btime (getindex($A, i, j)) setup=(i = rand(1:3); j = rand(1:3));

From d56060f5ec67f9d3714d4d97040d45725e012338 Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <Daniel.Karrasch@gmx.de>
Date: Fri, 3 Dec 2021 23:06:29 +0100
Subject: [PATCH 03/23] simplify compiler macros

---
 src/getindex.jl | 55 +++++++++++++++++++++----------------------------
 1 file changed, 24 insertions(+), 31 deletions(-)

diff --git a/src/getindex.jl b/src/getindex.jl
index 74934b9d..4794eff6 100644
--- a/src/getindex.jl
+++ b/src/getindex.jl
@@ -4,31 +4,28 @@
 #     ScaledMap, UniformScalingMap, WrappedMap
 
 # required in Base.to_indices for [:]-indexing
-Base.eachindex(::IndexLinear, A::LinearMap) = (Base.@_inline_meta; Base.OneTo(length(A)))
-Base.lastindex(A::LinearMap) = (Base.@_inline_meta; last(eachindex(IndexLinear(), A)))
-Base.firstindex(A::LinearMap) = (Base.@_inline_meta; first(eachindex(IndexLinear(), A)))
+Base.eachindex(::IndexLinear, A::LinearMap) = Base.OneTo(length(A))
+Base.lastindex(A::LinearMap) = last(eachindex(IndexLinear(), A))
+Base.firstindex(A::LinearMap) = first(eachindex(IndexLinear(), A))
 
 function Base.checkbounds(A::LinearMap, i, j)
-    Base.@_inline_meta
     Base.checkbounds_indices(Bool, axes(A), (i, j)) || throw(BoundsError(A, (i, j)))
     nothing
 end
 # Linear indexing is explicitly allowed when there is only one (non-cartesian) index
 function Base.checkbounds(A::LinearMap, i)
-    Base.@_inline_meta
     Base.checkindex(Bool, Base.OneTo(length(A)), i) || throw(BoundsError(A, i))
     nothing
 end
 
 # main entry point
-Base.@propagate_inbounds function Base.getindex(A::LinearMap, I...)
+function Base.getindex(A::LinearMap, I...)
     # TODO: introduce some sort of switch?
-    Base.@_inline_meta
     @boundscheck checkbounds(A, I...)
-    @inbounds _getindex(A, Base.to_indices(A, I)...)
+    _getindex(A, Base.to_indices(A, I)...)
 end
 # quick pass forward
-Base.@propagate_inbounds Base.getindex(A::ScaledMap, I...) = A.λ .* getindex(A.lmap, I...)
+Base.@propagate_inbounds Base.getindex(A::ScaledMap, I...) = A.λ .* A.lmap[I...]
 Base.@propagate_inbounds Base.getindex(A::WrappedMap, I...) = A.lmap[I...]
 Base.@propagate_inbounds Base.getindex(A::WrappedMap, i::Integer) = A.lmap[i]
 Base.@propagate_inbounds Base.getindex(A::WrappedMap, i::Integer, j::Integer) = A.lmap[i,j]
@@ -36,28 +33,26 @@ Base.@propagate_inbounds Base.getindex(A::WrappedMap, i::Integer, j::Integer) =
 ########################
 # linear indexing
 ########################
-Base.@propagate_inbounds function _getindex(A::LinearMap, i::Integer)
-    Base.@_inline_meta
+function _getindex(A::LinearMap, i::Integer)
     i1, i2 = Base._ind2sub(axes(A), i)
     return _getindex(A, i1, i2)
 end
-Base.@propagate_inbounds _getindex(A::LinearMap, I::AbstractVector{<:Integer}) =
-    [_getindex(A, i) for i in I]
+_getindex(A::LinearMap, I::AbstractVector{<:Integer}) = [_getindex(A, i) for i in I]
 _getindex(A::LinearMap, ::Base.Slice) = vec(Matrix(A))
 
 ########################
 # Cartesian indexing
 ########################
-Base.@propagate_inbounds _getindex(A::LinearMap, i::Integer, j::Integer) =
-    _getindex(A, Base.Slice(axes(A)[1]), j)[i]
-Base.@propagate_inbounds function _getindex(A::LinearMap, i::Integer, J::AbstractVector{<:Integer})
+_getindex(A::LinearMap, i::Integer, j::Integer) =
+    @inbounds _getindex(A, Base.Slice(axes(A)[1]), j)[i]
+function _getindex(A::LinearMap, i::Integer, J::AbstractVector{<:Integer})
     try
-        return (basevec(A, i)'A)[J]
+        return @inbounds (basevec(A, i)'A)[J]
     catch
         x = zeros(eltype(A), size(A, 2))
         y = similar(x, eltype(A), size(A, 1))
         r = similar(x, eltype(A), length(J))
-        for (ind, j) in enumerate(J)
+        @inbounds for (ind, j) in enumerate(J)
             x[j] = one(eltype(A))
             _unsafe_mul!(y, A, x)
             r[ind] = y[i]
@@ -73,10 +68,10 @@ function _getindex(A::LinearMap, i::Integer, J::Base.Slice)
         return vec(_getindex(A, i:i, J))
     end
 end
-Base.@propagate_inbounds _getindex(A::LinearMap, I::AbstractVector{<:Integer}, j::Integer) =
-    _getindex(A, Base.Slice(axes(A)[1]), j)[I] # = A[:,j][I] w/o bounds check
+_getindex(A::LinearMap, I::AbstractVector{<:Integer}, j::Integer) =
+    @inbounds _getindex(A, Base.Slice(axes(A)[1]), j)[I] # = A[:,j][I] w/o bounds check
 _getindex(A::LinearMap, ::Base.Slice, j::Integer) = A*basevec(A, j)
-Base.@propagate_inbounds function _getindex(A::LinearMap, Is::Vararg{AbstractVector{<:Integer},2})
+function _getindex(A::LinearMap, Is::Vararg{AbstractVector{<:Integer},2})
     shape = Base.index_shape(Is...)
     dest = zeros(eltype(A), shape)
     I, J = Is
@@ -86,11 +81,11 @@ Base.@propagate_inbounds function _getindex(A::LinearMap, Is::Vararg{AbstractVec
     end
     return dest
 end
-Base.@propagate_inbounds function _getindex(A::LinearMap, I::AbstractVector{<:Integer}, ::Base.Slice)
+function _getindex(A::LinearMap, I::AbstractVector{<:Integer}, ::Base.Slice)
     x = zeros(eltype(A), size(A, 2))
     y = similar(x, eltype(A), size(A, 1))
     r = similar(x, eltype(A), (length(I), size(A, 2)))
-    @views for j in axes(A)[2]
+    @inbounds @views for j in axes(A)[2]
         x[j] = one(eltype(A))
         _unsafe_mul!(y, A, x)
         r[:,j] .= y[I]
@@ -98,10 +93,10 @@ Base.@propagate_inbounds function _getindex(A::LinearMap, I::AbstractVector{<:In
     end
     return r
 end
-Base.@propagate_inbounds function _getindex(A::LinearMap, ::Base.Slice, J::AbstractVector{<:Integer})
+function _getindex(A::LinearMap, ::Base.Slice, J::AbstractVector{<:Integer})
     x = zeros(eltype(A), size(A, 2))
     y = similar(x, eltype(A), (size(A, 1), length(J)))
-    for (i, j) in enumerate(J)
+    @inbounds for (i, j) in enumerate(J)
         x[j] = one(eltype(A))
         _unsafe_mul!(selectdim(y, 2, i), A, x)
         x[j] = zero(eltype(A))
@@ -112,12 +107,10 @@ _getindex(A::LinearMap, ::Base.Slice, ::Base.Slice) = Matrix(A)
 
 # specialized methods
 _getindex(A::FillMap, ::Integer, ::Integer) = A.λ
-Base.@propagate_inbounds _getindex(A::LinearCombination, i::Integer, j::Integer) =
-    sum(a -> A.maps[a][i, j], eachindex(A.maps))
-Base.@propagate_inbounds _getindex(A::AdjointMap, i::Integer, j::Integer) =
-    adjoint(A.lmap[j, i])
-Base.@propagate_inbounds _getindex(A::TransposeMap, i::Integer, j::Integer) =
-    transpose(A.lmap[j, i])
+_getindex(A::LinearCombination, i::Integer, j::Integer) =
+    sum(a -> (@inbounds A.maps[a][i, j]), eachindex(A.maps))
+_getindex(A::AdjointMap, i::Integer, j::Integer) = @inbounds adjoint(A.lmap[j, i])
+_getindex(A::TransposeMap, i::Integer, j::Integer) = @inbounds transpose(A.lmap[j, i])
 _getindex(A::UniformScalingMap, i::Integer, j::Integer) = ifelse(i == j, A.λ, zero(eltype(A)))
 
 # helpers

From 0c9529c328489ddc4816e2b0a4ac1e4516bd2ae7 Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <Daniel.Karrasch@gmx.de>
Date: Tue, 7 Dec 2021 15:11:03 +0100
Subject: [PATCH 04/23] reduce code repetition

---
 src/getindex.jl | 74 +++++++++++++++++++++----------------------------
 1 file changed, 31 insertions(+), 43 deletions(-)

diff --git a/src/getindex.jl b/src/getindex.jl
index 4794eff6..9b176f11 100644
--- a/src/getindex.jl
+++ b/src/getindex.jl
@@ -3,6 +3,8 @@
 # using ..LinearMaps: LinearMap, AdjointMap, TransposeMap, FillMap, LinearCombination,
 #     ScaledMap, UniformScalingMap, WrappedMap
 
+const Indexer = AbstractVector{<:Integer}
+
 # required in Base.to_indices for [:]-indexing
 Base.eachindex(::IndexLinear, A::LinearMap) = Base.OneTo(length(A))
 Base.lastindex(A::LinearMap) = last(eachindex(IndexLinear(), A))
@@ -37,72 +39,58 @@ function _getindex(A::LinearMap, i::Integer)
     i1, i2 = Base._ind2sub(axes(A), i)
     return _getindex(A, i1, i2)
 end
-_getindex(A::LinearMap, I::AbstractVector{<:Integer}) = [_getindex(A, i) for i in I]
+_getindex(A::LinearMap, I::Indexer) = [_getindex(A, i) for i in I]
 _getindex(A::LinearMap, ::Base.Slice) = vec(Matrix(A))
 
 ########################
 # Cartesian indexing
 ########################
-_getindex(A::LinearMap, i::Integer, j::Integer) =
-    @inbounds _getindex(A, Base.Slice(axes(A)[1]), j)[i]
-function _getindex(A::LinearMap, i::Integer, J::AbstractVector{<:Integer})
+_getindex(A::LinearMap, i::Integer,   j::Integer) = (@inbounds (A*basevec(A, j))[i])
+_getindex(A::LinearMap, I::Indexer,   j::Integer) = (@inbounds (A*basevec(A, j))[I])
+_getindex(A::LinearMap, ::Base.Slice, j::Integer) = A*basevec(A, j)
+function _getindex(A::LinearMap, i::Integer, J::Indexer)
     try
+        # requires adjoint action to be defined
         return @inbounds (basevec(A, i)'A)[J]
     catch
-        x = zeros(eltype(A), size(A, 2))
-        y = similar(x, eltype(A), size(A, 1))
-        r = similar(x, eltype(A), length(J))
-        @inbounds for (ind, j) in enumerate(J)
-            x[j] = one(eltype(A))
-            _unsafe_mul!(y, A, x)
-            r[ind] = y[i]
-            x[j] = zero(eltype(A))
-        end
-        return r
+        return _getrows(A, i, J)
     end
 end
 function _getindex(A::LinearMap, i::Integer, J::Base.Slice)
     try
+        # requires adjoint action to be defined
         return vec(basevec(A, i)'A)
     catch
-        return vec(_getindex(A, i:i, J))
-    end
-end
-_getindex(A::LinearMap, I::AbstractVector{<:Integer}, j::Integer) =
-    @inbounds _getindex(A, Base.Slice(axes(A)[1]), j)[I] # = A[:,j][I] w/o bounds check
-_getindex(A::LinearMap, ::Base.Slice, j::Integer) = A*basevec(A, j)
-function _getindex(A::LinearMap, Is::Vararg{AbstractVector{<:Integer},2})
-    shape = Base.index_shape(Is...)
-    dest = zeros(eltype(A), shape)
-    I, J = Is
-    for (ind, ij) in zip(eachindex(dest), Iterators.product(I, J))
-        i, j = ij
-        dest[ind] = _getindex(A, i, j)
+        return _getrows(A, i, J)
     end
-    return dest
 end
-function _getindex(A::LinearMap, I::AbstractVector{<:Integer}, ::Base.Slice)
-    x = zeros(eltype(A), size(A, 2))
-    y = similar(x, eltype(A), size(A, 1))
-    r = similar(x, eltype(A), (length(I), size(A, 2)))
-    @inbounds @views for j in axes(A)[2]
-        x[j] = one(eltype(A))
-        _unsafe_mul!(y, A, x)
-        r[:,j] .= y[I]
-        x[j] = zero(eltype(A))
+function _getindex(A::LinearMap, I::Indexer, J::Indexer)
+    if length(I) <= length(J)
+        try
+            # requires adjoint action to be defined
+            return vcat(map(i -> (@inbounds (basevec(A, i)'A)[1:1,J]), I)...)
+        catch
+            return _getrows(A, I, J)
+        end
+    else
+        return _getrows(A, I, J)
     end
-    return r
 end
-function _getindex(A::LinearMap, ::Base.Slice, J::AbstractVector{<:Integer})
+_getrows(A::LinearMap, I, J) = _getrows!(zeros(eltype(A), Base.index_shape(I, J)), A, I, J)
+function _getrows!(dest, A, i, J)
     x = zeros(eltype(A), size(A, 2))
-    y = similar(x, eltype(A), (size(A, 1), length(J)))
-    @inbounds for (i, j) in enumerate(J)
+    temp = similar(x, eltype(A), size(A, 1))
+    @views @inbounds for (ind, j) in enumerate(J)
         x[j] = one(eltype(A))
-        _unsafe_mul!(selectdim(y, 2, i), A, x)
+        _unsafe_mul!(temp, A, x)
+        _copyto!(dest, ind, temp, i)
         x[j] = zero(eltype(A))
     end
-    return y
+    return dest
 end
+@inline _copyto!(dest, ind, temp, i::Integer) = (@inbounds dest[ind] = temp[i])
+@inline _copyto!(dest, ind, temp, I::Indexer) =
+    (@views @inbounds dest[:,ind] .= temp[I])
 _getindex(A::LinearMap, ::Base.Slice, ::Base.Slice) = Matrix(A)
 
 # specialized methods

From 743e4565c33eca27b357694f841350851f04af07 Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <Daniel.Karrasch@gmx.de>
Date: Wed, 8 Dec 2021 18:47:48 +0100
Subject: [PATCH 05/23] fix & test rectangular case, further optimizations

---
 src/getindex.jl  | 76 +++++++++++++++++++++++++++++++-----------------
 test/getindex.jl | 19 ++++++------
 2 files changed, 60 insertions(+), 35 deletions(-)

diff --git a/src/getindex.jl b/src/getindex.jl
index 9b176f11..2d94cbad 100644
--- a/src/getindex.jl
+++ b/src/getindex.jl
@@ -45,52 +45,38 @@ _getindex(A::LinearMap, ::Base.Slice) = vec(Matrix(A))
 ########################
 # Cartesian indexing
 ########################
-_getindex(A::LinearMap, i::Integer,   j::Integer) = (@inbounds (A*basevec(A, j))[i])
-_getindex(A::LinearMap, I::Indexer,   j::Integer) = (@inbounds (A*basevec(A, j))[I])
-_getindex(A::LinearMap, ::Base.Slice, j::Integer) = A*basevec(A, j)
+_getindex(A::LinearMap, i::Union{Integer,Indexer}, j::Integer) = (@inbounds (A*basevec(A, 2, j))[i])
+_getindex(A::LinearMap, ::Base.Slice, j::Integer) = A*basevec(A, 2, j)
 function _getindex(A::LinearMap, i::Integer, J::Indexer)
     try
         # requires adjoint action to be defined
-        return @inbounds (basevec(A, i)'A)[J]
+        return @inbounds (basevec(A, 1, i)'A)[J]
     catch
-        return _getrows(A, i, J)
+        return _fillbycols!(zeros(eltype(A), Base.index_shape(i, J)), A, i, J)
     end
 end
 function _getindex(A::LinearMap, i::Integer, J::Base.Slice)
     try
         # requires adjoint action to be defined
-        return vec(basevec(A, i)'A)
+        return vec(basevec(A, 1, i)'A)
     catch
-        return _getrows(A, i, J)
+        return _fillbycols!(zeros(eltype(A), Base.index_shape(i, J)), A, i, J)
     end
 end
 function _getindex(A::LinearMap, I::Indexer, J::Indexer)
+    dest = zeros(eltype(A), Base.index_shape(I, J))
     if length(I) <= length(J)
         try
             # requires adjoint action to be defined
-            return vcat(map(i -> (@inbounds (basevec(A, i)'A)[1:1,J]), I)...)
+            _fillbyrows!(dest, A, I, J)
         catch
-            return _getrows(A, I, J)
+            _fillbycols!(dest, A, I, J)
         end
     else
-        return _getrows(A, I, J)
-    end
-end
-_getrows(A::LinearMap, I, J) = _getrows!(zeros(eltype(A), Base.index_shape(I, J)), A, I, J)
-function _getrows!(dest, A, i, J)
-    x = zeros(eltype(A), size(A, 2))
-    temp = similar(x, eltype(A), size(A, 1))
-    @views @inbounds for (ind, j) in enumerate(J)
-        x[j] = one(eltype(A))
-        _unsafe_mul!(temp, A, x)
-        _copyto!(dest, ind, temp, i)
-        x[j] = zero(eltype(A))
+        _fillbycols!(dest, A, I, J)
     end
     return dest
 end
-@inline _copyto!(dest, ind, temp, i::Integer) = (@inbounds dest[ind] = temp[i])
-@inline _copyto!(dest, ind, temp, I::Indexer) =
-    (@views @inbounds dest[:,ind] .= temp[I])
 _getindex(A::LinearMap, ::Base.Slice, ::Base.Slice) = Matrix(A)
 
 # specialized methods
@@ -102,12 +88,50 @@ _getindex(A::TransposeMap, i::Integer, j::Integer) = @inbounds transpose(A.lmap[
 _getindex(A::UniformScalingMap, i::Integer, j::Integer) = ifelse(i == j, A.λ, zero(eltype(A)))
 
 # helpers
-function basevec(A, i::Integer)
-    x = zeros(eltype(A), size(A, 2))
+function basevec(A, dim, i::Integer)
+    x = zeros(eltype(A), size(A, dim))
     @inbounds x[i] = one(eltype(A))
     return x
 end
 
+function _fillbyrows!(dest, A, I, J)
+    x = zeros(eltype(A), size(A, 1))
+    temp = similar(x, eltype(A), size(A, 2))
+    @views @inbounds for (ind, i) in enumerate(I)
+        x[i] = one(eltype(A))
+        _unsafe_mul!(temp, A', x)
+        _copyrow!(dest, ind, temp, J)
+        x[i] = zero(eltype(A))
+    end
+    return dest
+end
+function _fillbycols!(dest, A, i, J)
+    x = zeros(eltype(A), size(A, 2))
+    temp = similar(x, eltype(A), size(A, 1))
+    @views @inbounds for (ind, j) in enumerate(J)
+        x[j] = one(eltype(A))
+        _unsafe_mul!(temp, A, x)
+        _copycol!(dest, ind, temp, i)
+        x[j] = zero(eltype(A))
+    end
+    return dest
+end
+function _fillbycols!(dest, A, ::Base.Slice, J)
+    x = zeros(eltype(A), size(A, 2))
+    @views @inbounds for (ind, j) in enumerate(J)
+        x[j] = one(eltype(A))
+        _unsafe_mul!(selectdim(dest, 2, ind), A, x)
+        x[j] = zero(eltype(A))
+    end
+    return dest
+end
+
+@inline _copycol!(dest, ind, temp, i::Integer) = (@inbounds dest[ind] = temp[i])
+@inline _copycol!(dest, ind, temp, I::Indexer) =
+    (@views @inbounds dest[:,ind] .= temp[I])
+@inline _copyrow!(dest, ind, temp, J::Indexer) =
+    (@views @inbounds dest[ind,:] .= adjoint.(temp[J]))
+
 # nogetindex_error() = error("indexing not allowed for LinearMaps; consider setting `LinearMaps.allowgetindex = true`")
 
 # end # module
diff --git a/test/getindex.jl b/test/getindex.jl
index eec4954f..b7e49316 100644
--- a/test/getindex.jl
+++ b/test/getindex.jl
@@ -22,29 +22,30 @@ function test_getindex(A::LinearMap, M::AbstractMatrix)
     @test_throws BoundsError A[firstindex(A)-1]
     @test_throws BoundsError A[lastindex(A)+1]
     @test_throws BoundsError A[6,1]
-    @test_throws BoundsError A[1,6]
-    @test_throws BoundsError A[2,1:6]
+    @test_throws BoundsError A[1,7]
+    @test_throws BoundsError A[2,1:7]
     @test_throws BoundsError A[1:6,2]
     return true
 end
 
 @testset "getindex" begin
-    A = rand(5,5)
+    A = rand(4,6)
     L = LinearMap(A)
     @test test_getindex(L, A)
-    # @btime getindex($A, i) setup=(i = rand(1:9));
-    # @btime getindex($L, i) setup=(i = rand(1:9));
-    # @btime (getindex($A, i, j)) setup=(i = rand(1:3); j = rand(1:3));
-    # @btime (getindex($L, i, j)) setup=(i = rand(1:3); j = rand(1:3));
+    # @btime getindex($A, i) setup=(i = rand(1:24));
+    # @btime getindex($L, i) setup=(i = rand(1:24));
+    # @btime (getindex($A, i, j)) setup=(i = rand(1:4); j = rand(1:6));
+    # @btime (getindex($L, i, j)) setup=(i = rand(1:4); j = rand(1:6));
 
     struct TwoMap <: LinearMaps.LinearMap{Float64} end
     Base.size(::TwoMap) = (5,5)
     LinearMaps._getindex(::TwoMap, i::Integer, j::Integer) = 2.0
     LinearMaps._unsafe_mul!(y::AbstractVector, ::TwoMap, x::AbstractVector) = fill!(y, 2.0*sum(x))
 
-    @test test_getindex(TwoMap(), fill(2.0, 5, 5))
+    T = TwoMap()
+    @test test_getindex(TwoMap(), fill(2.0, size(T)))
     Base.adjoint(A::TwoMap) = A
-    @test test_getindex(TwoMap(), fill(2.0, 5, 5))
+    @test test_getindex(TwoMap(), fill(2.0, size(T)))
 
     MA = rand(ComplexF64, 5, 5)
     for FA in (

From f61f2a6e8b8c148b859b753873c35366228a7e61 Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <Daniel.Karrasch@posteo.de>
Date: Tue, 25 Jan 2022 10:47:43 +0100
Subject: [PATCH 06/23] minor cleanup

---
 src/getindex.jl | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/src/getindex.jl b/src/getindex.jl
index 2d94cbad..c14ca882 100644
--- a/src/getindex.jl
+++ b/src/getindex.jl
@@ -100,7 +100,7 @@ function _fillbyrows!(dest, A, I, J)
     @views @inbounds for (ind, i) in enumerate(I)
         x[i] = one(eltype(A))
         _unsafe_mul!(temp, A', x)
-        _copyrow!(dest, ind, temp, J)
+        dest[ind,:] .= adjoint.(temp[J])
         x[i] = zero(eltype(A))
     end
     return dest
@@ -129,8 +129,6 @@ end
 @inline _copycol!(dest, ind, temp, i::Integer) = (@inbounds dest[ind] = temp[i])
 @inline _copycol!(dest, ind, temp, I::Indexer) =
     (@views @inbounds dest[:,ind] .= temp[I])
-@inline _copyrow!(dest, ind, temp, J::Indexer) =
-    (@views @inbounds dest[ind,:] .= adjoint.(temp[J]))
 
 # nogetindex_error() = error("indexing not allowed for LinearMaps; consider setting `LinearMaps.allowgetindex = true`")
 

From 090d928889d683d50488b67ff6fc05b74f00f12d Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <Daniel.Karrasch@posteo.de>
Date: Mon, 21 Mar 2022 12:12:36 +0100
Subject: [PATCH 07/23] update Aqua.jl badge

---
 README.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/README.md b/README.md
index 3e4a2008..856c6ea1 100644
--- a/README.md
+++ b/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

From f546fe1663dd5fd4c9383b56588ef67d74a8502a Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <Daniel.Karrasch@posteo.de>
Date: Mon, 21 Mar 2022 13:34:48 +0100
Subject: [PATCH 08/23] fix test

---
 test/getindex.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/test/getindex.jl b/test/getindex.jl
index b7e49316..80b0572e 100644
--- a/test/getindex.jl
+++ b/test/getindex.jl
@@ -1,5 +1,6 @@
-using BenchmarkTools, LinearAlgebra, LinearMaps, Test
+using LinearAlgebra, LinearMaps, Test
 # using LinearMaps.GetIndex
+# using BenchmarkTools
 
 function test_getindex(A::LinearMap, M::AbstractMatrix)
     @assert size(A) == size(M)

From 4e87153dab9a4d91c45461b36f953844eeaa5375 Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <Daniel.Karrasch@posteo.de>
Date: Mon, 23 May 2022 19:28:42 +0200
Subject: [PATCH 09/23] add logical and diagonal indexing

---
 src/getindex.jl  | 31 +++++++++++++++++++++++++++++++
 test/getindex.jl | 11 ++++++++---
 2 files changed, 39 insertions(+), 3 deletions(-)

diff --git a/src/getindex.jl b/src/getindex.jl
index c14ca882..97ff1e37 100644
--- a/src/getindex.jl
+++ b/src/getindex.jl
@@ -5,6 +5,7 @@
 
 const Indexer = AbstractVector{<:Integer}
 
+Base.IndexStyle(::LinearMap) = IndexCartesian()
 # required in Base.to_indices for [:]-indexing
 Base.eachindex(::IndexLinear, A::LinearMap) = Base.OneTo(length(A))
 Base.lastindex(A::LinearMap) = last(eachindex(IndexLinear(), A))
@@ -19,6 +20,10 @@ function Base.checkbounds(A::LinearMap, i)
     Base.checkindex(Bool, Base.OneTo(length(A)), i) || throw(BoundsError(A, i))
     nothing
 end
+# checkbounds in indexing via CartesianIndex
+Base.checkbounds(A::LinearMap, i::Union{CartesianIndex{2}, AbstractVecOrMat{CartesianIndex{2}}}) =
+    Base.checkbounds_indices(Bool, axes(A), (i,))
+Base.checkbounds(A::LinearMap, I::AbstractArray{Bool,2}) = axes(A) == axes(I)
 
 # main entry point
 function Base.getindex(A::LinearMap, I...)
@@ -41,6 +46,7 @@ function _getindex(A::LinearMap, i::Integer)
 end
 _getindex(A::LinearMap, I::Indexer) = [_getindex(A, i) for i in I]
 _getindex(A::LinearMap, ::Base.Slice) = vec(Matrix(A))
+_getindex(A::LinearMap, I::Vector{CartesianIndex{2}}) = [(@inbounds A[i]) for i in I]
 
 ########################
 # Cartesian indexing
@@ -130,6 +136,31 @@ end
 @inline _copycol!(dest, ind, temp, I::Indexer) =
     (@views @inbounds dest[:,ind] .= temp[I])
 
+# diagonal indexing
+function LinearAlgebra.diagind(A::LinearMap, k::Integer=0)
+    require_one_based_indexing(A)
+    diagind(size(A,1), size(A,2), k)
+end
+
+LinearAlgebra.diag(A::LinearMap, k::Integer=0) = A[diagind(A,k)]
+
+# logical indexing
+Base.getindex(A::LinearMap, mask::AbstractVecOrMat{Bool}) = A[findall(mask)]
+Base.getindex(A::LinearMap, i, mask::AbstractVector{Bool}) = A[i, findall(mask)]
+Base.getindex(A::LinearMap, mask::AbstractVector{Bool}, j) = A[findall(mask), j]
+Base.getindex(A::LinearMap, im::AbstractVector{Bool}, jm::AbstractVector{Bool}) =
+    A[findall(im), findall(jm)]
+# disambiguation
+for typ in (:WrappedMap, :ScaledMap)
+    @eval begin
+        Base.getindex(A::$typ, mask::AbstractVecOrMat{Bool}) = A[findall(mask)]
+        Base.getindex(A::$typ, i, mask::AbstractVector{Bool}) = A[i, findall(mask)]
+        Base.getindex(A::$typ, mask::AbstractVector{Bool}, j) = A[findall(mask), j]
+        Base.getindex(A::$typ, im::AbstractVector{Bool}, jm::AbstractVector{Bool}) =
+            A[findall(im), findall(jm)]       
+    end
+end
+
 # nogetindex_error() = error("indexing not allowed for LinearMaps; consider setting `LinearMaps.allowgetindex = true`")
 
 # end # module
diff --git a/test/getindex.jl b/test/getindex.jl
index 80b0572e..eca93682 100644
--- a/test/getindex.jl
+++ b/test/getindex.jl
@@ -4,6 +4,7 @@ using LinearAlgebra, LinearMaps, Test
 
 function test_getindex(A::LinearMap, M::AbstractMatrix)
     @assert size(A) == size(M)
+    mask = rand(Bool, size(A))
     @test all((A[i,j] == M[i,j] for i in axes(A, 1), j in axes(A, 2)))
     @test all((A[i] == M[i] for i in 1:length(A)))
     @test A[1,1] == M[1,1]
@@ -20,6 +21,10 @@ function test_getindex(A::LinearMap, M::AbstractMatrix)
     @test A[:,:] == M
     @test A[7] == M[7]
     @test A[3:7] == M[3:7]
+    @test (lastindex(A, 1), lastindex(A, 2)) == size(A)
+    @test diagind(A) == diagind(M)
+    for k in -1:1; @test diag(A, k) == diag(M, k) end
+    @test A[mask] == M[mask]
     @test_throws BoundsError A[firstindex(A)-1]
     @test_throws BoundsError A[lastindex(A)+1]
     @test_throws BoundsError A[6,1]
@@ -30,9 +35,9 @@ function test_getindex(A::LinearMap, M::AbstractMatrix)
 end
 
 @testset "getindex" begin
-    A = rand(4,6)
-    L = LinearMap(A)
-    @test test_getindex(L, A)
+    M = rand(4,6)
+    A = LinearMap(M)
+    @test test_getindex(A, M)
     # @btime getindex($A, i) setup=(i = rand(1:24));
     # @btime getindex($L, i) setup=(i = rand(1:24));
     # @btime (getindex($A, i, j)) setup=(i = rand(1:4); j = rand(1:6));

From 6582b3e6bf95cce6fcb7f7239d6b7b1f3284ba51 Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <Daniel.Karrasch@posteo.de>
Date: Mon, 23 May 2022 19:56:06 +0200
Subject: [PATCH 10/23] revert README change

---
 README.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/README.md b/README.md
index 856c6ea1..3e4a2008 100644
--- a/README.md
+++ b/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://raw.githubusercontent.com/JuliaTesting/Aqua.jl/master/badge.svg
+[aqua-img]: https://img.shields.io/badge/Aqua.jl-%F0%9F%8C%A2-aqua.svg
 [aqua-url]: https://github.com/JuliaTesting/Aqua.jl

From 04b86fc25a33d87aa763ffef7df9285d63544820 Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <Daniel.Karrasch@posteo.de>
Date: Mon, 23 May 2022 19:59:00 +0200
Subject: [PATCH 11/23] revert unrelated test changes

---
 test/linearmaps.jl | 64 +++++++++++++++++++++++-----------------------
 1 file changed, 32 insertions(+), 32 deletions(-)

diff --git a/test/linearmaps.jl b/test/linearmaps.jl
index 22c4efbf..2783cfe6 100644
--- a/test/linearmaps.jl
+++ b/test/linearmaps.jl
@@ -33,20 +33,20 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays
     end
 end
 
-@testset "new LinearMap type" begin
-    # new type
-    struct SimpleFunctionMap <: LinearMap{Float64}
-        f::Function
-        N::Int
-    end
-    struct SimpleComplexFunctionMap <: LinearMap{Complex{Float64}}
-        f::Function
-        N::Int
-    end
-    Base.size(A::Union{SimpleFunctionMap,SimpleComplexFunctionMap}) = (A.N, A.N)
-    Base.:(*)(A::Union{SimpleFunctionMap,SimpleComplexFunctionMap}, v::AbstractVector) = A.f(v)
-    LinearAlgebra.mul!(y::AbstractVector, A::Union{SimpleFunctionMap,SimpleComplexFunctionMap}, x::AbstractVector) = copyto!(y, *(A, x))
+# new type
+struct SimpleFunctionMap <: LinearMap{Float64}
+    f::Function
+    N::Int
+end
+struct SimpleComplexFunctionMap <: LinearMap{Complex{Float64}}
+    f::Function
+    N::Int
+end
+Base.size(A::Union{SimpleFunctionMap,SimpleComplexFunctionMap}) = (A.N, A.N)
+Base.:(*)(A::Union{SimpleFunctionMap,SimpleComplexFunctionMap}, v::AbstractVector) = A.f(v)
+LinearAlgebra.mul!(y::AbstractVector, A::Union{SimpleFunctionMap,SimpleComplexFunctionMap}, x::AbstractVector) = copyto!(y, *(A, x))
 
+@testset "new LinearMap type" begin
     F = SimpleFunctionMap(cumsum, 10)
     FC = SimpleComplexFunctionMap(cumsum, 10)
     @test @inferred ndims(F) == 2
@@ -84,27 +84,27 @@ end
     @test Fs isa SparseMatrixCSC
 end
 
-@testset "transpose of new LinearMap type" begin
-    struct MyFillMap{T} <: LinearMaps.LinearMap{T}
-        λ::T
-        size::Dims{2}
-        function MyFillMap(λ::T, dims::Dims{2}) where {T}
-            all(d -> d >= 0, dims) || throw(ArgumentError("dims of MyFillMap must be non-negative"))
-            promote_type(T, typeof(λ)) == T || throw(InexactError())
-            return new{T}(λ, dims)
-        end
-    end
-    Base.size(A::MyFillMap) = A.size
-    function LinearAlgebra.mul!(y::AbstractVecOrMat, A::MyFillMap, x::AbstractVector)
-        LinearMaps.check_dim_mul(y, A, x)
-        return fill!(y, iszero(A.λ) ? zero(eltype(y)) : A.λ*sum(x))
-    end
-    function LinearAlgebra.mul!(y::AbstractVecOrMat, transA::LinearMaps.TransposeMap{<:Any,<:MyFillMap}, x::AbstractVector)
-        LinearMaps.check_dim_mul(y, transA, x)
-        λ = transA.lmap.λ
-        return fill!(y, iszero(λ) ? zero(eltype(y)) : transpose(λ)*sum(x))
+struct MyFillMap{T} <: LinearMaps.LinearMap{T}
+    λ::T
+    size::Dims{2}
+    function MyFillMap(λ::T, dims::Dims{2}) where {T}
+        all(d -> d >= 0, dims) || throw(ArgumentError("dims of MyFillMap must be non-negative"))
+        promote_type(T, typeof(λ)) == T || throw(InexactError())
+        return new{T}(λ, dims)
     end
+end
+Base.size(A::MyFillMap) = A.size
+function LinearAlgebra.mul!(y::AbstractVecOrMat, A::MyFillMap, x::AbstractVector)
+    LinearMaps.check_dim_mul(y, A, x)
+    return fill!(y, iszero(A.λ) ? zero(eltype(y)) : A.λ*sum(x))
+end
+function LinearAlgebra.mul!(y::AbstractVecOrMat, transA::LinearMaps.TransposeMap{<:Any,<:MyFillMap}, x::AbstractVector)
+    LinearMaps.check_dim_mul(y, transA, x)
+    λ = transA.lmap.λ
+    return fill!(y, iszero(λ) ? zero(eltype(y)) : transpose(λ)*sum(x))
+end
 
+@testset "transpose of new LinearMap type" begin
     A = MyFillMap(5.0, (3, 3))
     x = ones(3)
     @test A * x == fill(15.0, 3)

From a205d15d7050c12bda587123773225c2ec5b4de0 Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <Daniel.Karrasch@posteo.de>
Date: Mon, 23 May 2022 22:16:40 +0200
Subject: [PATCH 12/23] add history entry

---
 docs/src/history.md | 13 +++++++++++++
 1 file changed, 13 insertions(+)

diff --git a/docs/src/history.md b/docs/src/history.md
index 27909a3e..ca68b6fc 100644
--- a/docs/src/history.md
+++ b/docs/src/history.md
@@ -24,6 +24,19 @@
   * `LinearMap(A::MapOrVecOrMat, dims::Dims{2}; offset::Dims{2})`, where the keyword
     argument `offset` determines the dimension of a virtual upper-left zero block, to which
     `A` gets (virtually) diagonally appended.
+* An often requested new feature has been added: indexing/slicing any `LinearMap` object via
+  `Base.getindex` overloads. Linear, Cartesian (both with two `Integer` indices and
+  (vectors of) `CartesianIndex{2}`), colon, and logical indexing as well as slicing all
+  work generically on any `LinearMap` subtype. For performance reasons, some slicing
+  operations require the adjoint operation of the `LinearMap` type to be defined, or will
+  fail otherwise. Important note: `LinearMap` objects are meant to model objects that act
+  on vectors efficiently, and are in general *not* backed up by storage-like types like
+  `Array`s. Therefore, indexing and/or slicing of `LinearMap`s is likely to be slow,
+  and it may require the (repeated) allocation of canonical unit vectors. As a consequence,
+  generic algorithms relying heavily on indexing and/or slicing are likely to work much
+  worse than expected for `AbstractArray`s. To avoid repeated indexing operations which
+  may involve redundant computations, it is strongly recommended to consider `convert`ing
+  `LinearMap`-typed objects to `Matrix` or `SparseMatrixCSC` first.
 
 ## What's new in v3.6
 

From b9bdfe832bb354493c14e0c6b72f3447d70e2be0 Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <Daniel.Karrasch@posteo.de>
Date: Tue, 24 May 2022 09:57:24 +0200
Subject: [PATCH 13/23] make adjoint a requirement for horizontal slicing

---
 src/getindex.jl | 28 ++++++++++++++--------------
 1 file changed, 14 insertions(+), 14 deletions(-)

diff --git a/src/getindex.jl b/src/getindex.jl
index 97ff1e37..81af821f 100644
--- a/src/getindex.jl
+++ b/src/getindex.jl
@@ -54,30 +54,30 @@ _getindex(A::LinearMap, I::Vector{CartesianIndex{2}}) = [(@inbounds A[i]) for i
 _getindex(A::LinearMap, i::Union{Integer,Indexer}, j::Integer) = (@inbounds (A*basevec(A, 2, j))[i])
 _getindex(A::LinearMap, ::Base.Slice, j::Integer) = A*basevec(A, 2, j)
 function _getindex(A::LinearMap, i::Integer, J::Indexer)
-    try
+    # try
         # requires adjoint action to be defined
         return @inbounds (basevec(A, 1, i)'A)[J]
-    catch
-        return _fillbycols!(zeros(eltype(A), Base.index_shape(i, J)), A, i, J)
-    end
+    # catch
+    #     return _fillbycols!(zeros(eltype(A), Base.index_shape(i, J)), A, i, J)
+    # end
 end
 function _getindex(A::LinearMap, i::Integer, J::Base.Slice)
-    try
+    # try
         # requires adjoint action to be defined
         return vec(basevec(A, 1, i)'A)
-    catch
-        return _fillbycols!(zeros(eltype(A), Base.index_shape(i, J)), A, i, J)
-    end
+    # catch
+    #     return _fillbycols!(zeros(eltype(A), Base.index_shape(i, J)), A, i, J)
+    # end
 end
 function _getindex(A::LinearMap, I::Indexer, J::Indexer)
     dest = zeros(eltype(A), Base.index_shape(I, J))
     if length(I) <= length(J)
-        try
+        # try
             # requires adjoint action to be defined
             _fillbyrows!(dest, A, I, J)
-        catch
-            _fillbycols!(dest, A, I, J)
-        end
+        # catch
+        #     _fillbycols!(dest, A, I, J)
+        # end
     else
         _fillbycols!(dest, A, I, J)
     end
@@ -103,10 +103,10 @@ end
 function _fillbyrows!(dest, A, I, J)
     x = zeros(eltype(A), size(A, 1))
     temp = similar(x, eltype(A), size(A, 2))
-    @views @inbounds for (ind, i) in enumerate(I)
+    @views @inbounds for (di, i) in zip(eachcol(dest), I)
         x[i] = one(eltype(A))
         _unsafe_mul!(temp, A', x)
-        dest[ind,:] .= adjoint.(temp[J])
+        di .= adjoint.(temp[J])
         x[i] = zero(eltype(A))
     end
     return dest

From 896d59c4031bd09c4ed453bfeabd7b5f1d97a491 Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <Daniel.Karrasch@posteo.de>
Date: Tue, 24 May 2022 11:34:58 +0200
Subject: [PATCH 14/23] fix logical indexing

---
 src/getindex.jl  | 23 ++++-------------------
 test/getindex.jl | 37 +++++++++++++++++++------------------
 2 files changed, 23 insertions(+), 37 deletions(-)

diff --git a/src/getindex.jl b/src/getindex.jl
index 81af821f..79ef0c09 100644
--- a/src/getindex.jl
+++ b/src/getindex.jl
@@ -23,7 +23,8 @@ end
 # checkbounds in indexing via CartesianIndex
 Base.checkbounds(A::LinearMap, i::Union{CartesianIndex{2}, AbstractVecOrMat{CartesianIndex{2}}}) =
     Base.checkbounds_indices(Bool, axes(A), (i,))
-Base.checkbounds(A::LinearMap, I::AbstractArray{Bool,2}) = axes(A) == axes(I)
+Base.checkbounds(A::LinearMap, I::AbstractMatrix{Bool}) =
+    axes(A) == axes(I) || Base.throw_boundserror(A, I)
 
 # main entry point
 function Base.getindex(A::LinearMap, I...)
@@ -47,6 +48,7 @@ end
 _getindex(A::LinearMap, I::Indexer) = [_getindex(A, i) for i in I]
 _getindex(A::LinearMap, ::Base.Slice) = vec(Matrix(A))
 _getindex(A::LinearMap, I::Vector{CartesianIndex{2}}) = [(@inbounds A[i]) for i in I]
+_getindex(A::LinearMap, I::Base.LogicalIndex) = [(@inbounds A[i]) for i in I]
 
 ########################
 # Cartesian indexing
@@ -103,7 +105,7 @@ end
 function _fillbyrows!(dest, A, I, J)
     x = zeros(eltype(A), size(A, 1))
     temp = similar(x, eltype(A), size(A, 2))
-    @views @inbounds for (di, i) in zip(eachcol(dest), I)
+    @views @inbounds for (di, i) in zip(eachrow(dest), I)
         x[i] = one(eltype(A))
         _unsafe_mul!(temp, A', x)
         di .= adjoint.(temp[J])
@@ -144,23 +146,6 @@ end
 
 LinearAlgebra.diag(A::LinearMap, k::Integer=0) = A[diagind(A,k)]
 
-# logical indexing
-Base.getindex(A::LinearMap, mask::AbstractVecOrMat{Bool}) = A[findall(mask)]
-Base.getindex(A::LinearMap, i, mask::AbstractVector{Bool}) = A[i, findall(mask)]
-Base.getindex(A::LinearMap, mask::AbstractVector{Bool}, j) = A[findall(mask), j]
-Base.getindex(A::LinearMap, im::AbstractVector{Bool}, jm::AbstractVector{Bool}) =
-    A[findall(im), findall(jm)]
-# disambiguation
-for typ in (:WrappedMap, :ScaledMap)
-    @eval begin
-        Base.getindex(A::$typ, mask::AbstractVecOrMat{Bool}) = A[findall(mask)]
-        Base.getindex(A::$typ, i, mask::AbstractVector{Bool}) = A[i, findall(mask)]
-        Base.getindex(A::$typ, mask::AbstractVector{Bool}, j) = A[findall(mask), j]
-        Base.getindex(A::$typ, im::AbstractVector{Bool}, jm::AbstractVector{Bool}) =
-            A[findall(im), findall(jm)]       
-    end
-end
-
 # nogetindex_error() = error("indexing not allowed for LinearMaps; consider setting `LinearMaps.allowgetindex = true`")
 
 # end # module
diff --git a/test/getindex.jl b/test/getindex.jl
index eca93682..837f8d7a 100644
--- a/test/getindex.jl
+++ b/test/getindex.jl
@@ -5,6 +5,8 @@ using LinearAlgebra, LinearMaps, Test
 function test_getindex(A::LinearMap, M::AbstractMatrix)
     @assert size(A) == size(M)
     mask = rand(Bool, size(A))
+    imask = rand(Bool, size(A, 1))
+    jmask = rand(Bool, size(A, 2))
     @test all((A[i,j] == M[i,j] for i in axes(A, 1), j in axes(A, 2)))
     @test all((A[i] == M[i] for i in 1:length(A)))
     @test A[1,1] == M[1,1]
@@ -25,12 +27,18 @@ function test_getindex(A::LinearMap, M::AbstractMatrix)
     @test diagind(A) == diagind(M)
     for k in -1:1; @test diag(A, k) == diag(M, k) end
     @test A[mask] == M[mask]
+    @test A[imask, 1] == M[imask, 1]
+    @test A[1, jmask] == M[1, jmask]
+    @test A[imask, jmask] == M[imask, jmask]
     @test_throws BoundsError A[firstindex(A)-1]
     @test_throws BoundsError A[lastindex(A)+1]
     @test_throws BoundsError A[6,1]
     @test_throws BoundsError A[1,7]
     @test_throws BoundsError A[2,1:7]
     @test_throws BoundsError A[1:6,2]
+    @test_throws BoundsError A[ones(Bool, 2, 2)]
+    @test_throws BoundsError A[[true, true], 1]
+    @test_throws BoundsError A[1, [true, true]]
     return true
 end
 
@@ -45,31 +53,24 @@ end
 
     struct TwoMap <: LinearMaps.LinearMap{Float64} end
     Base.size(::TwoMap) = (5,5)
+    Base.transpose(A::TwoMap) = A
     LinearMaps._getindex(::TwoMap, i::Integer, j::Integer) = 2.0
     LinearMaps._unsafe_mul!(y::AbstractVector, ::TwoMap, x::AbstractVector) = fill!(y, 2.0*sum(x))
 
     T = TwoMap()
     @test test_getindex(TwoMap(), fill(2.0, size(T)))
-    Base.adjoint(A::TwoMap) = A
-    @test test_getindex(TwoMap(), fill(2.0, size(T)))
 
     MA = rand(ComplexF64, 5, 5)
-    for FA in (
-        LinearMap{ComplexF64}((y, x) -> mul!(y, MA, x), (y, x) -> mul!(y, MA', x), 5, 5),
-        LinearMap{ComplexF64}((y, x) -> mul!(y, MA, x), 5, 5),
-    )
-        @test test_getindex(FA, MA)
-        @test test_getindex(3FA, 3MA)
-        @test test_getindex(FA + FA, 2MA)
-        if !isnothing(FA.fc)
-            @test test_getindex(transpose(FA), transpose(MA))
-            @test test_getindex(transpose(3FA), transpose(3MA))
-            @test test_getindex(3transpose(FA), transpose(3MA))
-            @test test_getindex(adjoint(FA), adjoint(MA))
-            @test test_getindex(adjoint(3FA), adjoint(3MA))
-            @test test_getindex(3adjoint(FA), adjoint(3MA))
-        end
-    end
+    FA = LinearMap{ComplexF64}((y, x) -> mul!(y, MA, x), (y, x) -> mul!(y, MA', x), 5, 5)
+    @test test_getindex(FA, MA)
+    @test test_getindex(3FA, 3MA)
+    @test test_getindex(FA + FA, 2MA)
+    @test test_getindex(transpose(FA), transpose(MA))
+    @test test_getindex(transpose(3FA), transpose(3MA))
+    @test test_getindex(3transpose(FA), transpose(3MA))
+    @test test_getindex(adjoint(FA), adjoint(MA))
+    @test test_getindex(adjoint(3FA), adjoint(3MA))
+    @test test_getindex(3adjoint(FA), adjoint(3MA))
 
     @test test_getindex(FillMap(0.5, (5, 5)), fill(0.5, (5, 5)))
     @test test_getindex(LinearMap(0.5I, 5), Matrix(0.5I, 5, 5))

From c6be0bca1bc9fc5b3f6a655887908cfaebb45cec Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <Daniel.Karrasch@posteo.de>
Date: Tue, 24 May 2022 21:32:40 +0200
Subject: [PATCH 15/23] minor fixes, more tests

---
 src/getindex.jl  | 16 ++++------------
 test/getindex.jl | 29 ++++++++++++++++-------------
 2 files changed, 20 insertions(+), 25 deletions(-)

diff --git a/src/getindex.jl b/src/getindex.jl
index 79ef0c09..ab13006d 100644
--- a/src/getindex.jl
+++ b/src/getindex.jl
@@ -1,8 +1,3 @@
-# module GetIndex
-
-# using ..LinearMaps: LinearMap, AdjointMap, TransposeMap, FillMap, LinearCombination,
-#     ScaledMap, UniformScalingMap, WrappedMap
-
 const Indexer = AbstractVector{<:Integer}
 
 Base.IndexStyle(::LinearMap) = IndexCartesian()
@@ -22,9 +17,9 @@ function Base.checkbounds(A::LinearMap, i)
 end
 # checkbounds in indexing via CartesianIndex
 Base.checkbounds(A::LinearMap, i::Union{CartesianIndex{2}, AbstractVecOrMat{CartesianIndex{2}}}) =
-    Base.checkbounds_indices(Bool, axes(A), (i,))
+    Base.checkbounds_indices(Bool, axes(A), (i,)) || throw(BoundsError(A, i))
 Base.checkbounds(A::LinearMap, I::AbstractMatrix{Bool}) =
-    axes(A) == axes(I) || Base.throw_boundserror(A, I)
+    axes(A) == axes(I) || throw(BoundsError(A, I))
 
 # main entry point
 function Base.getindex(A::LinearMap, I...)
@@ -134,7 +129,8 @@ function _fillbycols!(dest, A, ::Base.Slice, J)
     return dest
 end
 
-@inline _copycol!(dest, ind, temp, i::Integer) = (@inbounds dest[ind] = temp[i])
+# needed only if we accept the try-catch blocks above
+# @inline _copycol!(dest, ind, temp, i::Integer) = (@inbounds dest[ind] = temp[i])
 @inline _copycol!(dest, ind, temp, I::Indexer) =
     (@views @inbounds dest[:,ind] .= temp[I])
 
@@ -145,7 +141,3 @@ function LinearAlgebra.diagind(A::LinearMap, k::Integer=0)
 end
 
 LinearAlgebra.diag(A::LinearMap, k::Integer=0) = A[diagind(A,k)]
-
-# nogetindex_error() = error("indexing not allowed for LinearMaps; consider setting `LinearMaps.allowgetindex = true`")
-
-# end # module
diff --git a/test/getindex.jl b/test/getindex.jl
index 837f8d7a..b23f03e0 100644
--- a/test/getindex.jl
+++ b/test/getindex.jl
@@ -27,6 +27,8 @@ function test_getindex(A::LinearMap, M::AbstractMatrix)
     @test diagind(A) == diagind(M)
     for k in -1:1; @test diag(A, k) == diag(M, k) end
     @test A[mask] == M[mask]
+    @test A[findall(mask)] == M[findall(mask)]
+    @test A[CartesianIndex(1,1)] == M[CartesianIndex(1,1)]
     @test A[imask, 1] == M[imask, 1]
     @test A[1, jmask] == M[1, jmask]
     @test A[imask, jmask] == M[imask, jmask]
@@ -39,13 +41,14 @@ function test_getindex(A::LinearMap, M::AbstractMatrix)
     @test_throws BoundsError A[ones(Bool, 2, 2)]
     @test_throws BoundsError A[[true, true], 1]
     @test_throws BoundsError A[1, [true, true]]
+    @test_throws BoundsError A[CartesianIndex(0,1)]
     return true
 end
 
 @testset "getindex" begin
     M = rand(4,6)
     A = LinearMap(M)
-    @test test_getindex(A, M)
+    test_getindex(A, M)
     # @btime getindex($A, i) setup=(i = rand(1:24));
     # @btime getindex($L, i) setup=(i = rand(1:24));
     # @btime (getindex($A, i, j)) setup=(i = rand(1:4); j = rand(1:6));
@@ -58,20 +61,20 @@ end
     LinearMaps._unsafe_mul!(y::AbstractVector, ::TwoMap, x::AbstractVector) = fill!(y, 2.0*sum(x))
 
     T = TwoMap()
-    @test test_getindex(TwoMap(), fill(2.0, size(T)))
+    test_getindex(TwoMap(), fill(2.0, size(T)))
 
     MA = rand(ComplexF64, 5, 5)
     FA = LinearMap{ComplexF64}((y, x) -> mul!(y, MA, x), (y, x) -> mul!(y, MA', x), 5, 5)
-    @test test_getindex(FA, MA)
-    @test test_getindex(3FA, 3MA)
-    @test test_getindex(FA + FA, 2MA)
-    @test test_getindex(transpose(FA), transpose(MA))
-    @test test_getindex(transpose(3FA), transpose(3MA))
-    @test test_getindex(3transpose(FA), transpose(3MA))
-    @test test_getindex(adjoint(FA), adjoint(MA))
-    @test test_getindex(adjoint(3FA), adjoint(3MA))
-    @test test_getindex(3adjoint(FA), adjoint(3MA))
+    test_getindex(FA, MA)
+    test_getindex(3FA, 3MA)
+    test_getindex(FA + FA, 2MA)
+    test_getindex(transpose(FA), transpose(MA))
+    test_getindex(transpose(3FA), transpose(3MA))
+    test_getindex(3transpose(FA), transpose(3MA))
+    test_getindex(adjoint(FA), adjoint(MA))
+    test_getindex(adjoint(3FA), adjoint(3MA))
+    test_getindex(3adjoint(FA), adjoint(3MA))
 
-    @test test_getindex(FillMap(0.5, (5, 5)), fill(0.5, (5, 5)))
-    @test test_getindex(LinearMap(0.5I, 5), Matrix(0.5I, 5, 5))
+    test_getindex(FillMap(0.5, (5, 5)), fill(0.5, (5, 5)))
+    test_getindex(LinearMap(0.5I, 5), Matrix(0.5I, 5, 5))
 end

From 3e0f423f0636ec9753fea43146077b87440233df Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <Daniel.Karrasch@posteo.de>
Date: Fri, 3 Jun 2022 17:04:50 +0200
Subject: [PATCH 16/23] rename basevec -> unitvec

---
 src/getindex.jl | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/src/getindex.jl b/src/getindex.jl
index ab13006d..9e952fab 100644
--- a/src/getindex.jl
+++ b/src/getindex.jl
@@ -48,12 +48,12 @@ _getindex(A::LinearMap, I::Base.LogicalIndex) = [(@inbounds A[i]) for i in I]
 ########################
 # Cartesian indexing
 ########################
-_getindex(A::LinearMap, i::Union{Integer,Indexer}, j::Integer) = (@inbounds (A*basevec(A, 2, j))[i])
-_getindex(A::LinearMap, ::Base.Slice, j::Integer) = A*basevec(A, 2, j)
+_getindex(A::LinearMap, i::Union{Integer,Indexer}, j::Integer) = (@inbounds (A*unitvec(A, 2, j))[i])
+_getindex(A::LinearMap, ::Base.Slice, j::Integer) = A*unitvec(A, 2, j)
 function _getindex(A::LinearMap, i::Integer, J::Indexer)
     # try
         # requires adjoint action to be defined
-        return @inbounds (basevec(A, 1, i)'A)[J]
+        return @inbounds (unitvec(A, 1, i)'A)[J]
     # catch
     #     return _fillbycols!(zeros(eltype(A), Base.index_shape(i, J)), A, i, J)
     # end
@@ -61,7 +61,7 @@ end
 function _getindex(A::LinearMap, i::Integer, J::Base.Slice)
     # try
         # requires adjoint action to be defined
-        return vec(basevec(A, 1, i)'A)
+        return vec(unitvec(A, 1, i)'A)
     # catch
     #     return _fillbycols!(zeros(eltype(A), Base.index_shape(i, J)), A, i, J)
     # end
@@ -91,7 +91,7 @@ _getindex(A::TransposeMap, i::Integer, j::Integer) = @inbounds transpose(A.lmap[
 _getindex(A::UniformScalingMap, i::Integer, j::Integer) = ifelse(i == j, A.λ, zero(eltype(A)))
 
 # helpers
-function basevec(A, dim, i::Integer)
+function unitvec(A, dim, i::Integer)
     x = zeros(eltype(A), size(A, dim))
     @inbounds x[i] = one(eltype(A))
     return x

From 3fbd708ee87eb50a0a26fc4d8d0da752bf981396 Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <daniel.karrasch@posteo.de>
Date: Sun, 5 Jun 2022 18:02:08 +0200
Subject: [PATCH 17/23] Apply suggestions from code review

Co-authored-by: Jeff Fessler <JeffFessler@users.noreply.github.com>
---
 docs/src/history.md | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/docs/src/history.md b/docs/src/history.md
index ca68b6fc..33b4af96 100644
--- a/docs/src/history.md
+++ b/docs/src/history.md
@@ -33,10 +33,10 @@
   on vectors efficiently, and are in general *not* backed up by storage-like types like
   `Array`s. Therefore, indexing and/or slicing of `LinearMap`s is likely to be slow,
   and it may require the (repeated) allocation of canonical unit vectors. As a consequence,
-  generic algorithms relying heavily on indexing and/or slicing are likely to work much
-  worse than expected for `AbstractArray`s. To avoid repeated indexing operations which
+  generic algorithms relying heavily on indexing and/or slicing are likely to run much
+  slower than expected for `AbstractArray`s. To avoid repeated indexing operations which
   may involve redundant computations, it is strongly recommended to consider `convert`ing
-  `LinearMap`-typed objects to `Matrix` or `SparseMatrixCSC` first.
+  `LinearMap`-typed objects to `Matrix` or `SparseMatrixCSC` first, if memory permits.
 
 ## What's new in v3.6
 

From 84a95bbc10f4026231b31b003ab2b2ade713d1b9 Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <Daniel.Karrasch@posteo.de>
Date: Thu, 9 Jun 2022 12:46:50 +0200
Subject: [PATCH 18/23] overhaul after code review

---
 docs/src/history.md | 39 ++++++++++++++---------
 src/getindex.jl     | 76 ++++++++++++++-------------------------------
 test/getindex.jl    | 20 +++---------
 3 files changed, 53 insertions(+), 82 deletions(-)

diff --git a/docs/src/history.md b/docs/src/history.md
index 33b4af96..0d32dcb8 100644
--- a/docs/src/history.md
+++ b/docs/src/history.md
@@ -1,6 +1,7 @@
 # Version history
 
 ## What's new in v3.7
+
 * `mul!(M::AbstractMatrix, A::LinearMap, s::Number, a, b)` methods are provided, mimicking
   similar methods in `Base.LinearAlgebra`. This version allows for the memory efficient
   implementation of in-place addition and conversion of a `LinearMap` to `Matrix`.
@@ -9,7 +10,6 @@
   conversion, construction, and inplace addition will benefit from this supplied efficient
   implementation. If no specialisation is supplied, a generic fallback is used that is based
   on feeding the canonical basis of unit vectors to the linear map.
-
 * A new map type called `EmbeddedMap` is introduced. It is a wrapper of a "small" `LinearMap`
   (or a suitably converted `AbstractVecOrMat`) embedded into a "larger" zero map. Hence,
   the "small" map acts only on a subset of the coordinates and maps to another subset of
@@ -24,19 +24,30 @@
   * `LinearMap(A::MapOrVecOrMat, dims::Dims{2}; offset::Dims{2})`, where the keyword
     argument `offset` determines the dimension of a virtual upper-left zero block, to which
     `A` gets (virtually) diagonally appended.
-* An often requested new feature has been added: indexing/slicing any `LinearMap` object via
-  `Base.getindex` overloads. Linear, Cartesian (both with two `Integer` indices and
-  (vectors of) `CartesianIndex{2}`), colon, and logical indexing as well as slicing all
-  work generically on any `LinearMap` subtype. For performance reasons, some slicing
-  operations require the adjoint operation of the `LinearMap` type to be defined, or will
-  fail otherwise. Important note: `LinearMap` objects are meant to model objects that act
-  on vectors efficiently, and are in general *not* backed up by storage-like types like
-  `Array`s. Therefore, indexing and/or slicing of `LinearMap`s is likely to be slow,
-  and it may require the (repeated) allocation of canonical unit vectors. As a consequence,
-  generic algorithms relying heavily on indexing and/or slicing are likely to run much
-  slower than expected for `AbstractArray`s. To avoid repeated indexing operations which
-  may involve redundant computations, it is strongly recommended to consider `convert`ing
-  `LinearMap`-typed objects to `Matrix` or `SparseMatrixCSC` first, if memory permits.
+* An often requested new feature has been added: slicing (i.e., non-scalar indexing) any
+  `LinearMap` object via `Base.getindex` overloads. Note, however, that only rather
+  efficient slicing operations are implemented: `A[:,j]`, `A[:,J]`, `A[I,J]`, and `A[:,:]`,
+  where `j::Integer` and `J` is either of type `AbstractVector{<:Integer>}` or an
+  `AbstractVector{Bool}` of appropriate length ("logical slicing").
+
+  Scalar indexing `A[i::Integer,j::Integer]` as well as other indexing operations that fall
+  back on scalar indexing such as logical indexing by some `AbstractMatrix{Bool}`, or
+  indexing by vectors of (linear or Cartesian) indices are not supported; as an exception,
+  `getindex` calls on wrapped `AbstractVecOrMat`s is forwarded to corresponding `getindex`
+  methods from `Base` and therefore allow any type of usual indexing/slicing.
+  If scalar indexing is really required, consider using `A[:,j][i]` which is as efficient
+  as a reasonable generic implementation for `LinearMap`s can be.
+
+  Furthermore, (predominantly) horizontal slicing operations require the adjoint operation
+  of the `LinearMap` type to be defined, or will fail otherwise. Important note:
+  `LinearMap` objects are meant to model objects that act on vectors efficiently, and are
+  in general *not* backed up by storage-like types like `Array`s. Therefore, slicing of
+  `LinearMap`s is potentially slow, and it may require the (repeated) allocation of
+  standard unit vectors. As a consequence, generic algorithms relying heavily on indexing
+  and/or slicing are likely to run much slower than expected for `AbstractArray`s. To avoid
+  repeated indexing operations which may involve redundant computations, it is strongly
+  recommended to consider `convert`ing `LinearMap`-typed objects to `Matrix` or
+  `SparseMatrixCSC` first, if memory permits.
 
 ## What's new in v3.6
 
diff --git a/src/getindex.jl b/src/getindex.jl
index 9e952fab..714c59e8 100644
--- a/src/getindex.jl
+++ b/src/getindex.jl
@@ -10,71 +10,64 @@ function Base.checkbounds(A::LinearMap, i, j)
     Base.checkbounds_indices(Bool, axes(A), (i, j)) || throw(BoundsError(A, (i, j)))
     nothing
 end
-# Linear indexing is explicitly allowed when there is only one (non-cartesian) index
 function Base.checkbounds(A::LinearMap, i)
     Base.checkindex(Bool, Base.OneTo(length(A)), i) || throw(BoundsError(A, i))
     nothing
 end
 # checkbounds in indexing via CartesianIndex
-Base.checkbounds(A::LinearMap, i::Union{CartesianIndex{2}, AbstractVecOrMat{CartesianIndex{2}}}) =
+Base.checkbounds(A::LinearMap, i::Union{CartesianIndex{2}, AbstractArray{CartesianIndex{2}}}) =
     Base.checkbounds_indices(Bool, axes(A), (i,)) || throw(BoundsError(A, i))
 Base.checkbounds(A::LinearMap, I::AbstractMatrix{Bool}) =
     axes(A) == axes(I) || throw(BoundsError(A, I))
 
 # main entry point
 function Base.getindex(A::LinearMap, I...)
-    # TODO: introduce some sort of switch?
     @boundscheck checkbounds(A, I...)
     _getindex(A, Base.to_indices(A, I)...)
 end
 # quick pass forward
-Base.@propagate_inbounds Base.getindex(A::ScaledMap, I...) = A.λ .* A.lmap[I...]
+Base.@propagate_inbounds Base.getindex(A::ScaledMap, I...) = A.λ * A.lmap[I...]
 Base.@propagate_inbounds Base.getindex(A::WrappedMap, I...) = A.lmap[I...]
-Base.@propagate_inbounds Base.getindex(A::WrappedMap, i::Integer) = A.lmap[i]
-Base.@propagate_inbounds Base.getindex(A::WrappedMap, i::Integer, j::Integer) = A.lmap[i,j]
 
 ########################
 # linear indexing
 ########################
-function _getindex(A::LinearMap, i::Integer)
-    i1, i2 = Base._ind2sub(axes(A), i)
-    return _getindex(A, i1, i2)
-end
-_getindex(A::LinearMap, I::Indexer) = [_getindex(A, i) for i in I]
+_getindex(A::LinearMap, _) = error("linear indexing of LinearMaps is not supported")
 _getindex(A::LinearMap, ::Base.Slice) = vec(Matrix(A))
-_getindex(A::LinearMap, I::Vector{CartesianIndex{2}}) = [(@inbounds A[i]) for i in I]
-_getindex(A::LinearMap, I::Base.LogicalIndex) = [(@inbounds A[i]) for i in I]
 
 ########################
 # Cartesian indexing
 ########################
-_getindex(A::LinearMap, i::Union{Integer,Indexer}, j::Integer) = (@inbounds (A*unitvec(A, 2, j))[i])
+_getindex(A::LinearMap, i::Integer, j::Integer) =
+    error("scalar indexing of LinearMaps is not supported, consider using A[:,j][i] instead")
+_getindex(A::LinearMap, I::Indexer, j::Integer) = (@inbounds (A*unitvec(A, 2, j))[I])
 _getindex(A::LinearMap, ::Base.Slice, j::Integer) = A*unitvec(A, 2, j)
 function _getindex(A::LinearMap, i::Integer, J::Indexer)
-    # try
+    try
         # requires adjoint action to be defined
         return @inbounds (unitvec(A, 1, i)'A)[J]
-    # catch
-    #     return _fillbycols!(zeros(eltype(A), Base.index_shape(i, J)), A, i, J)
-    # end
+    catch
+        error("efficient horizontal slicing A[$i,$J] requires the adjoint of $(typeof(A)) to be defined")
+    end
 end
 function _getindex(A::LinearMap, i::Integer, J::Base.Slice)
-    # try
+    try
         # requires adjoint action to be defined
         return vec(unitvec(A, 1, i)'A)
-    # catch
-    #     return _fillbycols!(zeros(eltype(A), Base.index_shape(i, J)), A, i, J)
-    # end
+    catch
+        error("efficient horizontal slicing A[$i,:] requires the adjoint of $(typeof(A)) to be defined")
+    end
 end
 function _getindex(A::LinearMap, I::Indexer, J::Indexer)
     dest = zeros(eltype(A), Base.index_shape(I, J))
+    # choose whichever requires less map applications
     if length(I) <= length(J)
-        # try
+        try
             # requires adjoint action to be defined
             _fillbyrows!(dest, A, I, J)
-        # catch
-        #     _fillbycols!(dest, A, I, J)
-        # end
+        catch
+            error("efficient horizontal slicing A[I,J] with length(I) <= length(J) requires the adjoint of $(typeof(A)) to be defined")
+        end
     else
         _fillbycols!(dest, A, I, J)
     end
@@ -82,16 +75,8 @@ function _getindex(A::LinearMap, I::Indexer, J::Indexer)
 end
 _getindex(A::LinearMap, ::Base.Slice, ::Base.Slice) = Matrix(A)
 
-# specialized methods
-_getindex(A::FillMap, ::Integer, ::Integer) = A.λ
-_getindex(A::LinearCombination, i::Integer, j::Integer) =
-    sum(a -> (@inbounds A.maps[a][i, j]), eachindex(A.maps))
-_getindex(A::AdjointMap, i::Integer, j::Integer) = @inbounds adjoint(A.lmap[j, i])
-_getindex(A::TransposeMap, i::Integer, j::Integer) = @inbounds transpose(A.lmap[j, i])
-_getindex(A::UniformScalingMap, i::Integer, j::Integer) = ifelse(i == j, A.λ, zero(eltype(A)))
-
 # helpers
-function unitvec(A, dim, i::Integer)
+function unitvec(A, dim, i)
     x = zeros(eltype(A), size(A, dim))
     @inbounds x[i] = one(eltype(A))
     return x
@@ -108,36 +93,23 @@ function _fillbyrows!(dest, A, I, J)
     end
     return dest
 end
-function _fillbycols!(dest, A, i, J)
+function _fillbycols!(dest, A, I::Indexer, J)
     x = zeros(eltype(A), size(A, 2))
     temp = similar(x, eltype(A), size(A, 1))
-    @views @inbounds for (ind, j) in enumerate(J)
+    @inbounds for (ind, j) in enumerate(J)
         x[j] = one(eltype(A))
         _unsafe_mul!(temp, A, x)
-        _copycol!(dest, ind, temp, i)
+        dest[:,ind] .= temp[I]
         x[j] = zero(eltype(A))
     end
     return dest
 end
 function _fillbycols!(dest, A, ::Base.Slice, J)
     x = zeros(eltype(A), size(A, 2))
-    @views @inbounds for (ind, j) in enumerate(J)
+    @inbounds for (ind, j) in enumerate(J)
         x[j] = one(eltype(A))
         _unsafe_mul!(selectdim(dest, 2, ind), A, x)
         x[j] = zero(eltype(A))
     end
     return dest
 end
-
-# needed only if we accept the try-catch blocks above
-# @inline _copycol!(dest, ind, temp, i::Integer) = (@inbounds dest[ind] = temp[i])
-@inline _copycol!(dest, ind, temp, I::Indexer) =
-    (@views @inbounds dest[:,ind] .= temp[I])
-
-# diagonal indexing
-function LinearAlgebra.diagind(A::LinearMap, k::Integer=0)
-    require_one_based_indexing(A)
-    diagind(size(A,1), size(A,2), k)
-end
-
-LinearAlgebra.diag(A::LinearMap, k::Integer=0) = A[diagind(A,k)]
diff --git a/test/getindex.jl b/test/getindex.jl
index b23f03e0..da5efccd 100644
--- a/test/getindex.jl
+++ b/test/getindex.jl
@@ -1,5 +1,4 @@
 using LinearAlgebra, LinearMaps, Test
-# using LinearMaps.GetIndex
 # using BenchmarkTools
 
 function test_getindex(A::LinearMap, M::AbstractMatrix)
@@ -7,9 +6,6 @@ function test_getindex(A::LinearMap, M::AbstractMatrix)
     mask = rand(Bool, size(A))
     imask = rand(Bool, size(A, 1))
     jmask = rand(Bool, size(A, 2))
-    @test all((A[i,j] == M[i,j] for i in axes(A, 1), j in axes(A, 2)))
-    @test all((A[i] == M[i] for i in 1:length(A)))
-    @test A[1,1] == M[1,1]
     @test A[:] == M[:]
     @test A[1,:] == M[1,:]
     @test A[:,1] == M[:,1]
@@ -21,19 +17,10 @@ function test_getindex(A::LinearMap, M::AbstractMatrix)
     @test A[1:2,1:3] == M[1:2,1:3]
     @test A[[2,1],1:3] == M[[2,1],1:3]
     @test A[:,:] == M
-    @test A[7] == M[7]
-    @test A[3:7] == M[3:7]
     @test (lastindex(A, 1), lastindex(A, 2)) == size(A)
-    @test diagind(A) == diagind(M)
-    for k in -1:1; @test diag(A, k) == diag(M, k) end
-    @test A[mask] == M[mask]
-    @test A[findall(mask)] == M[findall(mask)]
-    @test A[CartesianIndex(1,1)] == M[CartesianIndex(1,1)]
     @test A[imask, 1] == M[imask, 1]
     @test A[1, jmask] == M[1, jmask]
     @test A[imask, jmask] == M[imask, jmask]
-    @test_throws BoundsError A[firstindex(A)-1]
-    @test_throws BoundsError A[lastindex(A)+1]
     @test_throws BoundsError A[6,1]
     @test_throws BoundsError A[1,7]
     @test_throws BoundsError A[2,1:7]
@@ -41,7 +28,6 @@ function test_getindex(A::LinearMap, M::AbstractMatrix)
     @test_throws BoundsError A[ones(Bool, 2, 2)]
     @test_throws BoundsError A[[true, true], 1]
     @test_throws BoundsError A[1, [true, true]]
-    @test_throws BoundsError A[CartesianIndex(0,1)]
     return true
 end
 
@@ -49,10 +35,10 @@ end
     M = rand(4,6)
     A = LinearMap(M)
     test_getindex(A, M)
+    # @btime getindex($M, i) setup=(i = rand(1:24));
     # @btime getindex($A, i) setup=(i = rand(1:24));
-    # @btime getindex($L, i) setup=(i = rand(1:24));
+    # @btime (getindex($M, i, j)) setup=(i = rand(1:4); j = rand(1:6));
     # @btime (getindex($A, i, j)) setup=(i = rand(1:4); j = rand(1:6));
-    # @btime (getindex($L, i, j)) setup=(i = rand(1:4); j = rand(1:6));
 
     struct TwoMap <: LinearMaps.LinearMap{Float64} end
     Base.size(::TwoMap) = (5,5)
@@ -65,7 +51,9 @@ end
 
     MA = rand(ComplexF64, 5, 5)
     FA = LinearMap{ComplexF64}((y, x) -> mul!(y, MA, x), (y, x) -> mul!(y, MA', x), 5, 5)
+    F = LinearMap{ComplexF64}(x -> MA*x, y -> MA'y, 5, 5)
     test_getindex(FA, MA)
+    test_getindex(F, MA)
     test_getindex(3FA, 3MA)
     test_getindex(FA + FA, 2MA)
     test_getindex(transpose(FA), transpose(MA))

From 46865b889dd04dcb01de73ae4f4675a38c7698ff Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <Daniel.Karrasch@posteo.de>
Date: Thu, 9 Jun 2022 16:23:29 +0200
Subject: [PATCH 19/23] throw more

---
 src/getindex.jl  | 26 +++++++++++++-------------
 test/getindex.jl | 35 ++++++++++++++++++++++++-----------
 2 files changed, 37 insertions(+), 24 deletions(-)

diff --git a/src/getindex.jl b/src/getindex.jl
index 714c59e8..fc9e363f 100644
--- a/src/getindex.jl
+++ b/src/getindex.jl
@@ -40,25 +40,26 @@ _getindex(A::LinearMap, ::Base.Slice) = vec(Matrix(A))
 ########################
 _getindex(A::LinearMap, i::Integer, j::Integer) =
     error("scalar indexing of LinearMaps is not supported, consider using A[:,j][i] instead")
-_getindex(A::LinearMap, I::Indexer, j::Integer) = (@inbounds (A*unitvec(A, 2, j))[I])
+_getindex(A::LinearMap, I::Indexer, j::Integer) =
+    error("partial vertical slicing of LinearMaps is not supported, consider using A[:,j][I] instead")
 _getindex(A::LinearMap, ::Base.Slice, j::Integer) = A*unitvec(A, 2, j)
-function _getindex(A::LinearMap, i::Integer, J::Indexer)
-    try
-        # requires adjoint action to be defined
-        return @inbounds (unitvec(A, 1, i)'A)[J]
-    catch
-        error("efficient horizontal slicing A[$i,$J] requires the adjoint of $(typeof(A)) to be defined")
-    end
-end
+_getindex(A::LinearMap, i::Integer, J::Indexer) =
+    error("partial horizontal slicing of LinearMaps is not supported, consider using A[i,:][J] instead")
 function _getindex(A::LinearMap, i::Integer, J::Base.Slice)
     try
         # requires adjoint action to be defined
         return vec(unitvec(A, 1, i)'A)
     catch
-        error("efficient horizontal slicing A[$i,:] requires the adjoint of $(typeof(A)) to be defined")
+        error("horizontal slicing A[$i,:] requires the adjoint of $(typeof(A)) to be defined")
     end
 end
-function _getindex(A::LinearMap, I::Indexer, J::Indexer)
+_getindex(A::LinearMap, I::Indexer, J::Indexer) =
+    error("partial two-dimensional slicing of LinearMaps is not supported, consider using A[:,J][I] or A[I,:][J] instead")
+_getindex(A::LinearMap, ::Base.Slice, ::Base.Slice) =
+    error("two-dimensional slicing of LinearMaps is not supported, consider using Matrix(A) or convert(Matrix, A)") 
+_getindex(A::LinearMap, I::Base.Slice, J::Indexer) = __getindex(A, I, J)
+_getindex(A::LinearMap, I::Indexer, J::Base.Slice) = __getindex(A, I, J)
+function __getindex(A, I, J)
     dest = zeros(eltype(A), Base.index_shape(I, J))
     # choose whichever requires less map applications
     if length(I) <= length(J)
@@ -66,14 +67,13 @@ function _getindex(A::LinearMap, I::Indexer, J::Indexer)
             # requires adjoint action to be defined
             _fillbyrows!(dest, A, I, J)
         catch
-            error("efficient horizontal slicing A[I,J] with length(I) <= length(J) requires the adjoint of $(typeof(A)) to be defined")
+            error("wide slicing A[I,J] with length(I) <= length(J) requires the adjoint of $(typeof(A)) to be defined")
         end
     else
         _fillbycols!(dest, A, I, J)
     end
     return dest
 end
-_getindex(A::LinearMap, ::Base.Slice, ::Base.Slice) = Matrix(A)
 
 # helpers
 function unitvec(A, dim, i)
diff --git a/test/getindex.jl b/test/getindex.jl
index da5efccd..a0ff8e9a 100644
--- a/test/getindex.jl
+++ b/test/getindex.jl
@@ -1,9 +1,9 @@
 using LinearAlgebra, LinearMaps, Test
+using LinearMaps: VecOrMatMap, ScaledMap
 # using BenchmarkTools
 
 function test_getindex(A::LinearMap, M::AbstractMatrix)
     @assert size(A) == size(M)
-    mask = rand(Bool, size(A))
     imask = rand(Bool, size(A, 1))
     jmask = rand(Bool, size(A, 2))
     @test A[:] == M[:]
@@ -11,16 +11,28 @@ function test_getindex(A::LinearMap, M::AbstractMatrix)
     @test A[:,1] == M[:,1]
     @test A[1:4,:] == M[1:4,:]
     @test A[:,1:4] == M[:,1:4]
-    @test A[1,1:3] == M[1,1:3]
-    @test A[1:3,1] == M[1:3,1]
-    @test A[2:end,1] == M[2:end,1]
-    @test A[1:2,1:3] == M[1:2,1:3]
-    @test A[[2,1],1:3] == M[[2,1],1:3]
-    @test A[:,:] == M
     @test (lastindex(A, 1), lastindex(A, 2)) == size(A)
-    @test A[imask, 1] == M[imask, 1]
-    @test A[1, jmask] == M[1, jmask]
-    @test A[imask, jmask] == M[imask, jmask]
+    if A isa VecOrMatMap || A isa ScaledMap{<:Any,<:Any,<:VecOrMatMap}
+        @test A[1,1:3] == M[1,1:3]
+        @test A[1:3,1] == M[1:3,1]
+        @test A[2:end,1] == M[2:end,1]
+        @test A[1:2,1:3] == M[1:2,1:3]
+        @test A[[2,1],1:3] == M[[2,1],1:3]
+        @test A[:,:] == M
+        @test A[imask, 1] == M[imask, 1]
+        @test A[1, jmask] == M[1, jmask]
+        @test A[imask, jmask] == M[imask, jmask]
+    else
+        @test_throws ErrorException A[1,1:3] == M[1,1:3]
+        @test_throws ErrorException A[1:3,1] == M[1:3,1]
+        @test_throws ErrorException A[2:end,1] == M[2:end,1]
+        @test_throws ErrorException A[1:2,1:3] == M[1:2,1:3]
+        @test_throws ErrorException A[[2,1],1:3] == M[[2,1],1:3]
+        @test_throws ErrorException A[:,:] == M
+        @test_throws ErrorException A[imask, 1] == M[imask, 1]
+        @test_throws ErrorException A[1, jmask] == M[1, jmask]
+        @test_throws ErrorException A[imask, jmask] == M[imask, jmask]
+    end
     @test_throws BoundsError A[6,1]
     @test_throws BoundsError A[1,7]
     @test_throws BoundsError A[2,1:7]
@@ -28,13 +40,14 @@ function test_getindex(A::LinearMap, M::AbstractMatrix)
     @test_throws BoundsError A[ones(Bool, 2, 2)]
     @test_throws BoundsError A[[true, true], 1]
     @test_throws BoundsError A[1, [true, true]]
-    return true
+    return nothing
 end
 
 @testset "getindex" begin
     M = rand(4,6)
     A = LinearMap(M)
     test_getindex(A, M)
+    test_getindex(2A, 2M)
     # @btime getindex($M, i) setup=(i = rand(1:24));
     # @btime getindex($A, i) setup=(i = rand(1:24));
     # @btime (getindex($M, i, j)) setup=(i = rand(1:4); j = rand(1:6));

From 8f17ed0208d0a9410517da7c417f450682621edf Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <Daniel.Karrasch@posteo.de>
Date: Thu, 9 Jun 2022 16:33:10 +0200
Subject: [PATCH 20/23] improve coverage

---
 test/getindex.jl | 32 ++++++++++++++++++++++----------
 1 file changed, 22 insertions(+), 10 deletions(-)

diff --git a/test/getindex.jl b/test/getindex.jl
index a0ff8e9a..8c873f6a 100644
--- a/test/getindex.jl
+++ b/test/getindex.jl
@@ -4,6 +4,7 @@ using LinearMaps: VecOrMatMap, ScaledMap
 
 function test_getindex(A::LinearMap, M::AbstractMatrix)
     @assert size(A) == size(M)
+    mask = rand(Bool, size(A))
     imask = rand(Bool, size(A, 1))
     jmask = rand(Bool, size(A, 2))
     @test A[:] == M[:]
@@ -13,25 +14,37 @@ function test_getindex(A::LinearMap, M::AbstractMatrix)
     @test A[:,1:4] == M[:,1:4]
     @test (lastindex(A, 1), lastindex(A, 2)) == size(A)
     if A isa VecOrMatMap || A isa ScaledMap{<:Any,<:Any,<:VecOrMatMap}
+        @test A[1,1] == M[1,1]
         @test A[1,1:3] == M[1,1:3]
         @test A[1:3,1] == M[1:3,1]
         @test A[2:end,1] == M[2:end,1]
         @test A[1:2,1:3] == M[1:2,1:3]
         @test A[[2,1],1:3] == M[[2,1],1:3]
         @test A[:,:] == M
+        @test A[7] == M[7]
+        @test A[3:7] == M[3:7]
+        @test A[mask] == M[mask]
+        @test A[findall(mask)] == M[findall(mask)]
+        @test A[CartesianIndex(1,1)] == M[CartesianIndex(1,1)]
         @test A[imask, 1] == M[imask, 1]
         @test A[1, jmask] == M[1, jmask]
         @test A[imask, jmask] == M[imask, jmask]
     else
-        @test_throws ErrorException A[1,1:3] == M[1,1:3]
-        @test_throws ErrorException A[1:3,1] == M[1:3,1]
-        @test_throws ErrorException A[2:end,1] == M[2:end,1]
-        @test_throws ErrorException A[1:2,1:3] == M[1:2,1:3]
-        @test_throws ErrorException A[[2,1],1:3] == M[[2,1],1:3]
-        @test_throws ErrorException A[:,:] == M
-        @test_throws ErrorException A[imask, 1] == M[imask, 1]
-        @test_throws ErrorException A[1, jmask] == M[1, jmask]
-        @test_throws ErrorException A[imask, jmask] == M[imask, jmask]
+        @test_throws ErrorException A[1,1]
+        @test_throws ErrorException A[1,1:3]
+        @test_throws ErrorException A[1:3,1]
+        @test_throws ErrorException A[2:end,1]
+        @test_throws ErrorException A[1:2,1:3]
+        @test_throws ErrorException A[[2,1],1:3]
+        @test_throws ErrorException A[:,:]
+        @test_throws ErrorException A[7]
+        @test_throws ErrorException A[3:7]
+        @test_throws ErrorException A[mask]
+        @test_throws ErrorException A[findall(mask)]
+        @test_throws ErrorException A[CartesianIndex(1,1)]
+        @test_throws ErrorException A[imask, 1]
+        @test_throws ErrorException A[1, jmask]
+        @test_throws ErrorException A[imask, jmask]
     end
     @test_throws BoundsError A[6,1]
     @test_throws BoundsError A[1,7]
@@ -56,7 +69,6 @@ end
     struct TwoMap <: LinearMaps.LinearMap{Float64} end
     Base.size(::TwoMap) = (5,5)
     Base.transpose(A::TwoMap) = A
-    LinearMaps._getindex(::TwoMap, i::Integer, j::Integer) = 2.0
     LinearMaps._unsafe_mul!(y::AbstractVector, ::TwoMap, x::AbstractVector) = fill!(y, 2.0*sum(x))
 
     T = TwoMap()

From c5c849a99a32c83bf7e3500610f4b2b369d484d6 Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <Daniel.Karrasch@posteo.de>
Date: Fri, 10 Jun 2022 12:13:58 +0200
Subject: [PATCH 21/23] allow double-colon

---
 docs/src/history.md |  5 +++--
 src/getindex.jl     |  4 +---
 test/getindex.jl    | 11 ++++++-----
 3 files changed, 10 insertions(+), 10 deletions(-)

diff --git a/docs/src/history.md b/docs/src/history.md
index 0d32dcb8..2f465cb8 100644
--- a/docs/src/history.md
+++ b/docs/src/history.md
@@ -26,9 +26,10 @@
     `A` gets (virtually) diagonally appended.
 * An often requested new feature has been added: slicing (i.e., non-scalar indexing) any
   `LinearMap` object via `Base.getindex` overloads. Note, however, that only rather
-  efficient slicing operations are implemented: `A[:,j]`, `A[:,J]`, `A[I,J]`, and `A[:,:]`,
+  efficient complete slicing operations are implemented: `A[:,j]`, `A[:,J]`, and `A[:,:]`,
   where `j::Integer` and `J` is either of type `AbstractVector{<:Integer>}` or an
-  `AbstractVector{Bool}` of appropriate length ("logical slicing").
+  `AbstractVector{Bool}` of appropriate length ("logical slicing"). Partial slicing
+  operations such as `A[I,j]` and `A[I,J]` where `I` is as `J` above are disallowed.
 
   Scalar indexing `A[i::Integer,j::Integer]` as well as other indexing operations that fall
   back on scalar indexing such as logical indexing by some `AbstractMatrix{Bool}`, or
diff --git a/src/getindex.jl b/src/getindex.jl
index fc9e363f..0dac204f 100644
--- a/src/getindex.jl
+++ b/src/getindex.jl
@@ -33,7 +33,6 @@ Base.@propagate_inbounds Base.getindex(A::WrappedMap, I...) = A.lmap[I...]
 # linear indexing
 ########################
 _getindex(A::LinearMap, _) = error("linear indexing of LinearMaps is not supported")
-_getindex(A::LinearMap, ::Base.Slice) = vec(Matrix(A))
 
 ########################
 # Cartesian indexing
@@ -55,8 +54,7 @@ function _getindex(A::LinearMap, i::Integer, J::Base.Slice)
 end
 _getindex(A::LinearMap, I::Indexer, J::Indexer) =
     error("partial two-dimensional slicing of LinearMaps is not supported, consider using A[:,J][I] or A[I,:][J] instead")
-_getindex(A::LinearMap, ::Base.Slice, ::Base.Slice) =
-    error("two-dimensional slicing of LinearMaps is not supported, consider using Matrix(A) or convert(Matrix, A)") 
+_getindex(A::LinearMap, ::Base.Slice, ::Base.Slice) = convert(AbstractMatrix, A)
 _getindex(A::LinearMap, I::Base.Slice, J::Indexer) = __getindex(A, I, J)
 _getindex(A::LinearMap, I::Indexer, J::Base.Slice) = __getindex(A, I, J)
 function __getindex(A, I, J)
diff --git a/test/getindex.jl b/test/getindex.jl
index 8c873f6a..7bf4d830 100644
--- a/test/getindex.jl
+++ b/test/getindex.jl
@@ -7,20 +7,20 @@ function test_getindex(A::LinearMap, M::AbstractMatrix)
     mask = rand(Bool, size(A))
     imask = rand(Bool, size(A, 1))
     jmask = rand(Bool, size(A, 2))
-    @test A[:] == M[:]
     @test A[1,:] == M[1,:]
     @test A[:,1] == M[:,1]
     @test A[1:4,:] == M[1:4,:]
     @test A[:,1:4] == M[:,1:4]
+    @test A[:,:] == M
     @test (lastindex(A, 1), lastindex(A, 2)) == size(A)
     if A isa VecOrMatMap || A isa ScaledMap{<:Any,<:Any,<:VecOrMatMap}
+        @test A[:] == M[:]
         @test A[1,1] == M[1,1]
         @test A[1,1:3] == M[1,1:3]
         @test A[1:3,1] == M[1:3,1]
         @test A[2:end,1] == M[2:end,1]
         @test A[1:2,1:3] == M[1:2,1:3]
         @test A[[2,1],1:3] == M[[2,1],1:3]
-        @test A[:,:] == M
         @test A[7] == M[7]
         @test A[3:7] == M[3:7]
         @test A[mask] == M[mask]
@@ -30,13 +30,13 @@ function test_getindex(A::LinearMap, M::AbstractMatrix)
         @test A[1, jmask] == M[1, jmask]
         @test A[imask, jmask] == M[imask, jmask]
     else
+        @test_throws ErrorException A[:]
         @test_throws ErrorException A[1,1]
         @test_throws ErrorException A[1,1:3]
         @test_throws ErrorException A[1:3,1]
         @test_throws ErrorException A[2:end,1]
         @test_throws ErrorException A[1:2,1:3]
         @test_throws ErrorException A[[2,1],1:3]
-        @test_throws ErrorException A[:,:]
         @test_throws ErrorException A[7]
         @test_throws ErrorException A[3:7]
         @test_throws ErrorException A[mask]
@@ -68,10 +68,11 @@ end
 
     struct TwoMap <: LinearMaps.LinearMap{Float64} end
     Base.size(::TwoMap) = (5,5)
-    Base.transpose(A::TwoMap) = A
     LinearMaps._unsafe_mul!(y::AbstractVector, ::TwoMap, x::AbstractVector) = fill!(y, 2.0*sum(x))
-
     T = TwoMap()
+    @test_throws ErrorException T[1,:]
+
+    Base.transpose(A::TwoMap) = A
     test_getindex(TwoMap(), fill(2.0, size(T)))
 
     MA = rand(ComplexF64, 5, 5)

From 31bfd43c3ab16b061d6714b3485b25f20c158b2b Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <Daniel.Karrasch@posteo.de>
Date: Fri, 10 Jun 2022 14:08:40 +0200
Subject: [PATCH 22/23] more (robust) tests

---
 test/getindex.jl | 12 ++++++++----
 1 file changed, 8 insertions(+), 4 deletions(-)

diff --git a/test/getindex.jl b/test/getindex.jl
index 7bf4d830..9a75af06 100644
--- a/test/getindex.jl
+++ b/test/getindex.jl
@@ -11,6 +11,8 @@ function test_getindex(A::LinearMap, M::AbstractMatrix)
     @test A[:,1] == M[:,1]
     @test A[1:4,:] == M[1:4,:]
     @test A[:,1:4] == M[:,1:4]
+    @test A[[2,1],:] == M[[2,1],:]
+    @test A[:,[2,1]] == M[:,[2,1]]
     @test A[:,:] == M
     @test (lastindex(A, 1), lastindex(A, 2)) == size(A)
     if A isa VecOrMatMap || A isa ScaledMap{<:Any,<:Any,<:VecOrMatMap}
@@ -46,10 +48,10 @@ function test_getindex(A::LinearMap, M::AbstractMatrix)
         @test_throws ErrorException A[1, jmask]
         @test_throws ErrorException A[imask, jmask]
     end
-    @test_throws BoundsError A[6,1]
-    @test_throws BoundsError A[1,7]
-    @test_throws BoundsError A[2,1:7]
-    @test_throws BoundsError A[1:6,2]
+    @test_throws BoundsError A[lastindex(A,1)+1,1]
+    @test_throws BoundsError A[1,lastindex(A,2)+1]
+    @test_throws BoundsError A[2,1:lastindex(A,2)+1]
+    @test_throws BoundsError A[1:lastindex(A,1)+1,2]
     @test_throws BoundsError A[ones(Bool, 2, 2)]
     @test_throws BoundsError A[[true, true], 1]
     @test_throws BoundsError A[1, [true, true]]
@@ -79,6 +81,8 @@ end
     FA = LinearMap{ComplexF64}((y, x) -> mul!(y, MA, x), (y, x) -> mul!(y, MA', x), 5, 5)
     F = LinearMap{ComplexF64}(x -> MA*x, y -> MA'y, 5, 5)
     test_getindex(FA, MA)
+    test_getindex([FA FA], [MA MA])
+    test_getindex([FA; FA], [MA; MA])
     test_getindex(F, MA)
     test_getindex(3FA, 3MA)
     test_getindex(FA + FA, 2MA)

From d378b04fdb2d92db7ee2cd43d6debe55b7f5f40e Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <Daniel.Karrasch@posteo.de>
Date: Sat, 11 Jun 2022 20:03:41 +0200
Subject: [PATCH 23/23] improve coverage

---
 src/getindex.jl  | 15 ++++++++-------
 test/getindex.jl |  2 +-
 2 files changed, 9 insertions(+), 8 deletions(-)

diff --git a/src/getindex.jl b/src/getindex.jl
index 0dac204f..bffacc45 100644
--- a/src/getindex.jl
+++ b/src/getindex.jl
@@ -1,10 +1,10 @@
 const Indexer = AbstractVector{<:Integer}
 
 Base.IndexStyle(::LinearMap) = IndexCartesian()
-# required in Base.to_indices for [:]-indexing
+# required in Base.to_indices for [:]-indexing (only size check)
 Base.eachindex(::IndexLinear, A::LinearMap) = Base.OneTo(length(A))
-Base.lastindex(A::LinearMap) = last(eachindex(IndexLinear(), A))
-Base.firstindex(A::LinearMap) = first(eachindex(IndexLinear(), A))
+# Base.lastindex(A::LinearMap) = last(eachindex(IndexLinear(), A))
+# Base.firstindex(A::LinearMap) = first(eachindex(IndexLinear(), A))
 
 function Base.checkbounds(A::LinearMap, i, j)
     Base.checkbounds_indices(Bool, axes(A), (i, j)) || throw(BoundsError(A, (i, j)))
@@ -35,15 +35,18 @@ Base.@propagate_inbounds Base.getindex(A::WrappedMap, I...) = A.lmap[I...]
 _getindex(A::LinearMap, _) = error("linear indexing of LinearMaps is not supported")
 
 ########################
-# Cartesian indexing
+# Cartesian indexing (partial slicing is not supported)
 ########################
 _getindex(A::LinearMap, i::Integer, j::Integer) =
     error("scalar indexing of LinearMaps is not supported, consider using A[:,j][i] instead")
 _getindex(A::LinearMap, I::Indexer, j::Integer) =
     error("partial vertical slicing of LinearMaps is not supported, consider using A[:,j][I] instead")
-_getindex(A::LinearMap, ::Base.Slice, j::Integer) = A*unitvec(A, 2, j)
 _getindex(A::LinearMap, i::Integer, J::Indexer) =
     error("partial horizontal slicing of LinearMaps is not supported, consider using A[i,:][J] instead")
+_getindex(A::LinearMap, I::Indexer, J::Indexer) =
+    error("partial two-dimensional slicing of LinearMaps is not supported, consider using A[:,J][I] or A[I,:][J] instead")
+
+_getindex(A::LinearMap, ::Base.Slice, j::Integer) = A*unitvec(A, 2, j)
 function _getindex(A::LinearMap, i::Integer, J::Base.Slice)
     try
         # requires adjoint action to be defined
@@ -52,8 +55,6 @@ function _getindex(A::LinearMap, i::Integer, J::Base.Slice)
         error("horizontal slicing A[$i,:] requires the adjoint of $(typeof(A)) to be defined")
     end
 end
-_getindex(A::LinearMap, I::Indexer, J::Indexer) =
-    error("partial two-dimensional slicing of LinearMaps is not supported, consider using A[:,J][I] or A[I,:][J] instead")
 _getindex(A::LinearMap, ::Base.Slice, ::Base.Slice) = convert(AbstractMatrix, A)
 _getindex(A::LinearMap, I::Base.Slice, J::Indexer) = __getindex(A, I, J)
 _getindex(A::LinearMap, I::Indexer, J::Base.Slice) = __getindex(A, I, J)
diff --git a/test/getindex.jl b/test/getindex.jl
index 9a75af06..7f1a9ba1 100644
--- a/test/getindex.jl
+++ b/test/getindex.jl
@@ -9,7 +9,7 @@ function test_getindex(A::LinearMap, M::AbstractMatrix)
     jmask = rand(Bool, size(A, 2))
     @test A[1,:] == M[1,:]
     @test A[:,1] == M[:,1]
-    @test A[1:4,:] == M[1:4,:]
+    @test A[1:lastindex(A,1)-2,:] == M[1:lastindex(A,1)-2,:]
     @test A[:,1:4] == M[:,1:4]
     @test A[[2,1],:] == M[[2,1],:]
     @test A[:,[2,1]] == M[:,[2,1]]