Move up the diagonal check in matrix power (#23158)

* move up diagonal check in matrix power * Tests for A^p for (secretly) Hermitian and diagonal
JuliaLang · Aug 7, 2017 · ed746ed · ed746ed
1 parent 9363b49
commit ed746ed
Show file tree

Hide file tree

Showing 2 changed files with 32 additions and 48 deletions.
diff --git 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
@@ -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