diff --git a/base/linalg/dense.jl b/base/linalg/dense.jl index fa36b882dc8ed..bde22d2c2d58a 100644 --- a/base/linalg/dense.jl +++ b/base/linalg/dense.jl @@ -397,6 +397,15 @@ end function (^)(A::AbstractMatrix{T}, p::Real) where T n = checksquare(A) + # Quicker return if A is diagonal + if isdiag(A) + retmat = copy(A) + for i in 1:n + retmat[i, i] = retmat[i, i] ^ p + end + return retmat + end + # For integer powers, use power_by_squaring isinteger(p) && return integerpow(A, p) @@ -408,15 +417,6 @@ function (^)(A::AbstractMatrix{T}, p::Real) where T return (Hermitian(A)^p) end - # Quicker return if A is diagonal - if isdiag(A) - retmat = copy(A) - for i in 1:n - retmat[i, i] = retmat[i, i] ^ p - end - return retmat - end - # Otherwise, use Schur decomposition return schurpow(A, p) end diff --git a/test/linalg/dense.jl b/test/linalg/dense.jl index 4b142dfc9891e..7eca457de6f2b 100644 --- a/test/linalg/dense.jl +++ b/test/linalg/dense.jl @@ -548,54 +548,38 @@ end #Aa : only positive real eigenvalues Aa = convert(Matrix{elty}, [5 4 2 1; 0 1 -1 -1; -1 -1 3 0; 1 1 -1 2]) - @test Aa^(1/2) ≈ sqrtm(Aa) - @test Aa^(-1/2) ≈ inv(sqrtm(Aa)) - @test Aa^(3/4) ≈ sqrtm(Aa) * sqrtm(sqrtm(Aa)) - @test Aa^(-3/4) ≈ inv(Aa) * sqrtm(sqrtm(Aa)) - @test Aa^(17/8) ≈ Aa^2 * sqrtm(sqrtm(sqrtm(Aa))) - @test Aa^(-17/8) ≈ inv(Aa^2 * sqrtm(sqrtm(sqrtm(Aa)))) - @test (Aa^0.2)^5 ≈ Aa - @test (Aa^(2/3))*(Aa^(1/3)) ≈ Aa - @test (Aa^im)^(-im) ≈ Aa - #Ab : both positive and negative real eigenvalues Ab = convert(Matrix{elty}, [1 2 3; 4 7 1; 2 1 4]) - @test Ab^(1/2) ≈ sqrtm(Ab) - @test Ab^(-1/2) ≈ inv(sqrtm(Ab)) - @test Ab^(3/4) ≈ sqrtm(Ab) * sqrtm(sqrtm(Ab)) - @test Ab^(-3/4) ≈ inv(Ab) * sqrtm(sqrtm(Ab)) - @test Ab^(17/8) ≈ Ab^2 * sqrtm(sqrtm(sqrtm(Ab))) - @test Ab^(-17/8) ≈ inv(Ab^2 * sqrtm(sqrtm(sqrtm(Ab)))) - @test (Ab^0.2)^5 ≈ Ab - @test (Ab^(2/3))*(Ab^(1/3)) ≈ Ab - @test (Ab^im)^(-im) ≈ Ab - #Ac : complex eigenvalues Ac = convert(Matrix{elty}, [5 4 2 1;0 1 -1 -1;-1 -1 3 6;1 1 -1 5]) - @test Ac^(1/2) ≈ sqrtm(Ac) - @test Ac^(-1/2) ≈ inv(sqrtm(Ac)) - @test Ac^(3/4) ≈ sqrtm(Ac) * sqrtm(sqrtm(Ac)) - @test Ac^(-3/4) ≈ inv(Ac) * sqrtm(sqrtm(Ac)) - @test Ac^(17/8) ≈ Ac^2 * sqrtm(sqrtm(sqrtm(Ac))) - @test Ac^(-17/8) ≈ inv(Ac^2 * sqrtm(sqrtm(sqrtm(Ac)))) - @test (Ac^0.2)^5 ≈ Ac - @test (Ac^(2/3))*(Ac^(1/3)) ≈ Ac - @test (Ac^im)^(-im) ≈ Ac - #Ad : defective Matrix Ad = convert(Matrix{elty}, [3 1; 0 3]) - @test Ad^(1/2) ≈ sqrtm(Ad) - @test Ad^(-1/2) ≈ inv(sqrtm(Ad)) - @test Ad^(3/4) ≈ sqrtm(Ad) * sqrtm(sqrtm(Ad)) - @test Ad^(-3/4) ≈ inv(Ad) * sqrtm(sqrtm(Ad)) - @test Ad^(17/8) ≈ Ad^2 * sqrtm(sqrtm(sqrtm(Ad))) - @test Ad^(-17/8) ≈ inv(Ad^2 * sqrtm(sqrtm(sqrtm(Ad)))) - @test (Ad^0.2)^5 ≈ Ad - @test (Ad^(2/3))*(Ad^(1/3)) ≈ Ad - @test (Ad^im)^(-im) ≈ Ad + #Ah : Hermitian Matrix + Ah = convert(Matrix{elty}, [3 1; 1 3]) + if elty <: Base.LinAlg.BlasComplex + Ah += [0 im; -im 0] + end + + #ADi : Diagonal Matrix + ADi = convert(Matrix{elty}, [3 0; 0 3]) + if elty <: Base.LinAlg.BlasComplex + ADi += [im 0; 0 im] + end + + for A in (Aa, Ab, Ac, Ad, Ah, ADi) + @test A^(1/2) ≈ sqrtm(A) + @test A^(-1/2) ≈ inv(sqrtm(A)) + @test A^(3/4) ≈ sqrtm(A) * sqrtm(sqrtm(A)) + @test A^(-3/4) ≈ inv(A) * sqrtm(sqrtm(A)) + @test A^(17/8) ≈ A^2 * sqrtm(sqrtm(sqrtm(A))) + @test A^(-17/8) ≈ inv(A^2 * sqrtm(sqrtm(sqrtm(A)))) + @test (A^0.2)^5 ≈ A + @test (A^(2/3))*(A^(1/3)) ≈ A + @test (A^im)^(-im) ≈ A + end end @testset "Least squares solutions" begin