Inline documentation and tests for real power of a matrix

Author: mfasi

base/linalg/triangular.jl

@@ -997,6 +997,12 @@ function ^(A::UpperTriangular, p::Real)
         return powm(A,p - ceil(p)) * A^ceil(p)
     end
 end
+# Complex matrix power for upper triangular factor, see:
+#   Higham and Lin, "A Schur-Padé algorithm for fractional powers of a Matrix",
+#     SIAM J. Matrix Anal. & Appl., 32 (3), (2011) 1056–1078.
+#   Higham and Lin, "An improved Schur-Padé algorithm for fractional powers of
+#     a matrix and their Fréchet derivatives", SIAM. J. Matrix Anal. & Appl.,
+#     34(3), (2013) 1341–1360.
 function powm{T<:Union{Float64,Complex{Float64}}}(A0::UpperTriangular{T}, p::Real)
 
     if abs(p) >= 1
@@ -1153,6 +1159,8 @@ end
 logm(A::LowerTriangular) = logm(A.').'
 
 # Auxiliary functions for logm and matrix power
+
+# Compute accurate diagonal of A = A0^s - I
 function sqrt_diag!(A0::UpperTriangular, A::UpperTriangular, s)
     n = chksquare(A0)
     for i = 1:n
@@ -1177,6 +1185,7 @@ function sqrt_diag!(A0::UpperTriangular, A::UpperTriangular, s)
     end
 end
 
+# Compute the repeated square root of A to match the conditions on theta
 function invsquaring{T}(A0::UpperTriangular{T}, theta)
     maxsqrt = 100
     tmax = size(theta, 1)
@@ -1274,6 +1283,7 @@ function invsquaring{T}(A0::UpperTriangular{T}, theta)
 
 end
 
+# Compute accurate diagonal and superdiagonal of A = A0^p
 function blockpower!(A0::UpperTriangular, A::UpperTriangular, p)
     n = chksquare(A0)
     for k = 1:n-1
@@ -1300,6 +1310,8 @@ function blockpower!(A0::UpperTriangular, A::UpperTriangular, p)
         end
     end
 end
+
+# Unwinding number
 unw(x::Real) = 0
 unw(x::Number) = ceil((imag(x) - pi) / (2 * pi))
 

test/linalg/dense.jl

@@ -387,6 +387,40 @@ eA10 = [ 1.0+0.0im   0.0+0.0im                 0.0+0.0im                0.0+0.0i
         0.0+0.0im   0.0+0.0im                 0.0+0.0im                1.0+0.0im]
 @test_approx_eq expm(A10) eA10
 
+for elty in (Float64, Complex{Float64})
+    A11  = convert(Matrix{elty}, [1 2 3; 4 7 1; 2 1 4])
+
+    OLD_STDERR = STDERR
+    rd,wr = redirect_stderr()
+
+    @test A11^(1/2) ≈ sqrtm(A11)
+    s = readline(rd)
+    @test contains(s, "WARNING: Matrix with nonpositive real eigenvalues, a nonprincipal matrix power will be returned.")
+
+    @test A11^(-1/2) ≈ inv(sqrtm(A11))
+    s = readline(rd)
+    @test contains(s, "WARNING: Matrix with nonpositive real eigenvalues, a nonprincipal matrix power will be returned.")
+
+    @test A11^(3/4) ≈ A11 * inv(sqrtm(sqrtm(A11)))
+    s = readline(rd)
+    @test contains(s, "WARNING: Matrix with nonpositive real eigenvalues, a nonprincipal matrix power will be returned.")
+
+    @test A11^(-3/4) ≈ inv(A11) * sqrtm(sqrtm(A11))
+    s = readline(rd)
+    @test contains(s, "WARNING: Matrix with nonpositive real eigenvalues, a nonprincipal matrix power will be returned.")
+
+    @test A11^(17/8) ≈ A11^2 * sqrtm(sqrtm(sqrtm(A11)))
+    s = readline(rd)
+    @test contains(s, "WARNING: Matrix with nonpositive real eigenvalues, a nonprincipal matrix power will be returned.")
+
+    @test A11^(-17/8) ≈ inv(A11^2 * sqrtm(sqrtm(sqrtm(A11))))
+    s = readline(rd)
+    @test contains(s, "WARNING: Matrix with nonpositive real eigenvalues, a nonprincipal matrix power will be returned.")
+
+    redirect_stderr(OLD_STDERR)
+end
+
+
 # issue #7181
 A = [ 1  5  9
       2  6 10