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,11 @@ 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)
+<<<<<<< 617b2a62223cd7cd0883a248a124b313bcf9d8b4
     for k = 1:n-1
+=======
+    @inbounds for k = 1:n-1
+>>>>>>> Make the algorithm for real powers of a matrix robust
         Ak = complex(A0[k,k])
         Akp1 = complex(A0[k+1,k+1])
         logAk = log(Ak)
@@ -1749,7 +1753,11 @@ logm(A::LowerTriangular) = logm(A.').'
 #      Numer. Algorithms, 59, (2012), 393–402.
 function sqrt_diag!(A0::UpperTriangular, A::UpperTriangular, s)
     n = checksquare(A0)
+<<<<<<< 617b2a62223cd7cd0883a248a124b313bcf9d8b4
     for i = 1:n
+=======
+    @inbounds for i = 1:n
+>>>>>>> Make the algorithm for real powers of a matrix robust
         a = complex(A0[i,i])
         if s == 0
             A[i,i] = a - 1
@@ -1806,25 +1814,31 @@ 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
+<<<<<<< 617b2a62223cd7cd0883a248a124b313bcf9d8b4
             if j <= 6
                 m = j
                 foundm = true
@@ -1845,21 +1859,45 @@ function invsquaring(A0::UpperTriangular, theta)
                     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
+>>>>>>> Make the algorithm for real powers of a matrix robust
                     end
+                    break
                 end
+<<<<<<< 617b2a62223cd7cd0883a248a124b313bcf9d8b4
                 foundm = true
                 break
+=======
+>>>>>>> Make the algorithm for real powers of a matrix robust
             end
-        end
 
+<<<<<<< 617b2a62223cd7cd0883a248a124b313bcf9d8b4
         if s == maxsqrt
             m = tmax
             foundm = true
             break
+=======
+            if s == maxsqrt
+                m = tmax
+                break
+            end
+            A = sqrtm(A)
+            AmI = A - I
+            s = s + 1
+>>>>>>> Make the algorithm for real powers of a matrix robust
         end
-        A = sqrtm(A)
-        AmI = A - I
-        s = s + 1
     end
 
     # Compute accurate superdiagonal of T