Make the algorithm for real powers of a matrix robust

Author: mfasi

base/linalg/diagonal.jl

@@ -251,6 +251,14 @@ expm(D::Diagonal) = Diagonal(exp.(D.diag))
 logm(D::Diagonal) = Diagonal(log.(D.diag))
 sqrtm(D::Diagonal) = Diagonal(sqrt.(D.diag))
 
+# Matrix functions
+^(D::Diagonal, p::Real) = Diagonal((D.diag).^p)
+for (funm, func) in ([:expm,:exp], [:sqrtm,:sqrt], [:logm,:log])
+    @eval begin
+        ($funm)(D::Diagonal) = Diagonal(($func)(D.diag))
+    end
+end
+
 #Linear solver
 function A_ldiv_B!(D::Diagonal, B::StridedVecOrMat)
     m, n = size(B, 1), size(B, 2)

base/linalg/triangular.jl

@@ -1672,7 +1672,7 @@ powm(A::LowerTriangular, p::Real) = powm(A.', p::Real).'
 # Based on the code available at http://eprints.ma.man.ac.uk/1851/02/logm.zip,
 # Copyright (c) 2011, Awad H. Al-Mohy and Nicholas J. Higham
 # Julia version relicensed with permission from original authors
-function logm{T<:Union{Float32,Float64,Complex{Float32},Complex{Float64}}}(A0::UpperTriangular{T})
+function logm{T<:Union{Float64,Complex{Float64}}}(A0::UpperTriangular{T})
     theta = [1.586970738772063e-005,
          2.313807884242979e-003,
          1.938179313533253e-002,
@@ -1718,7 +1718,7 @@ function logm{T<:Union{Float32,Float64,Complex{Float32},Complex{Float64}}}(A0::U
     scale!(2^s,Y.data)
 
     # Compute accurate diagonal and superdiagonal of log(A)
-    for k = 1:n-1
+    @inbounds for k = 1:n-1
         Ak = complex(A0[k,k])
         Akp1 = complex(A0[k+1,k+1])
         logAk = log(Ak)
@@ -1749,7 +1749,7 @@ logm(A::LowerTriangular) = logm(A.').'
 #      Numer. Algorithms, 59, (2012), 393–402.
 function sqrt_diag!(A0::UpperTriangular, A::UpperTriangular, s)
     n = checksquare(A0)
-    for i = 1:n
+    @inbounds for i = 1:n
         a = complex(A0[i,i])
         if s == 0
             A[i,i] = a - 1
@@ -1806,60 +1806,53 @@ function invsquaring(A0::UpperTriangular, theta)
     d2 = sqrt(norm(AmI^2, 1))
     d3 = cbrt(norm(AmI^3, 1))
     alpha2 = max(d2, d3)
-    foundm = false
     if alpha2 <= theta[2]
         m = alpha2<=theta[1]?1:2
-        foundm = true
-    end
-
-    while ~foundm
-        more = false
-        if s > s0
-            d3 = cbrt(norm(AmI^3, 1))
-        end
-        d4 = norm(AmI^4, 1)^(1/4)
-        alpha3 = max(d3, d4)
-        if alpha3 <= theta[tmax]
-            for j = 3:tmax
-                if alpha3 <= theta[j]
+    else
+        while true
+    	    more = false
+            if s > s0
+                d3 = cbrt(norm(AmI^3, 1))
+            end
+            d4 = norm(AmI^4, 1)^(1/4)
+            alpha3 = max(d3, d4)
+            if alpha3 <= theta[tmax]
+                for j = 3:tmax
+                    if alpha3 <= theta[j]
+                        break
+                    end
+                end
+                if j <= 6
+                    m = j
                     break
+                elseif alpha3 / 2 <= theta[5] && p < 2
+                    more = true
+                    p = p + 1
                 end
             end
-            if j <= 6
-                m = j
-                foundm = true
-                break
-            elseif alpha3 / 2 <= theta[5] && p < 2
-                more = true
-                p = p + 1
-           end
-        end
-
-        if ~more
-            d5 = norm(AmI^5, 1)^(1/5)
-            alpha4 = max(d4, d5)
-            eta = min(alpha3, alpha4)
-            if eta <= theta[tmax]
-                j = 0
-                for j = 6:tmax
-                    if eta <= theta[j]
-                        m = j
-                        break
+            if ~more
+                d5 = norm(AmI^5, 1)^(1/5)
+                alpha4 = max(d4, d5)
+                eta = min(alpha3, alpha4)
+                if eta <= theta[tmax]
+                    j = 0
+                    for j = 6:tmax
+                        if eta <= theta[j]
+                            m = j
+                            break
+                        end
                     end
+                    break
                 end
-                foundm = true
+            end
+            if s == maxsqrt
+                m = tmax
                 break
             end
+            A = sqrtm(A)
+            AmI = A - I
+            s = s + 1
         end
-
-        if s == maxsqrt
-            m = tmax
-            foundm = true
-            break
-        end
-        A = sqrtm(A)
-        AmI = A - I
-        s = s + 1
     end
 
     # Compute accurate superdiagonal of T