Add macros to generate matrix functions for symmetric and Hermitian matrices, fix type returned by logm and sqrtm

Author: mfasi

base/linalg.jl

@@ -66,7 +66,6 @@ export
     eigvals!,
     eigvecs,
     expm,
-    sqrtm,
     eye,
     factorize,
     givens,
@@ -106,6 +105,7 @@ export
     schur,
     schurfact!,
     schurfact,
+    sqrtm,
     svd,
     svdfact!,
     svdfact,

base/linalg/dense.jl

@@ -189,7 +189,9 @@ expm(x::Number) = exp(x)
 ## "Functions of Matrices: Theory and Computation", SIAM
 function expm!{T<:BlasFloat}(A::StridedMatrix{T})
     n = chksquare(A)
-    n<2 && return exp(A)
+    if ishermitian(A)
+        return full(expm(Hermitian(A)))
+    end
     ilo, ihi, scale = LAPACK.gebal!('B', A)    # modifies A
     nA   = norm(A, 1)
     I    = eye(T,n)
@@ -278,10 +280,12 @@ function rcswap!{T<:Number}(i::Integer, j::Integer, X::StridedMatrix{T})
     end
 end
 
+expm(x::Number) = exp(x)
+
 function logm(A::StridedMatrix)
     # If possible, use diagonalization
     if ishermitian(A)
-        return logm(Hermitian(A))
+        return full(logm(Hermitian(A)))
     end
 
     # Use Schur decomposition
@@ -313,12 +317,13 @@ function logm(A::StridedMatrix)
         return retmat
     end
 end
+
 logm(a::Number) = (b = log(complex(a)); imag(b) == 0 ? real(b) : b)
 logm(a::Complex) = log(a)
 
 function sqrtm{T<:Real}(A::StridedMatrix{T})
     if issym(A)
-        return sqrtm(Symmetric(A))
+        return full(sqrtm(Symmetric(A)))
     end
     n = chksquare(A)
     SchurF = schurfact(complex(A))
@@ -328,13 +333,14 @@ function sqrtm{T<:Real}(A::StridedMatrix{T})
 end
 function sqrtm{T<:Complex}(A::StridedMatrix{T})
     if ishermitian(A)
-        return sqrtm(Hermitian(A))
+        return full(sqrtm(Hermitian(A)))
     end
     n = chksquare(A)
     SchurF = schurfact(A)
     R = full(sqrtm(UpperTriangular(SchurF[:T])))
     SchurF[:vectors]*R*SchurF[:vectors]'
 end
+
 sqrtm(a::Number) = (b = sqrt(complex(a)); imag(b) == 0 ? real(b) : b)
 sqrtm(a::Complex) = sqrt(a)
 

base/linalg/symmetric.jl

@@ -126,14 +126,28 @@ function svdvals!{T<:Real,S}(A::Union{Hermitian{T,S}, Symmetric{T,S}, Hermitian{
 end
 
 #Matrix-valued functions
-expm{T<:Real}(A::RealHermSymComplexHerm{T}) = (F = eigfact(A); F.vectors*Diagonal(exp(F.values))*F.vectors')
-function logm{T<:Real}(A::RealHermSymComplexHerm{T})
-    F = eigfact(A)
-    isposdef(F) && return F.vectors*Diagonal(log(F.values))*F.vectors'
-    return F.vectors*Diagonal(log(complex(F.values)))*F.vectors'
-end
-function sqrtm{T<:Real}(A::RealHermSymComplexHerm{T})
-    F = eigfact(A)
-    isposdef(F) && return F.vectors*Diagonal(sqrt(F.values))*F.vectors'
-    return F.vectors*Diagonal(sqrt(complex(F.values)))*F.vectors'
+for (funm, func) in ([:expm, :exp], [:logm,:log], [:sqrtm,:sqrt])
+
+    @eval begin
+        function ($funm)(A::Symmetric)
+            F = eigfact(A)
+            if isposdef(F)
+                retmat = F.vectors * Diagonal(($func)(F.values)) * F.vectors'
+            else
+                retmat = F.vectors * Diagonal(($func)(complex(F.values))) * F.vectors'
+            end
+            return Symmetric(retmat)
+        end
+
+        function ($funm)(A::Hermitian)
+            F = eigfact(A)
+            if isposdef(F)
+                retmat = Hermitian(F.vectors * Diagonal(($func)(F.values)) * F.vectors')
+            else
+                retmat = F.vectors * Diagonal(($func)(complex(F.values))) * F.vectors'
+            end
+            return retmat
+        end
+    end
+
 end