Skip to content

Commit 57d78c4

Browse files
JeffFesslerJeffFessler
committed
Support mul! for left vectors
1 parent e45289d commit 57d78c4

13 files changed

+157
-86
lines changed

‎src/LinearMaps.jl

Lines changed: 19 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -4,6 +4,7 @@ export LinearMap
44
export ⊗, kronsum, ⊕
55

66
using LinearAlgebra
7+
import LinearAlgebra: AdjointAbsVec, TransposeAbsVec
78
using SparseArrays
89

910
if VERSION < v"1.2-"
@@ -18,6 +19,10 @@ abstract type LinearMap{T} end
1819
const MapOrMatrix{T} = Union{LinearMap{T},AbstractMatrix{T}}
1920
const RealOrComplex = Union{Real,Complex}
2021

22+
# valid types for left vector multiplication:
23+
const VecIn = Union{AdjointAbsVec, TransposeAbsVec}
24+
const VecOut = Union{AbstractVector, AdjointAbsVec, TransposeAbsVec}
25+
2126
Base.eltype(::LinearMap{T}) where {T} = T
2227

2328
abstract type MulStyle end
@@ -48,7 +53,7 @@ Base.size(A::LinearMap, n) = (n==1 || n==2 ? size(A)[n] : error("LinearMap objec
4853
Base.length(A::LinearMap) = size(A)[1] * size(A)[2]
4954

5055
# check dimension consistency for y = A*x and Y = A*X
51-
function check_dim_mul(y::AbstractVector, A::LinearMap, x::AbstractVector)
56+
function check_dim_mul(y::VecOut, A::LinearMap, x::AbstractVector)
5257
# @info "checked vector dimensions" # uncomment for testing
5358
m, n = size(A)
5459
(m == length(y) && n == length(x)) || throw(DimensionMismatch("mul!"))
@@ -61,6 +66,13 @@ function check_dim_mul(Y::AbstractMatrix, A::LinearMap, X::AbstractMatrix)
6166
return nothing
6267
end
6368

69+
# check dimension consistency for left multiplication x = y'*A
70+
function check_dim_mul(x::V, y::V, A::LinearMap) where {V <: VecIn}
71+
m, n = size(A)
72+
((1,m) == size(y) && (1,n) == size(x)) || throw(DimensionMismatch("left mul!"))
73+
return nothing
74+
end
75+
6476
# conversion of AbstractMatrix to LinearMap
6577
convert_to_lmaps_(A::AbstractMatrix) = LinearMap(A)
6678
convert_to_lmaps_(A::LinearMap) = A
@@ -73,11 +85,11 @@ function Base.:(*)(A::LinearMap, x::AbstractVector)
7385
size(A, 2) == length(x) || throw(DimensionMismatch("mul!"))
7486
return @inbounds A_mul_B!(similar(x, promote_type(eltype(A), eltype(x)), size(A, 1)), A, x)
7587
end
76-
function LinearAlgebra.mul!(y::AbstractVector, A::LinearMap, x::AbstractVector)
88+
function LinearAlgebra.mul!(y::VecOut, A::LinearMap, x::AbstractVector)
7789
@boundscheck check_dim_mul(y, A, x)
7890
return @inbounds A_mul_B!(y, A, x)
7991
end
80-
function LinearAlgebra.mul!(y::AbstractVector, A::LinearMap, x::AbstractVector, α::Number, β::Number)
92+
function LinearAlgebra.mul!(y::VecOut, A::LinearMap, x::AbstractVector, α::Number, β::Number)
8193
@boundscheck check_dim_mul(y, A, x)
8294
if isone(α)
8395
iszero(β) && (A_mul_B!(y, A, x); return y)
@@ -114,10 +126,11 @@ Base.@propagate_inbounds function LinearAlgebra.mul!(Y::AbstractMatrix, A::Linea
114126
return Y
115127
end
116128

117-
A_mul_B!(y::AbstractVector, A::AbstractMatrix, x::AbstractVector) = mul!(y, A, x)
118-
At_mul_B!(y::AbstractVector, A::AbstractMatrix, x::AbstractVector) = mul!(y, transpose(A), x)
119-
Ac_mul_B!(y::AbstractVector, A::AbstractMatrix, x::AbstractVector) = mul!(y, adjoint(A), x)
129+
A_mul_B!(y::VecOut, A::AbstractMatrix, x::AbstractVector) = mul!(y, A, x)
130+
At_mul_B!(y::VecOut, A::AbstractMatrix, x::AbstractVector) = mul!(y, transpose(A), x)
131+
Ac_mul_B!(y::VecOut, A::AbstractMatrix, x::AbstractVector) = mul!(y, adjoint(A), x)
120132

133+
include("left.jl") # left multiplication by a transpose or adjoint vector
121134
include("wrappedmap.jl") # wrap a matrix of linear map in a new type, thereby allowing to alter its properties
122135
include("uniformscalingmap.jl") # the uniform scaling map, to be able to make linear combinations of LinearMap objects and multiples of I
123136
include("transpose.jl") # transposing linear maps

‎src/blockmap.jl

Lines changed: 8 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -382,19 +382,19 @@ end
382382
# multiplication with vectors & matrices
383383
############
384384

385-
Base.@propagate_inbounds A_mul_B!(y::AbstractVector, A::BlockMap, x::AbstractVector) =
385+
Base.@propagate_inbounds A_mul_B!(y::VecOut, A::BlockMap, x::AbstractVector) =
386386
mul!(y, A, x, true, false)
387387

388-
Base.@propagate_inbounds A_mul_B!(y::AbstractVector, A::TransposeMap{<:Any,<:BlockMap}, x::AbstractVector) =
388+
Base.@propagate_inbounds A_mul_B!(y::VecOut, A::TransposeMap{<:Any,<:BlockMap}, x::AbstractVector) =
389389
mul!(y, A, x, true, false)
390390

391-
Base.@propagate_inbounds At_mul_B!(y::AbstractVector, A::BlockMap, x::AbstractVector) =
391+
Base.@propagate_inbounds At_mul_B!(y::VecOut, A::BlockMap, x::AbstractVector) =
392392
mul!(y, transpose(A), x, true, false)
393393

394-
Base.@propagate_inbounds A_mul_B!(y::AbstractVector, A::AdjointMap{<:Any,<:BlockMap}, x::AbstractVector) =
394+
Base.@propagate_inbounds A_mul_B!(y::VecOut, A::AdjointMap{<:Any,<:BlockMap}, x::AbstractVector) =
395395
mul!(y, A, x, true, false)
396396

397-
Base.@propagate_inbounds Ac_mul_B!(y::AbstractVector, A::BlockMap, x::AbstractVector) =
397+
Base.@propagate_inbounds Ac_mul_B!(y::VecOut, A::BlockMap, x::AbstractVector) =
398398
mul!(y, adjoint(A), x, true, false)
399399

400400
for Atype in (AbstractVector, AbstractMatrix)
@@ -495,13 +495,13 @@ LinearAlgebra.transpose(A::BlockDiagonalMap{T}) where {T} = BlockDiagonalMap{T}(
495495

496496
Base.:(==)(A::BlockDiagonalMap, B::BlockDiagonalMap) = (eltype(A) == eltype(B) && A.maps == B.maps)
497497

498-
Base.@propagate_inbounds A_mul_B!(y::AbstractVector, A::BlockDiagonalMap, x::AbstractVector) =
498+
Base.@propagate_inbounds A_mul_B!(y::VecOut, A::BlockDiagonalMap, x::AbstractVector) =
499499
mul!(y, A, x, true, false)
500500

501-
Base.@propagate_inbounds At_mul_B!(y::AbstractVector, A::BlockDiagonalMap, x::AbstractVector) =
501+
Base.@propagate_inbounds At_mul_B!(y::VecOut, A::BlockDiagonalMap, x::AbstractVector) =
502502
mul!(y, transpose(A), x, true, false)
503503

504-
Base.@propagate_inbounds Ac_mul_B!(y::AbstractVector, A::BlockDiagonalMap, x::AbstractVector) =
504+
Base.@propagate_inbounds Ac_mul_B!(y::VecOut, A::BlockDiagonalMap, x::AbstractVector) =
505505
mul!(y, adjoint(A), x, true, false)
506506

507507
for Atype in (AbstractVector, AbstractMatrix)

‎src/composition.jl

Lines changed: 7 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -129,23 +129,23 @@ LinearAlgebra.adjoint(A::CompositeMap{T}) where {T} = CompositeMap{T}(map(adjo
129129
Base.:(==)(A::CompositeMap, B::CompositeMap) = (eltype(A) == eltype(B) && A.maps == B.maps)
130130

131131
# multiplication with vectors
132-
function A_mul_B!(y::AbstractVector, A::CompositeMap{T,<:Tuple{LinearMap}}, x::AbstractVector) where {T}
132+
function A_mul_B!(y::VecOut, A::CompositeMap{T,<:Tuple{LinearMap}}, x::AbstractVector) where {T}
133133
return A_mul_B!(y, A.maps[1], x)
134134
end
135-
function A_mul_B!(y::AbstractVector, A::CompositeMap{T,<:Tuple{LinearMap,LinearMap}}, x::AbstractVector) where {T}
135+
function A_mul_B!(y::VecOut, A::CompositeMap{T,<:Tuple{LinearMap,LinearMap}}, x::AbstractVector) where {T}
136136
_compositemul!(y, A, x, similar(y, size(A.maps[1], 1)))
137137
end
138-
function A_mul_B!(y::AbstractVector, A::CompositeMap{T,<:Tuple{Vararg{LinearMap}}}, x::AbstractVector) where {T}
138+
function A_mul_B!(y::VecOut, A::CompositeMap{T,<:Tuple{Vararg{LinearMap}}}, x::AbstractVector) where {T}
139139
_compositemul!(y, A, x, similar(y, size(A.maps[1], 1)), similar(y, size(A.maps[2], 1)))
140140
end
141141

142-
function _compositemul!(y::AbstractVector, A::CompositeMap{T,<:Tuple{LinearMap,LinearMap}}, x::AbstractVector, z::AbstractVector) where {T}
142+
function _compositemul!(y::VecOut, A::CompositeMap{T,<:Tuple{LinearMap,LinearMap}}, x::AbstractVector, z::AbstractVector) where {T}
143143
# no size checking, will be done by individual maps
144144
A_mul_B!(z, A.maps[1], x)
145145
A_mul_B!(y, A.maps[2], z)
146146
return y
147147
end
148-
function _compositemul!(y::AbstractVector, A::CompositeMap{T,<:Tuple{Vararg{LinearMap}}}, x::AbstractVector, source::AbstractVector, dest::AbstractVector) where {T}
148+
function _compositemul!(y::VecOut, A::CompositeMap{T,<:Tuple{Vararg{LinearMap}}}, x::AbstractVector, source::AbstractVector, dest::AbstractVector) where {T}
149149
# no size checking, will be done by individual maps
150150
N = length(A.maps)
151151
A_mul_B!(source, A.maps[1], x)
@@ -166,6 +166,6 @@ function _compositemul!(y::AbstractVector, A::CompositeMap{T,<:Tuple{Vararg{Line
166166
return y
167167
end
168168

169-
At_mul_B!(y::AbstractVector, A::CompositeMap, x::AbstractVector) = A_mul_B!(y, transpose(A), x)
169+
At_mul_B!(y::VecOut, A::CompositeMap, x::AbstractVector) = A_mul_B!(y, transpose(A), x)
170170

171-
Ac_mul_B!(y::AbstractVector, A::CompositeMap, x::AbstractVector) = A_mul_B!(y, adjoint(A), x)
171+
Ac_mul_B!(y::VecOut, A::CompositeMap, x::AbstractVector) = A_mul_B!(y, adjoint(A), x)

‎src/functionmap.jl

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -107,13 +107,13 @@ function Base.:(*)(A::TransposeMap{<:Any,<:FunctionMap}, x::AbstractVector)
107107
end
108108
end
109109

110-
function A_mul_B!(y::AbstractVector, A::FunctionMap, x::AbstractVector)
110+
function A_mul_B!(y::VecOut, A::FunctionMap, x::AbstractVector)
111111
(length(x) == A.N && length(y) == A.M) || throw(DimensionMismatch("A_mul_B!"))
112112
ismutating(A) ? A.f(y, x) : copyto!(y, A.f(x))
113113
return y
114114
end
115115

116-
function At_mul_B!(y::AbstractVector, A::FunctionMap, x::AbstractVector)
116+
function At_mul_B!(y::VecOut, A::FunctionMap, x::AbstractVector)
117117
(issymmetric(A) || (isreal(A) && ishermitian(A))) && return A_mul_B!(y, A, x)
118118
(length(x) == A.M && length(y) == A.N) || throw(DimensionMismatch("At_mul_B!"))
119119
if A.fc !== nothing
@@ -135,7 +135,7 @@ function At_mul_B!(y::AbstractVector, A::FunctionMap, x::AbstractVector)
135135
end
136136
end
137137

138-
function Ac_mul_B!(y::AbstractVector, A::FunctionMap, x::AbstractVector)
138+
function Ac_mul_B!(y::VecOut, A::FunctionMap, x::AbstractVector)
139139
ishermitian(A) && return A_mul_B!(y, A, x)
140140
(length(x) == A.M && length(y) == A.N) || throw(DimensionMismatch("Ac_mul_B!"))
141141
if A.fc !== nothing

‎src/kronecker.jl

Lines changed: 9 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -131,15 +131,15 @@ end
131131
# multiplication with vectors
132132
#################
133133

134-
Base.@propagate_inbounds function A_mul_B!(y::AbstractVector, L::KroneckerMap{T,<:NTuple{2,LinearMap}}, x::AbstractVector) where {T}
134+
Base.@propagate_inbounds function A_mul_B!(y::VecOut, L::KroneckerMap{T,<:NTuple{2,LinearMap}}, x::AbstractVector) where {T}
135135
require_one_based_indexing(y)
136136
@boundscheck check_dim_mul(y, L, x)
137137
A, B = L.maps
138138
X = LinearMap(reshape(x, (size(B, 2), size(A, 2))); issymmetric=false, ishermitian=false, isposdef=false)
139139
_kronmul!(y, B, X, transpose(A), T)
140140
return y
141141
end
142-
Base.@propagate_inbounds function A_mul_B!(y::AbstractVector, L::KroneckerMap{T}, x::AbstractVector) where {T}
142+
Base.@propagate_inbounds function A_mul_B!(y::VecOut, L::KroneckerMap{T}, x::AbstractVector) where {T}
143143
require_one_based_indexing(y)
144144
@boundscheck check_dim_mul(y, L, x)
145145
A = first(L.maps)
@@ -150,7 +150,7 @@ Base.@propagate_inbounds function A_mul_B!(y::AbstractVector, L::KroneckerMap{T}
150150
end
151151
# mixed-product rule, prefer the right if possible
152152
# (A₁ ⊗ A₂ ⊗ ... ⊗ Aᵣ) * (B₁ ⊗ B₂ ⊗ ... ⊗ Bᵣ) = (A₁B₁) ⊗ (A₂B₂) ⊗ ... ⊗ (AᵣBᵣ)
153-
Base.@propagate_inbounds function A_mul_B!(y::AbstractVector, L::CompositeMap{<:Any,<:Tuple{KroneckerMap,KroneckerMap}}, x::AbstractVector)
153+
Base.@propagate_inbounds function A_mul_B!(y::VecOut, L::CompositeMap{<:Any,<:Tuple{KroneckerMap,KroneckerMap}}, x::AbstractVector)
154154
B, A = L.maps
155155
if length(A.maps) == length(B.maps) && all(M -> check_dim_mul(M[1], M[2]), zip(A.maps, B.maps))
156156
A_mul_B!(y, kron(map(*, A.maps, B.maps)...), x)
@@ -160,7 +160,7 @@ Base.@propagate_inbounds function A_mul_B!(y::AbstractVector, L::CompositeMap{<:
160160
end
161161
# mixed-product rule, prefer the right if possible
162162
# (A₁ ⊗ B₁)*(A₂⊗B₂)*...*(Aᵣ⊗Bᵣ) = (A₁*A₂*...*Aᵣ) ⊗ (B₁*B₂*...*Bᵣ)
163-
Base.@propagate_inbounds function A_mul_B!(y::AbstractVector, L::CompositeMap{T,<:Tuple{Vararg{KroneckerMap{<:Any,<:Tuple{LinearMap,LinearMap}}}}}, x::AbstractVector) where {T}
163+
Base.@propagate_inbounds function A_mul_B!(y::VecOut, L::CompositeMap{T,<:Tuple{Vararg{KroneckerMap{<:Any,<:Tuple{LinearMap,LinearMap}}}}}, x::AbstractVector) where {T}
164164
As = map(AB -> AB.maps[1], L.maps)
165165
Bs = map(AB -> AB.maps[2], L.maps)
166166
As1, As2 = Base.front(As), Base.tail(As)
@@ -173,10 +173,10 @@ Base.@propagate_inbounds function A_mul_B!(y::AbstractVector, L::CompositeMap{T,
173173
end
174174
end
175175

176-
Base.@propagate_inbounds At_mul_B!(y::AbstractVector, A::KroneckerMap, x::AbstractVector) =
176+
Base.@propagate_inbounds At_mul_B!(y::VecOut, A::KroneckerMap, x::AbstractVector) =
177177
A_mul_B!(y, transpose(A), x)
178178

179-
Base.@propagate_inbounds Ac_mul_B!(y::AbstractVector, A::KroneckerMap, x::AbstractVector) =
179+
Base.@propagate_inbounds Ac_mul_B!(y::VecOut, A::KroneckerMap, x::AbstractVector) =
180180
A_mul_B!(y, adjoint(A), x)
181181

182182
###############
@@ -261,7 +261,7 @@ LinearAlgebra.transpose(A::KroneckerSumMap{T}) where {T} = KroneckerSumMap{T}(ma
261261

262262
Base.:(==)(A::KroneckerSumMap, B::KroneckerSumMap) = (eltype(A) == eltype(B) && A.maps == B.maps)
263263

264-
Base.@propagate_inbounds function A_mul_B!(y::AbstractVector, L::KroneckerSumMap, x::AbstractVector)
264+
Base.@propagate_inbounds function A_mul_B!(y::VecOut, L::KroneckerSumMap, x::AbstractVector)
265265
@boundscheck check_dim_mul(y, L, x)
266266
A, B = L.maps
267267
ma, na = size(A)
@@ -273,8 +273,8 @@ Base.@propagate_inbounds function A_mul_B!(y::AbstractVector, L::KroneckerSumMap
273273
return y
274274
end
275275

276-
Base.@propagate_inbounds At_mul_B!(y::AbstractVector, A::KroneckerSumMap, x::AbstractVector) =
276+
Base.@propagate_inbounds At_mul_B!(y::VecOut, A::KroneckerSumMap, x::AbstractVector) =
277277
A_mul_B!(y, transpose(A), x)
278278

279-
Base.@propagate_inbounds Ac_mul_B!(y::AbstractVector, A::KroneckerSumMap, x::AbstractVector) =
279+
Base.@propagate_inbounds Ac_mul_B!(y::VecOut, A::KroneckerSumMap, x::AbstractVector) =
280280
A_mul_B!(y, adjoint(A), x)

‎src/left.jl

Lines changed: 35 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,35 @@
1+
# left.jl
2+
# "special cases" related left vector multiplication like x = y'*A
3+
# The subtlety is that y' is a AdjointAbsVec or a TransposeAbsVec
4+
# which is a subtype of AbstractMatrix but it is really "vector like"
5+
# so we want handle it like a vector (Issue#99)
6+
# So this is an exception to the left multiplication rule by a AbstractMatrix
7+
# that usually makes a WrappedMap.
8+
# The "transpose" versions may be of dubious use, but the adjoint versions
9+
# are useful for ensuring that (y'A)*x ≈ y'*(A*x) are both scalars.
10+
11+
import LinearAlgebra: AdjointAbsVec, TransposeAbsVec
12+
13+
Base.:(*)(y::AdjointAbsVec, A::LinearMap) = adjoint(*(A', y'))
14+
Base.:(*)(y::TransposeAbsVec, A::LinearMap) = transpose(transpose(A) * transpose(y))
15+
16+
# mul!(x, y', A)
17+
LinearAlgebra.mul!(x::AdjointAbsVec, y::AdjointAbsVec, A::LinearMap) =
18+
mul!(x, y, A, true, false)
19+
20+
# not sure if we need bounds checks and propagate inbounds stuff here
21+
22+
# mul!(x, y', A, α, β)
23+
function LinearAlgebra.mul!(x::AdjointAbsVec, y::AdjointAbsVec, A::LinearMap,
24+
α::Number, β::Number)
25+
check_dim_mul(x, y, A)
26+
mul!(x, A', y', α, β)
27+
return conj!(x)
28+
end
29+
30+
# mul!(x, transpose(y), A, α, β)
31+
function LinearAlgebra.mul!(x::TransposeAbsVec, y::TransposeAbsVec, A::LinearMap,
32+
α::Number, β::Number)
33+
check_dim_mul(x, y, A)
34+
return mul!(x, transpose(A), transpose(y), α, β)
35+
end

‎src/linearcombination.jl

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -117,8 +117,8 @@ end
117117
return y
118118
end
119119

120-
A_mul_B!(y::AbstractVector, A::LinearCombination, x::AbstractVector) = mul!(y, A, x, true, false)
120+
A_mul_B!(y::VecOut, A::LinearCombination, x::AbstractVector) = mul!(y, A, x, true, false)
121121

122-
At_mul_B!(y::AbstractVector, A::LinearCombination, x::AbstractVector) = mul!(y, transpose(A), x, true, false)
122+
At_mul_B!(y::VecOut, A::LinearCombination, x::AbstractVector) = mul!(y, transpose(A), x, true, false)
123123

124-
Ac_mul_B!(y::AbstractVector, A::LinearCombination, x::AbstractVector) = mul!(y, adjoint(A), x, true, false)
124+
Ac_mul_B!(y::VecOut, A::LinearCombination, x::AbstractVector) = mul!(y, adjoint(A), x, true, false)

‎src/scaledmap.jl

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -54,15 +54,15 @@ Base.:(*)(A::ScaledMap, B::LinearMap) = A.λ * (A.lmap * B)
5454
Base.:(*)(A::LinearMap, B::ScaledMap) = (A * B.lmap) * B.λ
5555

5656
# multiplication with vectors
57-
function A_mul_B!(y::AbstractVector, A::ScaledMap, x::AbstractVector)
57+
function A_mul_B!(y::VecOut, A::ScaledMap, x::AbstractVector)
5858
# no size checking, will be done by map
5959
mul!(y, A.lmap, x, A.λ, false)
6060
end
6161

62-
function LinearAlgebra.mul!(y::AbstractVector, A::ScaledMap, x::AbstractVector, α::Number, β::Number)
62+
function LinearAlgebra.mul!(y::VecOut, A::ScaledMap, x::AbstractVector, α::Number, β::Number)
6363
# no size checking, will be done by map
6464
mul!(y, A.lmap, x, A.λ * α, β)
6565
end
6666

67-
At_mul_B!(y::AbstractVector, A::ScaledMap, x::AbstractVector) = A_mul_B!(y, transpose(A), x)
68-
Ac_mul_B!(y::AbstractVector, A::ScaledMap, x::AbstractVector) = A_mul_B!(y, adjoint(A), x)
67+
At_mul_B!(y::VecOut, A::ScaledMap, x::AbstractVector) = A_mul_B!(y, transpose(A), x)
68+
Ac_mul_B!(y::VecOut, A::ScaledMap, x::AbstractVector) = A_mul_B!(y, adjoint(A), x)

‎src/transpose.jl

Lines changed: 6 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -32,18 +32,18 @@ Base.:(==)(A::LinearMap, B::TransposeMap) = issymmetric(A) && B.lmap == A
3232
Base.:(==)(A::LinearMap, B::AdjointMap) = ishermitian(A) && B.lmap == A
3333

3434
# multiplication with vector
35-
A_mul_B!(y::AbstractVector, A::TransposeMap, x::AbstractVector) =
35+
A_mul_B!(y::VecOut, A::TransposeMap, x::AbstractVector) =
3636
issymmetric(A.lmap) ? A_mul_B!(y, A.lmap, x) : At_mul_B!(y, A.lmap, x)
3737

38-
At_mul_B!(y::AbstractVector, A::TransposeMap, x::AbstractVector) = A_mul_B!(y, A.lmap, x)
38+
At_mul_B!(y::VecOut, A::TransposeMap, x::AbstractVector) = A_mul_B!(y, A.lmap, x)
3939

40-
Ac_mul_B!(y::AbstractVector, A::TransposeMap, x::AbstractVector) =
40+
Ac_mul_B!(y::VecOut, A::TransposeMap, x::AbstractVector) =
4141
isreal(A.lmap) ? A_mul_B!(y, A.lmap, x) : (A_mul_B!(y, A.lmap, conj(x)); conj!(y))
4242

43-
A_mul_B!(y::AbstractVector, A::AdjointMap, x::AbstractVector) =
43+
A_mul_B!(y::VecOut, A::AdjointMap, x::AbstractVector) =
4444
ishermitian(A.lmap) ? A_mul_B!(y, A.lmap, x) : Ac_mul_B!(y, A.lmap, x)
4545

46-
At_mul_B!(y::AbstractVector, A::AdjointMap, x::AbstractVector) =
46+
At_mul_B!(y::VecOut, A::AdjointMap, x::AbstractVector) =
4747
isreal(A.lmap) ? A_mul_B!(y, A.lmap, x) : (A_mul_B!(y, A.lmap, conj(x)); conj!(y))
4848

49-
Ac_mul_B!(y::AbstractVector, A::AdjointMap, x::AbstractVector) = A_mul_B!(y, A.lmap, x)
49+
Ac_mul_B!(y::VecOut, A::AdjointMap, x::AbstractVector) = A_mul_B!(y, A.lmap, x)

‎src/wrappedmap.jl

Lines changed: 3 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,3 @@
1-
import LinearAlgebra: AdjointAbsVec, TransposeAbsVec
2-
31
struct WrappedMap{T, A<:MapOrMatrix} <: LinearMap{T}
42
lmap::A
53
_issymmetric::Bool
@@ -39,13 +37,13 @@ LinearAlgebra.ishermitian(A::WrappedMap) = A._ishermitian
3937
LinearAlgebra.isposdef(A::WrappedMap) = A._isposdef
4038

4139
# multiplication with vectors & matrices
42-
A_mul_B!(y::AbstractVector, A::WrappedMap, x::AbstractVector) = A_mul_B!(y, A.lmap, x)
40+
A_mul_B!(y::VecOut, A::WrappedMap, x::AbstractVector) = A_mul_B!(y, A.lmap, x)
4341
Base.:(*)(A::WrappedMap, x::AbstractVector) = *(A.lmap, x)
4442

45-
At_mul_B!(y::AbstractVector, A::WrappedMap, x::AbstractVector) =
43+
At_mul_B!(y::VecOut, A::WrappedMap, x::AbstractVector) =
4644
(issymmetric(A) || (isreal(A) && ishermitian(A))) ? A_mul_B!(y, A.lmap, x) : At_mul_B!(y, A.lmap, x)
4745

48-
Ac_mul_B!(y::AbstractVector, A::WrappedMap, x::AbstractVector) =
46+
Ac_mul_B!(y::VecOut, A::WrappedMap, x::AbstractVector) =
4947
ishermitian(A) ? A_mul_B!(y, A.lmap, x) : Ac_mul_B!(y, A.lmap, x)
5048

5149
if VERSION ≥ v"1.3.0-alpha.115"
@@ -73,8 +71,3 @@ Base.:(-)(A₁::AbstractMatrix, A₂::LinearMap) = -(WrappedMap(A₁), A₂)
7371

7472
Base.:(*)(A₁::LinearMap, A₂::AbstractMatrix) = *(A₁, WrappedMap(A₂))
7573
Base.:(*)(A₁::AbstractMatrix, A₂::LinearMap) = *(WrappedMap(A₁), A₂)
76-
77-
# An AdjointAbsVec isa AbstractMatrix but we handle it like a vector (Issue#99)
78-
# which is an exception to the left multiplication rule that makes a WrappedMap
79-
Base.:(*)(y::AdjointAbsVec, A::LinearMap) = adjoint(*(A', y'))
80-
Base.:(*)(y::TransposeAbsVec, A::LinearMap) = transpose(transpose(A) * transpose(y))

0 commit comments

Comments
 (0)