allow direct multiplication with out-of-place transpose/adjoint FunctionMaps (#83)

Author: dkarrasch

.travis.yml

@@ -8,7 +8,7 @@ os:
 
 julia:
   - 1.0
-  - 1.3
+  - 1
   - nightly
 
 matrix:

Project.toml

@@ -1,6 +1,6 @@
 name = "LinearMaps"
 uuid = "7a12625a-238d-50fd-b39a-03d52299707e"
-version = "2.6"
+version = "2.6.1"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

src/functionmap.jl

@@ -48,6 +48,64 @@ function Base.:(*)(A::FunctionMap, x::AbstractVector)
     end
     return y
 end
+function Base.:(*)(A::AdjointMap{<:Any,<:FunctionMap}, x::AbstractVector)
+    Afun = A.lmap
+    ishermitian(Afun) && return Afun*x
+    length(x) == size(A, 2) || throw(DimensionMismatch())
+    if Afun.fc !== nothing
+        if ismutating(Afun)
+            y = similar(x, promote_type(eltype(A), eltype(x)), size(A, 1))
+            Afun.fc(y, x)
+        else
+            y = Afun.fc(x)
+        end
+        return y
+    elseif issymmetric(Afun) # but !isreal(A), Afun.f can be used
+        x = conj(x)
+        if ismutating(Afun)
+            y = similar(x, promote_type(eltype(A), eltype(x)), size(A, 1))
+            Afun.f(y, x)
+        else
+            y = Afun.f(x)
+        end
+        conj!(y)
+        return y
+    else
+        error("adjoint not implemented for $A")
+    end
+end
+function Base.:(*)(A::TransposeMap{<:Any,<:FunctionMap}, x::AbstractVector)
+    Afun = A.lmap
+    (issymmetric(Afun) || (isreal(A) && ishermitian(Afun))) && return Afun*x
+    length(x) == size(A, 2) || throw(DimensionMismatch())
+    if Afun.fc !== nothing
+        if !isreal(A)
+            x = conj(x)
+        end
+        if ismutating(Afun)
+            y = similar(x, promote_type(eltype(A), eltype(x)), size(A, 1))
+            Afun.fc(y, x)
+        else
+            y = Afun.fc(x)
+        end
+        if !isreal(A)
+            conj!(y)
+        end
+        return y
+    elseif ishermitian(Afun) # but !isreal(A), Afun.f can be used
+        x = conj(x)
+        if ismutating(Afun)
+            y = similar(x, promote_type(eltype(A), eltype(x)), size(A, 1))
+            Afun.f(y, x)
+        else
+            y = Afun.f(x)
+        end
+        conj!(y)
+        return y
+    else
+        error("transpose not implemented for $A")
+    end
+end
 
 function A_mul_B!(y::AbstractVector, A::FunctionMap, x::AbstractVector)
     (length(x) == A.N && length(y) == A.M) || throw(DimensionMismatch("A_mul_B!"))
@@ -78,7 +136,7 @@ function At_mul_B!(y::AbstractVector, A::FunctionMap, x::AbstractVector)
 end
 
 function Ac_mul_B!(y::AbstractVector, A::FunctionMap, x::AbstractVector)
-    (ishermitian(A) || (isreal(A) && issymmetric(A))) && return A_mul_B!(y, A, x)
+    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
         ismutating(A) ? A.fc(y, x) : copyto!(y, A.fc(x))

test/functionmap.jl

@@ -56,6 +56,8 @@ using Test, LinearMaps, LinearAlgebra
     @test @inferred mul!(similar(v), CS!, v) == cv
     @test_throws ErrorException CS!'v
     @test_throws ErrorException adjoint(CS!) * v
+    @test_throws ErrorException mul!(similar(v), CS!', v)
+    @test_throws ErrorException mul!(similar(v), transpose(CS!), v)
     CS! = LinearMap{ComplexF64}(cumsum!, (y, x) -> (copyto!(y, x); reverse!(y); cumsum!(y, y)), 10; ismutating=true)
     @inferred adjoint(CS!)
     @test @inferred LinearMaps.ismutating(CS!)
@@ -73,15 +75,37 @@ using Test, LinearMaps, LinearAlgebra
     @test @inferred transpose(2 * L) * v ≈ 2 * v
     L = @inferred LinearMap{ComplexF64}(x -> x, x -> x, 10)
     v = rand(ComplexF64, 10)
+    w = similar(v)
     @test @inferred (2 * L)' * v ≈ 2 * v
     @test @inferred transpose(2 * L) * v ≈ 2 * v
 
-    L = LinearMap{ComplexF64}(identity, 10; issymmetric=true)
-    @test L * v == v
-    @test adjoint(L) * v == v
-    @test transpose(L) * v == v
-    L = LinearMap{ComplexF64}(identity, 10; ishermitian=true)
-    @test L * v == v
-    @test adjoint(L) * v == v
-    @test transpose(L) * v == v
+    A = rand(ComplexF64, 10, 10)
+    L = LinearMap{ComplexF64}(x -> A*x, 10)
+    @test L * v == A * v == mul!(w, L, v)
+    L = LinearMap{ComplexF64}((y, x) -> mul!(y, A, x), 10)
+    @test L * v == A * v == mul!(w, L, v)
+    L = LinearMap{ComplexF64}((y, x) -> mul!(y, A, x), (y, x) -> mul!(y, A', x), 10)
+    @test L * v == A * v == mul!(w, L, v)
+    @test adjoint(L) * v ≈ A'v ≈ mul!(w, L', v)
+    @test transpose(L) * v ≈ transpose(A)*v ≈ mul!(w, transpose(L), v)
+
+    A = Symmetric(rand(ComplexF64, 10, 10))
+    L = LinearMap{ComplexF64}(x -> A*x, 10; issymmetric=true)
+    @test L * v == A * v == mul!(w, L, v)
+    @test adjoint(L) * v ≈ A'v ≈ mul!(w, L', v)
+    @test transpose(L) * v ≈ transpose(A)*v ≈ mul!(w, transpose(L), v)
+    L = LinearMap{ComplexF64}((y, x) -> mul!(y, A, x), 10; issymmetric=true)
+    @test L * v == A * v == mul!(w, L, v)
+    @test adjoint(L) * v ≈ A'v ≈ mul!(w, L', v)
+    @test transpose(L) * v ≈ transpose(A)*v ≈ mul!(w, transpose(L), v)
+
+    A = Hermitian(rand(ComplexF64, 10, 10))
+    L = LinearMap{ComplexF64}(x -> A*x, 10; ishermitian=true)
+    @test L * v == A * v == mul!(w, L, v)
+    @test adjoint(L) * v ≈ A'v ≈ mul!(w, L', v)
+    @test transpose(L) * v ≈ transpose(A)*v ≈ mul!(w, transpose(L), v)
+    L = LinearMap{ComplexF64}((y, x) -> mul!(y, A, x), 10; ishermitian=true)
+    @test L * v == A * v == mul!(w, L, v)
+    @test adjoint(L) * v ≈ A'v ≈ mul!(w, L', v)
+    @test transpose(L) * v ≈ transpose(A)*v ≈ mul!(w, transpose(L), v)
 end