From 98791b1fde6e54d6ba2192b65e9593db362bf262 Mon Sep 17 00:00:00 2001 From: jishnub Date: Tue, 8 Feb 2022 15:21:32 +0400 Subject: [PATCH 1/3] Optimize real matrix * complex matrix by using real matmul --- stdlib/LinearAlgebra/src/matmul.jl | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/stdlib/LinearAlgebra/src/matmul.jl b/stdlib/LinearAlgebra/src/matmul.jl index 9338b5e10d73d..b7e13713d5222 100644 --- a/stdlib/LinearAlgebra/src/matmul.jl +++ b/stdlib/LinearAlgebra/src/matmul.jl @@ -174,12 +174,8 @@ function (*)(A::AdjOrTransStridedMat{<:BlasComplex}, B::StridedMaybeAdjOrTransMa wrapperop(B)(convert(AbstractArray{real(TS)}, _parent(B)))) end # the following case doesn't seem to benefit from the translation A*B = (B' * A')' -function (*)(A::StridedMatrix{<:BlasReal}, B::StridedMatrix{<:BlasComplex}) - TS = promote_type(eltype(A), eltype(B)) - mul!(similar(B, TS, (size(A, 1), size(B, 2))), - convert(AbstractArray{TS}, A), - convert(AbstractArray{TS}, B)) -end +(*)(A::StridedMatrix{<:BlasReal}, B::StridedMatrix{<:BlasComplex}) = + (A * real(B)) .+ im .* (A * imag(B)) (*)(A::AdjOrTransStridedMat{<:BlasReal}, B::StridedMatrix{<:BlasComplex}) = copy(transpose(transpose(B) * parent(A))) (*)(A::StridedMaybeAdjOrTransMat{<:BlasReal}, B::AdjOrTransStridedMat{<:BlasComplex}) = copy(wrapperop(B)(parent(B) * transpose(A))) From 3364366be3119368969563217b01b1cfed2af84e Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Tue, 8 Feb 2022 16:25:07 +0400 Subject: [PATCH 2/3] Broadcast Complex instead of + Co-authored-by: Daniel Karrasch --- stdlib/LinearAlgebra/src/matmul.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stdlib/LinearAlgebra/src/matmul.jl b/stdlib/LinearAlgebra/src/matmul.jl index b7e13713d5222..e3d74d6cf46a2 100644 --- a/stdlib/LinearAlgebra/src/matmul.jl +++ b/stdlib/LinearAlgebra/src/matmul.jl @@ -175,7 +175,7 @@ function (*)(A::AdjOrTransStridedMat{<:BlasComplex}, B::StridedMaybeAdjOrTransMa end # the following case doesn't seem to benefit from the translation A*B = (B' * A')' (*)(A::StridedMatrix{<:BlasReal}, B::StridedMatrix{<:BlasComplex}) = - (A * real(B)) .+ im .* (A * imag(B)) + Complex.(A * real(B), A * imag(B)) (*)(A::AdjOrTransStridedMat{<:BlasReal}, B::StridedMatrix{<:BlasComplex}) = copy(transpose(transpose(B) * parent(A))) (*)(A::StridedMaybeAdjOrTransMat{<:BlasReal}, B::AdjOrTransStridedMat{<:BlasComplex}) = copy(wrapperop(B)(parent(B) * transpose(A))) From 84375b2a10f8cb8dcc0fbe0de7656bfc1ed35262 Mon Sep 17 00:00:00 2001 From: jishnub Date: Tue, 8 Feb 2022 16:38:29 +0400 Subject: [PATCH 3/3] reduce allocations --- stdlib/LinearAlgebra/src/matmul.jl | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/stdlib/LinearAlgebra/src/matmul.jl b/stdlib/LinearAlgebra/src/matmul.jl index e3d74d6cf46a2..c1c60bfeb72bd 100644 --- a/stdlib/LinearAlgebra/src/matmul.jl +++ b/stdlib/LinearAlgebra/src/matmul.jl @@ -174,8 +174,13 @@ function (*)(A::AdjOrTransStridedMat{<:BlasComplex}, B::StridedMaybeAdjOrTransMa wrapperop(B)(convert(AbstractArray{real(TS)}, _parent(B)))) end # the following case doesn't seem to benefit from the translation A*B = (B' * A')' -(*)(A::StridedMatrix{<:BlasReal}, B::StridedMatrix{<:BlasComplex}) = - Complex.(A * real(B), A * imag(B)) +function (*)(A::StridedMatrix{<:BlasReal}, B::StridedMatrix{<:BlasComplex}) + temp = real(B) + R = A * temp + temp .= imag.(B) + I = A * temp + Complex.(R, I) +end (*)(A::AdjOrTransStridedMat{<:BlasReal}, B::StridedMatrix{<:BlasComplex}) = copy(transpose(transpose(B) * parent(A))) (*)(A::StridedMaybeAdjOrTransMat{<:BlasReal}, B::AdjOrTransStridedMat{<:BlasComplex}) = copy(wrapperop(B)(parent(B) * transpose(A)))