Merge pull request #441 from schmrlng/fix438

Cast expm Pade approximant coefficients to result eltype, fix #438

Author: SimonDanisch

src/expm.jl

@@ -68,48 +68,50 @@ end
 
 # Adapted from implementation in Base; algorithm from
 # Higham, "Functions of Matrices: Theory and Computation", SIAM, 2008
-function _exp(::Size, A::StaticMatrix{<:Any,<:Any,T}) where T
+function _exp(::Size, _A::StaticMatrix{<:Any,<:Any,T}) where T
+    S = typeof((zero(T)*zero(T) + zero(T)*zero(T))/one(T))
+    A = S.(_A)
     # omitted: matrix balancing, i.e., LAPACK.gebal!
     nA = maximum(sum(abs.(A), Val{1}))    # marginally more performant than norm(A, 1)
     ## For sufficiently small nA, use lower order Padé-Approximations
     if (nA <= 2.1)
         A2 = A*A
         if nA > 0.95
-            U = @evalpoly(A2, T(8821612800)*I, T(302702400)*I, T(2162160)*I, T(3960)*I, T(1)*I)
+            U = @evalpoly(A2, S(8821612800)*I, S(302702400)*I, S(2162160)*I, S(3960)*I, S(1)*I)
             U = A*U
-            V = @evalpoly(A2, T(17643225600)*I, T(2075673600)*I, T(30270240)*I, T(110880)*I, T(90)*I)
+            V = @evalpoly(A2, S(17643225600)*I, S(2075673600)*I, S(30270240)*I, S(110880)*I, S(90)*I)
         elseif nA > 0.25
-            U = @evalpoly(A2, T(8648640)*I, T(277200)*I, T(1512)*I, T(1)*I)
+            U = @evalpoly(A2, S(8648640)*I, S(277200)*I, S(1512)*I, S(1)*I)
             U = A*U
-            V = @evalpoly(A2, T(17297280)*I, T(1995840)*I, T(25200)*I, T(56)*I)
+            V = @evalpoly(A2, S(17297280)*I, S(1995840)*I, S(25200)*I, S(56)*I)
         elseif nA > 0.015
-            U = @evalpoly(A2, T(15120)*I, T(420)*I, T(1)*I)
+            U = @evalpoly(A2, S(15120)*I, S(420)*I, S(1)*I)
             U = A*U
-            V = @evalpoly(A2, T(30240)*I, T(3360)*I, T(30)*I)
+            V = @evalpoly(A2, S(30240)*I, S(3360)*I, S(30)*I)
         else
-            U = @evalpoly(A2, T(60)*I, T(1)*I)
+            U = @evalpoly(A2, S(60)*I, S(1)*I)
             U = A*U
-            V = @evalpoly(A2, T(120)*I, T(12)*I)
+            V = @evalpoly(A2, S(120)*I, S(12)*I)
         end
         expA = (V - U) \ (V + U)
     else
         s  = log2(nA/5.4)               # power of 2 later reversed by squaring
         if s > 0
             si = ceil(Int,s)
-            A = A / T(2^si)
+            A = A / S(2^si)
         end
 
         A2 = A*A
         A4 = A2*A2
         A6 = A2*A4
 
-        U = A6*(T(1)*A6 + T(16380)*A4 + T(40840800)*A2) +
-            (T(33522128640)*A6 + T(10559470521600)*A4 + T(1187353796428800)*A2) +
-            T(32382376266240000)*I
+        U = A6*(S(1)*A6 + S(16380)*A4 + S(40840800)*A2) +
+            (S(33522128640)*A6 + S(10559470521600)*A4 + S(1187353796428800)*A2) +
+            S(32382376266240000)*I
         U = A*U
-        V = A6*(T(182)*A6 + T(960960)*A4 + T(1323241920)*A2) +
-            (T(670442572800)*A6 + T(129060195264000)*A4 + T(7771770303897600)*A2) +
-            T(64764752532480000)*I
+        V = A6*(S(182)*A6 + S(960960)*A4 + S(1323241920)*A2) +
+            (S(670442572800)*A6 + S(129060195264000)*A4 + S(7771770303897600)*A2) +
+            S(64764752532480000)*I
         expA = (V - U) \ (V + U)
 
         if s > 0            # squaring to reverse dividing by power of 2