Skip to content

Left vector multiply #100

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 18 commits into from
Aug 24, 2020
Merged
Show file tree
Hide file tree
Changes from 7 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion Project.toml
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
name = "LinearMaps"
uuid = "7a12625a-238d-50fd-b39a-03d52299707e"
version = "2.7.0"
version = "2.8.0"

[deps]
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
Expand Down
4 changes: 4 additions & 0 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,10 @@ A Julia package for defining and working with linear maps, also known as linear
transformations or linear operators acting on vectors. The only requirement for
a LinearMap is that it can act on a vector (by multiplication) efficiently.

## What's new in v2.8
* Left multiplying by a transpose or adjoint vector (e.g., `y'*A`)
produces a transpose or adjoint vector output, rather than a `LinearMap`.

## What's new in v2.7
* Potential reduction of memory allocations in multiplication of
`LinearCombination`s, `BlockMap`s, and real- or complex-scaled `LinearMap`s.
Expand Down
25 changes: 19 additions & 6 deletions src/LinearMaps.jl
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@ export LinearMap
export ⊗, kronsum, ⊕

using LinearAlgebra
import LinearAlgebra: AdjointAbsVec, TransposeAbsVec
using SparseArrays

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

# valid types for left vector multiplication:
const VecIn = Union{AdjointAbsVec, TransposeAbsVec}
const VecOut = Union{AbstractVector, AdjointAbsVec, TransposeAbsVec}

Base.eltype(::LinearMap{T}) where {T} = T

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

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

# check dimension consistency for left multiplication x = y'*A
function check_dim_mul(x::V, y::V, A::LinearMap) where {V <: VecIn}
m, n = size(A)
((1,m) == size(y) && (1,n) == size(x)) || throw(DimensionMismatch("left mul!"))
return nothing
end

# conversion of AbstractMatrix to LinearMap
convert_to_lmaps_(A::AbstractMatrix) = LinearMap(A)
convert_to_lmaps_(A::LinearMap) = A
Expand All @@ -73,11 +85,11 @@ function Base.:(*)(A::LinearMap, x::AbstractVector)
size(A, 2) == length(x) || throw(DimensionMismatch("mul!"))
return @inbounds A_mul_B!(similar(x, promote_type(eltype(A), eltype(x)), size(A, 1)), A, x)
end
function LinearAlgebra.mul!(y::AbstractVector, A::LinearMap, x::AbstractVector)
function LinearAlgebra.mul!(y::VecOut, A::LinearMap, x::AbstractVector)
@boundscheck check_dim_mul(y, A, x)
return @inbounds A_mul_B!(y, A, x)
end
function LinearAlgebra.mul!(y::AbstractVector, A::LinearMap, x::AbstractVector, α::Number, β::Number)
function LinearAlgebra.mul!(y::VecOut, A::LinearMap, x::AbstractVector, α::Number, β::Number)
@boundscheck check_dim_mul(y, A, x)
if isone(α)
iszero(β) && (A_mul_B!(y, A, x); return y)
Expand Down Expand Up @@ -114,10 +126,11 @@ Base.@propagate_inbounds function LinearAlgebra.mul!(Y::AbstractMatrix, A::Linea
return Y
end

A_mul_B!(y::AbstractVector, A::AbstractMatrix, x::AbstractVector) = mul!(y, A, x)
At_mul_B!(y::AbstractVector, A::AbstractMatrix, x::AbstractVector) = mul!(y, transpose(A), x)
Ac_mul_B!(y::AbstractVector, A::AbstractMatrix, x::AbstractVector) = mul!(y, adjoint(A), x)
A_mul_B!(y::VecOut, A::AbstractMatrix, x::AbstractVector) = mul!(y, A, x)
At_mul_B!(y::VecOut, A::AbstractMatrix, x::AbstractVector) = mul!(y, transpose(A), x)
Ac_mul_B!(y::VecOut, A::AbstractMatrix, x::AbstractVector) = mul!(y, adjoint(A), x)

include("left.jl") # left multiplication by a transpose or adjoint vector
include("wrappedmap.jl") # wrap a matrix of linear map in a new type, thereby allowing to alter its properties
include("uniformscalingmap.jl") # the uniform scaling map, to be able to make linear combinations of LinearMap objects and multiples of I
include("transpose.jl") # transposing linear maps
Expand Down
16 changes: 8 additions & 8 deletions src/blockmap.jl
Original file line number Diff line number Diff line change
Expand Up @@ -382,19 +382,19 @@ end
# multiplication with vectors & matrices
############

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

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

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

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

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

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

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

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

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

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

for Atype in (AbstractVector, AbstractMatrix)
Expand Down
14 changes: 7 additions & 7 deletions src/composition.jl
Original file line number Diff line number Diff line change
Expand Up @@ -129,23 +129,23 @@ LinearAlgebra.adjoint(A::CompositeMap{T}) where {T} = CompositeMap{T}(map(adjo
Base.:(==)(A::CompositeMap, B::CompositeMap) = (eltype(A) == eltype(B) && A.maps == B.maps)

# multiplication with vectors
function A_mul_B!(y::AbstractVector, A::CompositeMap{T,<:Tuple{LinearMap}}, x::AbstractVector) where {T}
function A_mul_B!(y::VecOut, A::CompositeMap{T,<:Tuple{LinearMap}}, x::AbstractVector) where {T}
return A_mul_B!(y, A.maps[1], x)
end
function A_mul_B!(y::AbstractVector, A::CompositeMap{T,<:Tuple{LinearMap,LinearMap}}, x::AbstractVector) where {T}
function A_mul_B!(y::VecOut, A::CompositeMap{T,<:Tuple{LinearMap,LinearMap}}, x::AbstractVector) where {T}
_compositemul!(y, A, x, similar(y, size(A.maps[1], 1)))
end
function A_mul_B!(y::AbstractVector, A::CompositeMap{T,<:Tuple{Vararg{LinearMap}}}, x::AbstractVector) where {T}
function A_mul_B!(y::VecOut, A::CompositeMap{T,<:Tuple{Vararg{LinearMap}}}, x::AbstractVector) where {T}
_compositemul!(y, A, x, similar(y, size(A.maps[1], 1)), similar(y, size(A.maps[2], 1)))
end

function _compositemul!(y::AbstractVector, A::CompositeMap{T,<:Tuple{LinearMap,LinearMap}}, x::AbstractVector, z::AbstractVector) where {T}
function _compositemul!(y::VecOut, A::CompositeMap{T,<:Tuple{LinearMap,LinearMap}}, x::AbstractVector, z::AbstractVector) where {T}
# no size checking, will be done by individual maps
A_mul_B!(z, A.maps[1], x)
A_mul_B!(y, A.maps[2], z)
return y
end
function _compositemul!(y::AbstractVector, A::CompositeMap{T,<:Tuple{Vararg{LinearMap}}}, x::AbstractVector, source::AbstractVector, dest::AbstractVector) where {T}
function _compositemul!(y::VecOut, A::CompositeMap{T,<:Tuple{Vararg{LinearMap}}}, x::AbstractVector, source::AbstractVector, dest::AbstractVector) where {T}
# no size checking, will be done by individual maps
N = length(A.maps)
A_mul_B!(source, A.maps[1], x)
Expand All @@ -166,6 +166,6 @@ function _compositemul!(y::AbstractVector, A::CompositeMap{T,<:Tuple{Vararg{Line
return y
end

At_mul_B!(y::AbstractVector, A::CompositeMap, x::AbstractVector) = A_mul_B!(y, transpose(A), x)
At_mul_B!(y::VecOut, A::CompositeMap, x::AbstractVector) = A_mul_B!(y, transpose(A), x)

Ac_mul_B!(y::AbstractVector, A::CompositeMap, x::AbstractVector) = A_mul_B!(y, adjoint(A), x)
Ac_mul_B!(y::VecOut, A::CompositeMap, x::AbstractVector) = A_mul_B!(y, adjoint(A), x)
6 changes: 3 additions & 3 deletions src/functionmap.jl
Original file line number Diff line number Diff line change
Expand Up @@ -107,13 +107,13 @@ function Base.:(*)(A::TransposeMap{<:Any,<:FunctionMap}, x::AbstractVector)
end
end

function A_mul_B!(y::AbstractVector, A::FunctionMap, x::AbstractVector)
function A_mul_B!(y::VecOut, A::FunctionMap, x::AbstractVector)
(length(x) == A.N && length(y) == A.M) || throw(DimensionMismatch("A_mul_B!"))
ismutating(A) ? A.f(y, x) : copyto!(y, A.f(x))
return y
end

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

function Ac_mul_B!(y::AbstractVector, A::FunctionMap, x::AbstractVector)
function Ac_mul_B!(y::VecOut, A::FunctionMap, x::AbstractVector)
ishermitian(A) && return A_mul_B!(y, A, x)
(length(x) == A.M && length(y) == A.N) || throw(DimensionMismatch("Ac_mul_B!"))
if A.fc !== nothing
Expand Down
18 changes: 9 additions & 9 deletions src/kronecker.jl
Original file line number Diff line number Diff line change
Expand Up @@ -131,15 +131,15 @@ end
# multiplication with vectors
#################

Base.@propagate_inbounds function A_mul_B!(y::AbstractVector, L::KroneckerMap{T,<:NTuple{2,LinearMap}}, x::AbstractVector) where {T}
Base.@propagate_inbounds function A_mul_B!(y::VecOut, L::KroneckerMap{T,<:NTuple{2,LinearMap}}, x::AbstractVector) where {T}
require_one_based_indexing(y)
@boundscheck check_dim_mul(y, L, x)
A, B = L.maps
X = LinearMap(reshape(x, (size(B, 2), size(A, 2))); issymmetric=false, ishermitian=false, isposdef=false)
_kronmul!(y, B, X, transpose(A), T)
return y
end
Base.@propagate_inbounds function A_mul_B!(y::AbstractVector, L::KroneckerMap{T}, x::AbstractVector) where {T}
Base.@propagate_inbounds function A_mul_B!(y::VecOut, L::KroneckerMap{T}, x::AbstractVector) where {T}
require_one_based_indexing(y)
@boundscheck check_dim_mul(y, L, x)
A = first(L.maps)
Expand All @@ -150,7 +150,7 @@ Base.@propagate_inbounds function A_mul_B!(y::AbstractVector, L::KroneckerMap{T}
end
# mixed-product rule, prefer the right if possible
# (A₁ ⊗ A₂ ⊗ ... ⊗ Aᵣ) * (B₁ ⊗ B₂ ⊗ ... ⊗ Bᵣ) = (A₁B₁) ⊗ (A₂B₂) ⊗ ... ⊗ (AᵣBᵣ)
Base.@propagate_inbounds function A_mul_B!(y::AbstractVector, L::CompositeMap{<:Any,<:Tuple{KroneckerMap,KroneckerMap}}, x::AbstractVector)
Base.@propagate_inbounds function A_mul_B!(y::VecOut, L::CompositeMap{<:Any,<:Tuple{KroneckerMap,KroneckerMap}}, x::AbstractVector)
B, A = L.maps
if length(A.maps) == length(B.maps) && all(M -> check_dim_mul(M[1], M[2]), zip(A.maps, B.maps))
A_mul_B!(y, kron(map(*, A.maps, B.maps)...), x)
Expand All @@ -160,7 +160,7 @@ Base.@propagate_inbounds function A_mul_B!(y::AbstractVector, L::CompositeMap{<:
end
# mixed-product rule, prefer the right if possible
# (A₁ ⊗ B₁)*(A₂⊗B₂)*...*(Aᵣ⊗Bᵣ) = (A₁*A₂*...*Aᵣ) ⊗ (B₁*B₂*...*Bᵣ)
Base.@propagate_inbounds function A_mul_B!(y::AbstractVector, L::CompositeMap{T,<:Tuple{Vararg{KroneckerMap{<:Any,<:Tuple{LinearMap,LinearMap}}}}}, x::AbstractVector) where {T}
Base.@propagate_inbounds function A_mul_B!(y::VecOut, L::CompositeMap{T,<:Tuple{Vararg{KroneckerMap{<:Any,<:Tuple{LinearMap,LinearMap}}}}}, x::AbstractVector) where {T}
As = map(AB -> AB.maps[1], L.maps)
Bs = map(AB -> AB.maps[2], L.maps)
As1, As2 = Base.front(As), Base.tail(As)
Expand All @@ -173,10 +173,10 @@ Base.@propagate_inbounds function A_mul_B!(y::AbstractVector, L::CompositeMap{T,
end
end

Base.@propagate_inbounds At_mul_B!(y::AbstractVector, A::KroneckerMap, x::AbstractVector) =
Base.@propagate_inbounds At_mul_B!(y::VecOut, A::KroneckerMap, x::AbstractVector) =
A_mul_B!(y, transpose(A), x)

Base.@propagate_inbounds Ac_mul_B!(y::AbstractVector, A::KroneckerMap, x::AbstractVector) =
Base.@propagate_inbounds Ac_mul_B!(y::VecOut, A::KroneckerMap, x::AbstractVector) =
A_mul_B!(y, adjoint(A), x)

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

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

Base.@propagate_inbounds function A_mul_B!(y::AbstractVector, L::KroneckerSumMap, x::AbstractVector)
Base.@propagate_inbounds function A_mul_B!(y::VecOut, L::KroneckerSumMap, x::AbstractVector)
@boundscheck check_dim_mul(y, L, x)
A, B = L.maps
ma, na = size(A)
Expand All @@ -273,8 +273,8 @@ Base.@propagate_inbounds function A_mul_B!(y::AbstractVector, L::KroneckerSumMap
return y
end

Base.@propagate_inbounds At_mul_B!(y::AbstractVector, A::KroneckerSumMap, x::AbstractVector) =
Base.@propagate_inbounds At_mul_B!(y::VecOut, A::KroneckerSumMap, x::AbstractVector) =
A_mul_B!(y, transpose(A), x)

Base.@propagate_inbounds Ac_mul_B!(y::AbstractVector, A::KroneckerSumMap, x::AbstractVector) =
Base.@propagate_inbounds Ac_mul_B!(y::VecOut, A::KroneckerSumMap, x::AbstractVector) =
A_mul_B!(y, adjoint(A), x)
35 changes: 35 additions & 0 deletions src/left.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,35 @@
# left.jl
# "special cases" related left vector multiplication like x = y'*A
# The subtlety is that y' is a AdjointAbsVec or a TransposeAbsVec
# which is a subtype of AbstractMatrix but it is really "vector like"
# so we want handle it like a vector (Issue#99)
# So this is an exception to the left multiplication rule by a AbstractMatrix
# that usually makes a WrappedMap.
# The "transpose" versions may be of dubious use, but the adjoint versions
# are useful for ensuring that (y'A)*x ≈ y'*(A*x) are both scalars.

import LinearAlgebra: AdjointAbsVec, TransposeAbsVec

Base.:(*)(y::AdjointAbsVec, A::LinearMap) = adjoint(*(A', y'))
Base.:(*)(y::TransposeAbsVec, A::LinearMap) = transpose(transpose(A) * transpose(y))

# mul!(x, y', A)
LinearAlgebra.mul!(x::AdjointAbsVec, y::AdjointAbsVec, A::LinearMap) =
mul!(x, y, A, true, false)

# not sure if we need bounds checks and propagate inbounds stuff here

# mul!(x, y', A, α, β)
function LinearAlgebra.mul!(x::AdjointAbsVec, y::AdjointAbsVec, A::LinearMap,
α::Number, β::Number)
check_dim_mul(x, y, A)
mul!(x, A', y', α, β)
return conj!(x)
end

# mul!(x, transpose(y), A, α, β)
function LinearAlgebra.mul!(x::TransposeAbsVec, y::TransposeAbsVec, A::LinearMap,
α::Number, β::Number)
check_dim_mul(x, y, A)
return mul!(x, transpose(A), transpose(y), α, β)
end
6 changes: 3 additions & 3 deletions src/linearcombination.jl
Original file line number Diff line number Diff line change
Expand Up @@ -117,8 +117,8 @@ end
return y
end

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

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

Ac_mul_B!(y::AbstractVector, A::LinearCombination, x::AbstractVector) = mul!(y, adjoint(A), x, true, false)
Ac_mul_B!(y::VecOut, A::LinearCombination, x::AbstractVector) = mul!(y, adjoint(A), x, true, false)
8 changes: 4 additions & 4 deletions src/scaledmap.jl
Original file line number Diff line number Diff line change
Expand Up @@ -54,15 +54,15 @@ Base.:(*)(A::ScaledMap, B::LinearMap) = A.λ * (A.lmap * B)
Base.:(*)(A::LinearMap, B::ScaledMap) = (A * B.lmap) * B.λ

# multiplication with vectors
function A_mul_B!(y::AbstractVector, A::ScaledMap, x::AbstractVector)
function A_mul_B!(y::VecOut, A::ScaledMap, x::AbstractVector)
# no size checking, will be done by map
mul!(y, A.lmap, x, A.λ, false)
end

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

At_mul_B!(y::AbstractVector, A::ScaledMap, x::AbstractVector) = A_mul_B!(y, transpose(A), x)
Ac_mul_B!(y::AbstractVector, A::ScaledMap, x::AbstractVector) = A_mul_B!(y, adjoint(A), x)
At_mul_B!(y::VecOut, A::ScaledMap, x::AbstractVector) = A_mul_B!(y, transpose(A), x)
Ac_mul_B!(y::VecOut, A::ScaledMap, x::AbstractVector) = A_mul_B!(y, adjoint(A), x)
12 changes: 6 additions & 6 deletions src/transpose.jl
Original file line number Diff line number Diff line change
Expand Up @@ -32,18 +32,18 @@ Base.:(==)(A::LinearMap, B::TransposeMap) = issymmetric(A) && B.lmap == A
Base.:(==)(A::LinearMap, B::AdjointMap) = ishermitian(A) && B.lmap == A

# multiplication with vector
A_mul_B!(y::AbstractVector, A::TransposeMap, x::AbstractVector) =
A_mul_B!(y::VecOut, A::TransposeMap, x::AbstractVector) =
issymmetric(A.lmap) ? A_mul_B!(y, A.lmap, x) : At_mul_B!(y, A.lmap, x)

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

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

A_mul_B!(y::AbstractVector, A::AdjointMap, x::AbstractVector) =
A_mul_B!(y::VecOut, A::AdjointMap, x::AbstractVector) =
ishermitian(A.lmap) ? A_mul_B!(y, A.lmap, x) : Ac_mul_B!(y, A.lmap, x)

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

Ac_mul_B!(y::AbstractVector, A::AdjointMap, x::AbstractVector) = A_mul_B!(y, A.lmap, x)
Ac_mul_B!(y::VecOut, A::AdjointMap, x::AbstractVector) = A_mul_B!(y, A.lmap, x)
6 changes: 3 additions & 3 deletions src/wrappedmap.jl
Original file line number Diff line number Diff line change
Expand Up @@ -37,13 +37,13 @@ LinearAlgebra.ishermitian(A::WrappedMap) = A._ishermitian
LinearAlgebra.isposdef(A::WrappedMap) = A._isposdef

# multiplication with vectors & matrices
A_mul_B!(y::AbstractVector, A::WrappedMap, x::AbstractVector) = A_mul_B!(y, A.lmap, x)
A_mul_B!(y::VecOut, A::WrappedMap, x::AbstractVector) = A_mul_B!(y, A.lmap, x)
Base.:(*)(A::WrappedMap, x::AbstractVector) = *(A.lmap, x)

At_mul_B!(y::AbstractVector, A::WrappedMap, x::AbstractVector) =
At_mul_B!(y::VecOut, A::WrappedMap, x::AbstractVector) =
(issymmetric(A) || (isreal(A) && ishermitian(A))) ? A_mul_B!(y, A.lmap, x) : At_mul_B!(y, A.lmap, x)

Ac_mul_B!(y::AbstractVector, A::WrappedMap, x::AbstractVector) =
Ac_mul_B!(y::VecOut, A::WrappedMap, x::AbstractVector) =
ishermitian(A) ? A_mul_B!(y, A.lmap, x) : Ac_mul_B!(y, A.lmap, x)

if VERSION ≥ v"1.3.0-alpha.115"
Expand Down
Loading