From e59019597ff1998883b6c13add6b2a81c357b7b3 Mon Sep 17 00:00:00 2001 From: Elias Jarlebring Date: Tue, 4 May 2021 21:18:50 +0200 Subject: [PATCH 01/10] first mul! impl for uniformscaling --- stdlib/LinearAlgebra/src/uniformscaling.jl | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/stdlib/LinearAlgebra/src/uniformscaling.jl b/stdlib/LinearAlgebra/src/uniformscaling.jl index 1b2e3010ff27d..9e624d88c50f9 100644 --- a/stdlib/LinearAlgebra/src/uniformscaling.jl +++ b/stdlib/LinearAlgebra/src/uniformscaling.jl @@ -284,6 +284,22 @@ end mul!(C, A, J.λ, alpha, beta) @inline mul!(C::AbstractVecOrMat, J::UniformScaling, B::AbstractVecOrMat, alpha::Number, beta::Number) = mul!(C, J.λ, B, alpha, beta) + +function mul!(out::AbstractMatrix{T},a::Bool,B::UniformScaling{T},α::Bool,β::Bool) where {T} + if !β # zero contribution of the out matrix + fill!(out,zero(T)); + end + if (a*α) + # Add B.λ*I + m = min(size(out)...) + s = B.λ + @inbounds for i = 1:m; + out[i,i] += s; + end + end + return out +end + rmul!(A::AbstractMatrix, J::UniformScaling) = rmul!(A, J.λ) lmul!(J::UniformScaling, B::AbstractVecOrMat) = lmul!(J.λ, B) rdiv!(A::AbstractMatrix, J::UniformScaling) = rdiv!(A, J.λ) From bb76b95985aacc79b75520d9d67bf395e67e10f8 Mon Sep 17 00:00:00 2001 From: Elias Jarlebring Date: Thu, 6 May 2021 12:43:35 +0200 Subject: [PATCH 02/10] use five arg mul! for in-place I adds --- stdlib/LinearAlgebra/src/dense.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/stdlib/LinearAlgebra/src/dense.jl b/stdlib/LinearAlgebra/src/dense.jl index 94926805bb387..8020d86c8c4a1 100644 --- a/stdlib/LinearAlgebra/src/dense.jl +++ b/stdlib/LinearAlgebra/src/dense.jl @@ -636,8 +636,8 @@ function exp!(A::StridedMatrix{T}) where T<:BlasFloat # Compute U and V: Even/odd terms in Padé numerator & denom # Expansion of k=1 in for loop P = A2 - U = C[2]*I + C[4]*P - V = C[1]*I + C[3]*P + U = mul!(C[4]*P, true, C[2]*I, true, true) #U = C[2]*I + C[4]*P + V = mul!(C[3]*P, true, C[1]*I, true, true) #V = C[1]*I + C[3]*P for k in 2:(div(size(C, 1), 2) - 1) k2 = 2 * k P *= A2 @@ -664,7 +664,7 @@ function exp!(A::StridedMatrix{T}) where T<:BlasFloat A4 = A2 * A2 A6 = A2 * A4 Ut = CC[4]*A2 - Ut[diagind(Ut)] .+= CC[2] + Ut = mul!(CC[4]*A2,true,CC[2]*I,true,true); # Ut = CC[4]*A2+CC[2]*I # Allocation economical version of: #U = A * (A6 * (CC[14].*A6 .+ CC[12].*A4 .+ CC[10].*A2) .+ # CC[8].*A6 .+ CC[6].*A4 .+ Ut) @@ -676,7 +676,7 @@ function exp!(A::StridedMatrix{T}) where T<:BlasFloat # Allocation economical version of: Vt = CC[3]*A2 (recycle Ut) Vt = mul!(Ut, CC[3], A2, true, false) - Vt[diagind(Vt)] .+= CC[1] + mul!(Vt,true,CC[1]*I,true,true); # Vt += CC[1]*I # Allocation economical version of: #V = A6 * (CC[13].*A6 .+ CC[11].*A4 .+ CC[9].*A2) .+ # CC[7].*A6 .+ CC[5].*A4 .+ Vt From 3d127dce3fbad4cfce60992c1e67c5ebec20cfb9 Mon Sep 17 00:00:00 2001 From: Elias Jarlebring Date: Thu, 6 May 2021 12:51:55 +0200 Subject: [PATCH 03/10] Typos --- stdlib/LinearAlgebra/src/dense.jl | 5 ++--- stdlib/LinearAlgebra/src/uniformscaling.jl | 2 +- 2 files changed, 3 insertions(+), 4 deletions(-) diff --git a/stdlib/LinearAlgebra/src/dense.jl b/stdlib/LinearAlgebra/src/dense.jl index 8020d86c8c4a1..e0fd0a6d532bc 100644 --- a/stdlib/LinearAlgebra/src/dense.jl +++ b/stdlib/LinearAlgebra/src/dense.jl @@ -663,8 +663,7 @@ function exp!(A::StridedMatrix{T}) where T<:BlasFloat A2 = A * A A4 = A2 * A2 A6 = A2 * A4 - Ut = CC[4]*A2 - Ut = mul!(CC[4]*A2,true,CC[2]*I,true,true); # Ut = CC[4]*A2+CC[2]*I + Ut = mul!(CC[4]*A2, true,CC[2]*I, true, true); # Ut = CC[4]*A2+CC[2]*I # Allocation economical version of: #U = A * (A6 * (CC[14].*A6 .+ CC[12].*A4 .+ CC[10].*A2) .+ # CC[8].*A6 .+ CC[6].*A4 .+ Ut) @@ -676,7 +675,7 @@ function exp!(A::StridedMatrix{T}) where T<:BlasFloat # Allocation economical version of: Vt = CC[3]*A2 (recycle Ut) Vt = mul!(Ut, CC[3], A2, true, false) - mul!(Vt,true,CC[1]*I,true,true); # Vt += CC[1]*I + mul!(Vt, true, CC[1]*I, true, true); # Vt += CC[1]*I # Allocation economical version of: #V = A6 * (CC[13].*A6 .+ CC[11].*A4 .+ CC[9].*A2) .+ # CC[7].*A6 .+ CC[5].*A4 .+ Vt diff --git a/stdlib/LinearAlgebra/src/uniformscaling.jl b/stdlib/LinearAlgebra/src/uniformscaling.jl index 9e624d88c50f9..4fa2f0d8b1aeb 100644 --- a/stdlib/LinearAlgebra/src/uniformscaling.jl +++ b/stdlib/LinearAlgebra/src/uniformscaling.jl @@ -285,7 +285,7 @@ end @inline mul!(C::AbstractVecOrMat, J::UniformScaling, B::AbstractVecOrMat, alpha::Number, beta::Number) = mul!(C, J.λ, B, alpha, beta) -function mul!(out::AbstractMatrix{T},a::Bool,B::UniformScaling{T},α::Bool,β::Bool) where {T} +function mul!(out::AbstractMatrix{T}, a::Bool, B::UniformScaling{T}, α::Bool, β::Bool) where {T} if !β # zero contribution of the out matrix fill!(out,zero(T)); end From 6c06c63a3049d00436958b168c977df1e4c17c32 Mon Sep 17 00:00:00 2001 From: Elias Jarlebring Date: Thu, 6 May 2021 17:44:27 +0200 Subject: [PATCH 04/10] More general definition bool -> number Co-authored-by: Daniel Karrasch --- stdlib/LinearAlgebra/src/uniformscaling.jl | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/stdlib/LinearAlgebra/src/uniformscaling.jl b/stdlib/LinearAlgebra/src/uniformscaling.jl index 4fa2f0d8b1aeb..663049c767ede 100644 --- a/stdlib/LinearAlgebra/src/uniformscaling.jl +++ b/stdlib/LinearAlgebra/src/uniformscaling.jl @@ -285,16 +285,17 @@ end @inline mul!(C::AbstractVecOrMat, J::UniformScaling, B::AbstractVecOrMat, alpha::Number, beta::Number) = mul!(C, J.λ, B, alpha, beta) -function mul!(out::AbstractMatrix{T}, a::Bool, B::UniformScaling{T}, α::Bool, β::Bool) where {T} - if !β # zero contribution of the out matrix - fill!(out,zero(T)); +function mul!(out::AbstractMatrix{T}, a::Number, B::UniformScaling, α::Number, β::Number) where {T} + checksquare(out) + if iszero(β) # zero contribution of the out matrix + fill!(out, zero(T)) + elseif !isone(β) + rmul!(out, β) end - if (a*α) - # Add B.λ*I - m = min(size(out)...) - s = B.λ - @inbounds for i = 1:m; - out[i,i] += s; + s = convert(T, a*B.λ*α) + if !iszero(s) + @inbounds for i in diagind(out) + out[i] += s end end return out From 02eb439d37d97d9316d04cfc25aa81f053d1ccc4 Mon Sep 17 00:00:00 2001 From: Elias Jarlebring Date: Thu, 6 May 2021 18:12:35 +0200 Subject: [PATCH 05/10] unit test for five arg mul! uniformscaling --- stdlib/LinearAlgebra/test/uniformscaling.jl | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/stdlib/LinearAlgebra/test/uniformscaling.jl b/stdlib/LinearAlgebra/test/uniformscaling.jl index 33b2ba7032734..068d002c18206 100644 --- a/stdlib/LinearAlgebra/test/uniformscaling.jl +++ b/stdlib/LinearAlgebra/test/uniformscaling.jl @@ -454,6 +454,11 @@ end target = J * A * alpha + C * beta @test mul!(copy(C), J, A, alpha, beta) ≈ target @test mul!(copy(C), A, J, alpha, beta) ≈ target + + a = randn(); + target_5mul = a*alpha*J+beta*C; + @test mul!(copy(C), a, J, alpha, beta) ≈ target_5mul + end @testset "Construct Diagonal from UniformScaling" begin From aab88fd0d772fae95c24b2f9300abea7cc62581a Mon Sep 17 00:00:00 2001 From: Elias Jarlebring Date: Thu, 6 May 2021 18:15:40 +0200 Subject: [PATCH 06/10] unit test for five arg mul! uniformscaling --- stdlib/LinearAlgebra/test/uniformscaling.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/stdlib/LinearAlgebra/test/uniformscaling.jl b/stdlib/LinearAlgebra/test/uniformscaling.jl index 068d002c18206..8fb3edadfc179 100644 --- a/stdlib/LinearAlgebra/test/uniformscaling.jl +++ b/stdlib/LinearAlgebra/test/uniformscaling.jl @@ -455,8 +455,9 @@ end @test mul!(copy(C), J, A, alpha, beta) ≈ target @test mul!(copy(C), A, J, alpha, beta) ≈ target - a = randn(); - target_5mul = a*alpha*J+beta*C; + a = randn() + C = randn(3, 3) + target_5mul = a*alpha*J + beta*C @test mul!(copy(C), a, J, alpha, beta) ≈ target_5mul end From 5cea2baabf889eabb42d186fb14d9a2ef4a1ef9f Mon Sep 17 00:00:00 2001 From: Elias Jarlebring Date: Thu, 6 May 2021 20:13:35 +0200 Subject: [PATCH 07/10] alpha=0, beta=0, unit test for five arg mul! uniformscaling --- stdlib/LinearAlgebra/test/uniformscaling.jl | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/stdlib/LinearAlgebra/test/uniformscaling.jl b/stdlib/LinearAlgebra/test/uniformscaling.jl index 8fb3edadfc179..b9639d75c1d6b 100644 --- a/stdlib/LinearAlgebra/test/uniformscaling.jl +++ b/stdlib/LinearAlgebra/test/uniformscaling.jl @@ -459,6 +459,10 @@ end C = randn(3, 3) target_5mul = a*alpha*J + beta*C @test mul!(copy(C), a, J, alpha, beta) ≈ target_5mul + target_5mul = beta*C # alpha = 0 + @test mul!(copy(C), a, J, 0, beta) ≈ target_5mul + target_5mul = a*alpha*Matrix(J,n,n) # beta = 0 + @test mul!(copy(C), a, J, alpha, 0) ≈ target_5mul end From bafaaf603f0697f3805e1dc9bf59c2bbaf45c3df Mon Sep 17 00:00:00 2001 From: Elias Jarlebring Date: Thu, 6 May 2021 20:38:10 +0200 Subject: [PATCH 08/10] flip argument 2 and 3 --- stdlib/LinearAlgebra/src/uniformscaling.jl | 3 ++- stdlib/LinearAlgebra/test/uniformscaling.jl | 1 + 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/stdlib/LinearAlgebra/src/uniformscaling.jl b/stdlib/LinearAlgebra/src/uniformscaling.jl index 663049c767ede..728d12f6359a4 100644 --- a/stdlib/LinearAlgebra/src/uniformscaling.jl +++ b/stdlib/LinearAlgebra/src/uniformscaling.jl @@ -300,7 +300,8 @@ function mul!(out::AbstractMatrix{T}, a::Number, B::UniformScaling, α::Number, end return out end - +@inline mul!(out::AbstractMatrix, A::UniformScaling, b::Number, α::Number, β::Number)= + mul!(out, b, A, α, β) rmul!(A::AbstractMatrix, J::UniformScaling) = rmul!(A, J.λ) lmul!(J::UniformScaling, B::AbstractVecOrMat) = lmul!(J.λ, B) rdiv!(A::AbstractMatrix, J::UniformScaling) = rdiv!(A, J.λ) diff --git a/stdlib/LinearAlgebra/test/uniformscaling.jl b/stdlib/LinearAlgebra/test/uniformscaling.jl index b9639d75c1d6b..aa2b3117ef456 100644 --- a/stdlib/LinearAlgebra/test/uniformscaling.jl +++ b/stdlib/LinearAlgebra/test/uniformscaling.jl @@ -459,6 +459,7 @@ end C = randn(3, 3) target_5mul = a*alpha*J + beta*C @test mul!(copy(C), a, J, alpha, beta) ≈ target_5mul + @test mul!(copy(C), J, a, alpha, beta) ≈ target_5mul target_5mul = beta*C # alpha = 0 @test mul!(copy(C), a, J, 0, beta) ≈ target_5mul target_5mul = a*alpha*Matrix(J,n,n) # beta = 0 From 78fb45bbb885fecac9cf27cbe169603a84b5c748 Mon Sep 17 00:00:00 2001 From: Elias Jarlebring Date: Fri, 7 May 2021 08:14:56 +0200 Subject: [PATCH 09/10] non-commutative respecting argument flipping --- stdlib/LinearAlgebra/src/uniformscaling.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stdlib/LinearAlgebra/src/uniformscaling.jl b/stdlib/LinearAlgebra/src/uniformscaling.jl index 728d12f6359a4..44ad5d98e99a4 100644 --- a/stdlib/LinearAlgebra/src/uniformscaling.jl +++ b/stdlib/LinearAlgebra/src/uniformscaling.jl @@ -301,7 +301,7 @@ function mul!(out::AbstractMatrix{T}, a::Number, B::UniformScaling, α::Number, return out end @inline mul!(out::AbstractMatrix, A::UniformScaling, b::Number, α::Number, β::Number)= - mul!(out, b, A, α, β) + mul!(out, A.λ, UniformScaling(b), α, β) rmul!(A::AbstractMatrix, J::UniformScaling) = rmul!(A, J.λ) lmul!(J::UniformScaling, B::AbstractVecOrMat) = lmul!(J.λ, B) rdiv!(A::AbstractMatrix, J::UniformScaling) = rdiv!(A, J.λ) From ab1b1d45c2ce856ab1f09a3ce9c0786227d88c0a Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Fri, 7 May 2021 12:13:05 +0200 Subject: [PATCH 10/10] Update stdlib/LinearAlgebra/test/uniformscaling.jl --- stdlib/LinearAlgebra/test/uniformscaling.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stdlib/LinearAlgebra/test/uniformscaling.jl b/stdlib/LinearAlgebra/test/uniformscaling.jl index aa2b3117ef456..b7b2e5c81cf88 100644 --- a/stdlib/LinearAlgebra/test/uniformscaling.jl +++ b/stdlib/LinearAlgebra/test/uniformscaling.jl @@ -462,7 +462,7 @@ end @test mul!(copy(C), J, a, alpha, beta) ≈ target_5mul target_5mul = beta*C # alpha = 0 @test mul!(copy(C), a, J, 0, beta) ≈ target_5mul - target_5mul = a*alpha*Matrix(J,n,n) # beta = 0 + target_5mul = a*alpha*Matrix(J, 3, 3) # beta = 0 @test mul!(copy(C), a, J, alpha, 0) ≈ target_5mul end