Restrict ScaledMap to RealOrComplex (#127)

Author: dkarrasch

src/scaledmap.jl

@@ -1,15 +1,20 @@
-struct ScaledMap{T, S<:RealOrComplex, A<:LinearMap} <: LinearMap{T}
+"""
+Lazy representation of a scaled map `λ * A = A * λ` with real or complex map
+`A <: LinearMap{RealOrComplex}` and real or complex scaling factor
+`λ <: RealOrComplex`.
+"""
+struct ScaledMap{T, S<:RealOrComplex, L<:LinearMap} <: LinearMap{T}
     λ::S
-    lmap::A
-    function ScaledMap{T}(λ::S, lmap::A) where {T, S <: RealOrComplex, A <: LinearMap}
-        @assert Base.promote_op(*, S, eltype(lmap)) == T "target type $T cannot hold products of $S and $(eltype(lmap)) objects"
-        new{T,S,A}(λ, lmap)
+    lmap::L
+    function ScaledMap{T}(λ::S, A::L) where {T, S <: RealOrComplex, L <: LinearMap{<:RealOrComplex}}
+        @assert Base.promote_op(*, S, eltype(A)) == T "target type $T cannot hold products of $S and $(eltype(A)) objects"
+        new{T,S,L}(λ, A)
     end
 end
 
 # constructor
-ScaledMap(λ::S, lmap::A) where {S<:RealOrComplex,A<:LinearMap} =
-    ScaledMap{Base.promote_op(*, S, eltype(lmap))}(λ, lmap)
+ScaledMap(λ::RealOrComplex, lmap::LinearMap{<:RealOrComplex}) =
+    ScaledMap{Base.promote_op(*, typeof(λ), eltype(lmap))}(λ, lmap)
 
 # basic methods
 Base.size(A::ScaledMap) = size(A.lmap)
@@ -26,8 +31,8 @@ Base.:(==)(A::ScaledMap, B::ScaledMap) =
     eltype(A) == eltype(B) && A.lmap == B.lmap && A.λ == B.λ
 
 # scalar multiplication and division
-Base.:(*)(α::RealOrComplex, A::LinearMap) = ScaledMap(α, A)
-Base.:(*)(A::LinearMap, α::RealOrComplex) = ScaledMap(α, A)
+Base.:(*)(α::RealOrComplex, A::LinearMap{<:RealOrComplex}) = ScaledMap(α, A)
+Base.:(*)(A::LinearMap{<:RealOrComplex}, α::RealOrComplex) = ScaledMap(α, A)
 
 Base.:(*)(α::Number, A::ScaledMap) = (α * A.λ) * A.lmap
 Base.:(*)(A::ScaledMap, α::Number) = A.lmap * (A.λ * α)
@@ -42,13 +47,19 @@ Base.:(*)(A::ScaledMap, B::LinearMap) = A.λ * (A.lmap * B)
 Base.:(*)(A::LinearMap, B::ScaledMap) = (A * B.lmap) * B.λ
 
 # multiplication with vectors/matrices
-for (In, Out) in ((AbstractVector, AbstractVecOrMat), (AbstractMatrix, AbstractMatrix))
+for (In, Out) in ((AbstractVector, AbstractVecOrMat),
+                  (AbstractMatrix, AbstractMatrix))
     @eval begin
-        function _unsafe_mul!(y::$Out, A::ScaledMap, x::$In)
+        # commutative case
+        function _unsafe_mul!(y::$Out, A::ScaledMap, x::$In{<:RealOrComplex})
             return _unsafe_mul!(y, A.lmap, x, A.λ, false)
         end
-        function _unsafe_mul!(y::$Out, A::ScaledMap, x::$In, α::Number, β::Number)
+        function _unsafe_mul!(y::$Out, A::ScaledMap, x::$In{<:RealOrComplex}, α::Number, β::Number)
             return _unsafe_mul!(y, A.lmap, x, A.λ * α, β)
         end
+        # non-commutative case
+        function _unsafe_mul!(y::$Out, A::ScaledMap, x::$In)
+            return lmul!(A.λ, _unsafe_mul!(y, A.lmap, x))
+        end
     end
 end

test/numbertypes.jl

@@ -1,12 +1,18 @@
 using Test, LinearMaps, LinearAlgebra, Quaternions
 
+# type piracy because Quaternions.jl doesn't have it right
+Base.:(*)(z::Complex{T}, q::Quaternion{T}) where {T<:Real} = quat(z) * q
+Base.:(*)(q::Quaternion{T}, z::Complex{T}) where {T<:Real} = q * quat(z)
+
 @testset "noncommutative number type" begin
     x = Quaternion.(rand(10), rand(10), rand(10), rand(10))
     v = rand(10)
     A = Quaternion.(rand(10,10), rand(10,10), rand(10,10), rand(10,10))
+    B = rand(ComplexF64, 10, 10)
     γ = Quaternion.(rand(4)...) # "Number"
     α = UniformScaling(γ)
     β = UniformScaling(Quaternion.(rand(4)...))
+    λ = rand(ComplexF64)
     L = LinearMap(A)
     @test Array(L) == A
     @test Array(L') == A'
@@ -21,8 +27,14 @@ using Test, LinearMaps, LinearAlgebra, Quaternions
     @test L' * x ≈ A' * x
     @test α * (L * x) ≈ α * (A * x)
     @test α * L * x ≈ α * A * x
+    @test L * α * x ≈ A * α * x
     @test 3L * x ≈ 3A * x
     @test 3L' * x ≈ 3A' * x
+    @test λ*L isa LinearMaps.CompositeMap
+    @test γ * (λ * LinearMap(B)) isa LinearMaps.CompositeMap
+    @test (λ * LinearMap(B)) * γ isa LinearMaps.CompositeMap
+    @test λ*L * x ≈ λ*A * x
+    @test λ*L' * x ≈ λ*A' * x
     @test α * (3L * x) ≈ α * (3A * x)
     @test (@inferred α * 3L) * x ≈ α * 3A * x
     @test (@inferred 3L * α) * x ≈ 3A * α * x
@@ -53,6 +65,13 @@ using Test, LinearMaps, LinearAlgebra, Quaternions
     @test Array(-L) == -A
     @test Array(γ \ L) ≈ γ \ A
     @test Array(L / γ) ≈ A / γ
+    M = rand(ComplexF64, 10, 10); α = rand(ComplexF64);
+    y = α * M * x; Y = α * M * A
+    @test (α * LinearMap(M)) * x ≈ (quat(α) * LinearMap(M)) * x ≈ y
+    @test mul!(copy(y), α * LinearMap(M), x, α, false) ≈ α * M * x * α
+    @test mul!(copy(y), α * LinearMap(M), x, quat(α), false) ≈ α * M * x * α
+    @test mul!(copy(Y), α * LinearMap(M), A) ≈ α * M * A
+    @test mul!(copy(Y), α * LinearMap(M), A, α, false) ≈ α * M * A * α
 end
 
 @testset "nonassociative number type" begin