From b1720283917939ca5675db0517f2f31bab05e7ae Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Sun, 22 Dec 2024 19:55:39 +0530 Subject: [PATCH 1/3] `sqrt`, `cbrt` and `log` for dense diagonal matrices --- src/dense.jl | 12 +++++++++--- test/dense.jl | 9 +++++++++ 2 files changed, 18 insertions(+), 3 deletions(-) diff --git a/src/dense.jl b/src/dense.jl index 5e479841..f7fc3461 100644 --- a/src/dense.jl +++ b/src/dense.jl @@ -890,7 +890,9 @@ julia> log(A) """ function log(A::AbstractMatrix) # If possible, use diagonalization - if ishermitian(A) + if isdiag(A) + return applydiagonal(log, A) + elseif ishermitian(A) logHermA = log(Hermitian(A)) return ishermitian(logHermA) ? copytri!(parent(logHermA), 'U', true) : parent(logHermA) elseif istriu(A) @@ -969,7 +971,9 @@ sqrt(::AbstractMatrix) function sqrt(A::AbstractMatrix{T}) where {T<:Union{Real,Complex}} if checksquare(A) == 0 - return copy(A) + return copy(float(A)) + elseif isdiag(A) + return applydiagonal(sqrt, A) elseif ishermitian(A) sqrtHermA = sqrt(Hermitian(A)) return ishermitian(sqrtHermA) ? copytri!(parent(sqrtHermA), 'U', true) : parent(sqrtHermA) @@ -1035,7 +1039,9 @@ true """ function cbrt(A::AbstractMatrix{<:Real}) if checksquare(A) == 0 - return copy(A) + return copy(float(A)) + elseif isdiag(A) + return applydiagonal(cbrt, A) elseif issymmetric(A) return cbrt(Symmetric(A, :U)) else diff --git a/test/dense.jl b/test/dense.jl index a7616e2f..82ea68c3 100644 --- a/test/dense.jl +++ b/test/dense.jl @@ -814,6 +814,7 @@ end A13 = convert(Matrix{elty}, [2 0; 0 2]) @test typeof(log(A13)) == Array{elty, 2} + @test exp(log(A13)) ≈ log(exp(A13)) ≈ A13 T = elty == Float64 ? Symmetric : Hermitian @test typeof(log(T(A13))) == T{elty, Array{elty, 2}} @@ -965,6 +966,10 @@ end @test typeof(sqrt(A8)) == Matrix{elty} end end +@testset "sqrt for diagonal" begin + A = diagm(0 => [1, 2, 3]) + @test sqrt(A)^2 ≈ A +end @testset "issue #40141" begin x = [-1 -eps() 0 0; eps() -1 0 0; 0 0 -1 -eps(); 0 0 eps() -1] @@ -1297,6 +1302,10 @@ end T = cbrt(A) @test T*T*T ≈ A @test eltype(A) == eltype(T) + @testset "diagonal" begin + A = diagm(0 => [1, 2, 3]) + @test cbrt(A)^3 ≈ A + end end @testset "tr" begin From 454302156e0a45dd410d9a930484db3142d7f36c Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Sun, 22 Dec 2024 23:17:59 +0530 Subject: [PATCH 2/3] Tests for `cbrt` of empty matrix --- test/dense.jl | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/test/dense.jl b/test/dense.jl index 82ea68c3..34f7f5ee 100644 --- a/test/dense.jl +++ b/test/dense.jl @@ -1306,6 +1306,12 @@ end A = diagm(0 => [1, 2, 3]) @test cbrt(A)^3 ≈ A end + @testset "empty" begin + A = Matrix{Float64}(undef, 0, 0) + @test cbrt(A) == A + A = Matrix{Int}(undef, 0, 0) + @test cbrt(A) isa Matrix{Float64} + end end @testset "tr" begin From 4f02f26b45694f7d9688bf46d680c18c648ce39d Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Mon, 23 Dec 2024 11:00:53 +0530 Subject: [PATCH 3/3] Test for a dense symmetric matrix --- test/dense.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/dense.jl b/test/dense.jl index 34f7f5ee..e773a0bd 100644 --- a/test/dense.jl +++ b/test/dense.jl @@ -1282,6 +1282,7 @@ end T = cbrt(Symmetric(S,:U)) @test T*T*T ≈ S @test eltype(S) == eltype(T) + @test cbrt(Array(Symmetric(S,:U))) == T # Real valued symmetric S = (A -> (A+A')/2)(randn(N,N)) T = cbrt(Symmetric(S,:L))