Improve specific kron case (#139)

Author: dkarrasch

Project.toml

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

src/kronecker.jl

@@ -141,7 +141,12 @@ end
     mb, nb = size(B)
     X = reshape(x, (nb, na))
     Y = reshape(y, (mb, ma))
-    if nb*ma < mb*na
+    if B isa UniformScalingMap
+        # the following is (maybe due to the reshape?) faster than
+        # _unsafe_mul!(Y, B * X, At.lmap)
+        _unsafe_mul!(Y, X, At.lmap)
+        lmul!(B.λ, y)
+    elseif nb*ma <= mb*na
         _unsafe_mul!(Y, B, X * At.lmap)
     else
         _unsafe_mul!(Y, Matrix(B*X), At.lmap)

src/uniformscalingmap.jl

@@ -39,8 +39,16 @@ Base.:(*)(α::RealOrComplex, J::UniformScalingMap) = UniformScalingMap(α * J.λ
 Base.:(*)(J::UniformScalingMap, α::RealOrComplex) = UniformScalingMap(J.λ * α, size(J))
 
 # multiplication with vector
-Base.:(*)(A::UniformScalingMap, x::AbstractVector) =
-    length(x) == A.M ? A.λ * x : throw(DimensionMismatch("*"))
+Base.:(*)(J::UniformScalingMap, x::AbstractVector) =
+    length(x) == J.M ? J.λ * x : throw(DimensionMismatch("*"))
+# multiplication with matrix
+Base.:(*)(J::UniformScalingMap, B::AbstractMatrix) =
+    size(B, 1) == J.M ? J.λ * LinearMap(B) : throw(DimensionMismatch("*"))
+Base.:(*)(A::AbstractMatrix, J::UniformScalingMap) =
+    size(A, 2) == J.M ? LinearMap(A) * J.λ : throw(DimensionMismatch("*"))
+# disambiguation
+Base.:(*)(xc::LinearAlgebra.AdjointAbsVec, J::UniformScalingMap) = xc * J.λ
+Base.:(*)(xt::LinearAlgebra.TransposeAbsVec, J::UniformScalingMap) = xt * J.λ
 
 # multiplication with vector/matrix
 for (In, Out) in ((AbstractVector, AbstractVecOrMat), (AbstractMatrix, AbstractMatrix))
@@ -49,6 +57,11 @@ for (In, Out) in ((AbstractVector, AbstractVecOrMat), (AbstractMatrix, AbstractM
             _scaling!(y, J.λ, x, true, false)
             return y
         end
+        function _unsafe_mul!(y::$Out, J::UniformScalingMap{<:RealOrComplex}, x::$In{<:RealOrComplex},
+                    α::RealOrComplex, β::Number)
+            _scaling!(y, J.λ * α, x, true, β)
+            return y
+        end
         function _unsafe_mul!(y::$Out, J::UniformScalingMap, x::$In,
                     α::Number, β::Number)
             _scaling!(y, J.λ, x, α, β)

test/kronecker.jl

@@ -66,6 +66,18 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays
         @test Matrix(@inferred K*K) ≈ kron(A, B)*kron(A, B)
         A = rand(3, 2); B = rand(4, 3)
         @test Matrix(kron(LinearMap(A), B, [A A])*kron(LinearMap(A), B, A')) ≈ kron(A, B, [A A])*kron(A, B, A')
+
+        m = 3
+        A = rand(m, m)
+        S = sparse(I, m, m)
+        J = LinearMap(I, m)
+        v = rand(m^3)
+        for (K, M) in ((⊗(A, J, J), kron(A, S, S)),
+                       (⊗(J, A, J), kron(S, A, S)),
+                       (⊗(J, J, A), kron(S, S, A)))
+            @test K * v ≈ M * v
+            @test Matrix(K) ≈ M
+        end
     end
 
     @testset "Kronecker sum" begin

test/left.jl

@@ -70,4 +70,7 @@ end
     @test left_tester(W)
     @test left_tester(W')
     @test left_tester(transpose(W))
+    
+    J = LinearMap(1.0I, 5) # UniformScalingMap
+    @test left_tester(J)
 end

test/numbertypes.jl

@@ -35,30 +35,30 @@ Base.:(+)(q::Quaternion, z::Complex) = q + quat(z)
         @test mul!(copy(C), transpose(A), M, γ, λ) ≈ transpose(A)*A*γ + C*λ
         @test mul!(copy(C), adjoint(A), M, γ, λ) ≈ A'*A*γ + C*λ
     end
-    @test Array((α * F')') ≈ (α * A')' ≈ A * conj(α)
+    @test Array((α * F')') ≈ (γ * A')' ≈ A * conj(γ)
     @test L * x ≈ A * x
     @test L' * x ≈ A' * x
-    @test α * (L * x) ≈ α * (A * x)
-    @test α * L * x ≈ α * A * x
-    @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
-    @test (α * L') * x ≈ (α * A') * x
-    @test (α * L')' * x ≈ (α * A')' * x
-    @test (α * L')' * v ≈ (α * A')' * v
-    @test Array(@inferred adjoint(α * L * β)) ≈ conj(β) * A' * conj(α)
-    @test Array(@inferred transpose(α * L * β)) ≈ β * transpose(A) * α
+    @test α * (3L * x) ≈ γ * (3A * x)
+    @test (@inferred α * 3L) * x ≈ γ * 3A * x
+    @test (@inferred 3L * α) * x ≈ 3A * γ * x
+    @test (α * L') * x ≈ (γ * A') * x
+    @test (α * L')' * x ≈ (γ * A')' * x
+    @test (α * L')' * v ≈ (γ * A')' * v
+    @test Array(@inferred adjoint(α * L * β)) ≈ conj(β) * A' * conj(γ)
+    @test Array(@inferred transpose(α * L * β)) ≈ β * transpose(A) * γ
     J = LinearMap(α, 10)
-    @test (β * J) * x ≈ LinearMap(β*α, 10) * x ≈ β*α*x
-    @test (J * β) * x ≈ LinearMap(α*β, 10) * x ≈ α*β*x
+    @test (β * J) * x ≈ LinearMap(β*α, 10) * x ≈ β*γ*x
+    @test (J * β) * x ≈ LinearMap(α*β, 10) * x ≈ γ*β*x
     M = β.λ * (γ * L * L)
     @test M == β * (γ * L * L)
     @test length(M.maps) == 3
@@ -71,7 +71,11 @@ Base.:(+)(q::Quaternion, z::Complex) = q + quat(z)
     @test length(M.maps) == 3
     @test M.maps[1].λ == γ*β.λ
     @test γ*FillMap(γ, (3, 4)) == FillMap(γ^2, (3, 4)) == FillMap(γ, (3, 4))*γ
-
+    U = LinearMap(quat(1.0)*I, 10)
+    for β in (0, 1, rand())
+        @test mul!(copy(x), J, x, γ, β) == γ * x * γ + x * β
+        @test mul!(copy(x), U, x, γ, β) == x * γ + x * β
+    end
     # exercise non-RealOrComplex scalar operations
     @test Array(γ * (L'*L)) ≈ γ * (A'*A) # CompositeMap
     @test Array((L'*L) * γ) ≈ (A'*A) * γ

test/uniformscalingmap.jl

@@ -48,4 +48,8 @@ using Test, LinearMaps, LinearAlgebra, BenchmarkTools
     end
     X = rand(10, 10); Y = similar(X)
     @test mul!(Y, Id, X) == X
+    @test Id*X isa LinearMap
+    @test X*Id isa LinearMap
+    @test Matrix(Id*X) == X
+    @test Matrix(X*Id) == X
 end