Consolidate Hermitian and Symmetric routines

- Reimplement expm and sqrtm (ref #4006)
JuliaLang · Jan 26, 2014 · de6c1f6 · de6c1f6 · andreasnoack · Jan 26, 2014
1 parent 5b8e611
commit de6c1f6
Show file tree

Hide file tree

Showing 3 changed files with 63 additions and 117 deletions.
diff --git a/base/linalg.jl b/base/linalg.jl
@@ -197,7 +197,6 @@ include("linalg/factorization.jl")
 
 include("linalg/bunchkaufman.jl")
 include("linalg/triangular.jl")
-include("linalg/hermitian.jl")
 include("linalg/symmetric.jl")
 include("linalg/woodbury.jl")
 include("linalg/tridiag.jl")

diff --git a/base/linalg/hermitian.jl b/base/linalg/hermitian.jl
diff --git a/base/linalg/symmetric.jl b/base/linalg/symmetric.jl
@@ -1,66 +1,79 @@
-# Symmetric matrices
+## Hermitian matrices
 
-type Symmetric{T<:Number} <: AbstractMatrix{T}
-    S::Matrix{T}
-    uplo::Char
-end
-function Symmetric{T<:Number}(S::Matrix{T}, uplo::Symbol)
-    chksquare(S)
-    Symmetric(S, string(uplo)[1])
-end
-Symmetric(A::StridedMatrix) = Symmetric(A, :U)
+for ty in (:Hermitian, :Symmetric)
+    @eval begin
+        type $ty{T} <: AbstractMatrix{T}
+            S::Matrix{T}
+            uplo::Char
+        end
+        function $ty(S::Matrix, uplo::Symbol=:U)
+            chksquare(S)
+            $ty(S, string(uplo)[1])
+        end
 
-function copy!(A::Symmetric, B::Symmetric)
-    copy!(A.S, B.S)
-    A.uplo = B.uplo
-    A
+        function copy!(A::$ty, B::$ty)
+            copy!(A.S, B.S)
+            A.uplo = B.uplo
+            A
+        end
+        similar(A::$ty, args...) = $ty(similar(A.S, args...), A.uplo)
+    end
 end
-size(A::Symmetric, args...) = size(A.S, args...)
-getindex(A::Symmetric, i::Integer, j::Integer) = (A.uplo == 'U') == (i < j) ? getindex(A.S, i, j) : getindex(A.S, j, i)
+
+typealias HermOrSym Union(Hermitian, Symmetric)
+
+size(A::HermOrSym, args...) = size(A.S, args...)
+getindex(A::HermOrSym, i::Integer, j::Integer) = (A.uplo == 'U') == (i < j) ? getindex(A.S, i, j) : conj(getindex(A.S, j, i))
+full(A::Hermitian) = copytri!(A.S, A.uplo, true)
 full(A::Symmetric) = copytri!(A.S, A.uplo)
+ishermitian(A::Hermitian) = true
 ishermitian{T<:Real}(A::Symmetric{T}) = true
 ishermitian{T<:Complex}(A::Symmetric{T}) = all(imag(A.S) .== 0)
+issym{T<:Real}(A::Hermitian{T}) = true
+issym{T<:Complex}(A::Hermitian{T}) = all(imag(A.S) .== 0)
 issym(A::Symmetric) = true
 transpose(A::Symmetric) = A
-similar(A::Symmetric, args...) = Symmetric(similar(A.S, args...), A.uplo)
+ctranspose(A::Hermitian) = A
 
-*(A::Symmetric, B::Symmetric) = *(full(A), full(B))
-*(A::Symmetric, B::StridedMatrix) = *(full(A), B)
-*(A::StridedMatrix, B::Symmetric) = *(A, full(B))
+*(A::HermOrSym, B::HermOrSym) = full(A)*full(B)
+*(A::HermOrSym, B::StridedMatrix) = full(A)*B
+*(A::StridedMatrix, B::HermOrSym) = A*full(B)
 
-factorize!{T<:Real}(A::Symmetric{T}) = bkfact!(A.S, symbol(A.uplo))
-factorize!{T<:Complex}(A::Symmetric{T}) = bkfact!(A.S, symbol(A.uplo), true)
-\(A::Symmetric, B::StridedVecOrMat) = \(bkfact(A.S, symbol(A.uplo), true), B)
+factorize!(A::HermOrSym) = bkfact!(A.S, symbol(A.uplo), issym(A))
+\(A::HermOrSym, B::StridedVecOrMat) = \(bkfact(A.S, symbol(A.uplo), issym(A)), B)
 
-eigfact!{T<:BlasReal}(A::Symmetric{T}) = Eigen(LAPACK.syevr!('V', 'A', A.uplo, A.S, 0.0, 0.0, 0, 0, -1.0)...)
-eigfact(A::Symmetric) = eigfact!(copy(A))
-eigvals!{T<:BlasReal}(A::Symmetric{T}, il::Int, ih::Int) = LAPACK.syevr!('N', 'I', A.uplo, A.S, 0.0, 0.0, il, ih, -1.0)[1]
-eigvals!{T<:BlasReal}(A::Symmetric{T}, vl::Real, vh::Real) = LAPACK.syevr!('N', 'V', A.uplo, A.S, vl, vh, 0, 0, -1.0)[1]
-# eigvals!(A::Symmetric, args...) = eigvals!(float(A), args...)
-eigvals!(A::Symmetric) = eigvals!(A, 1, size(A, 1))
-eigmax(A::Symmetric) = eigvals(A, size(A, 1), size(A, 1))[1]
-eigmin(A::Symmetric) = eigvals(A, 1, 1)[1]
+# eigvals!(A::HermOrSym, args...) = eigvals!(float(A), args...)
+eigvals!(A::HermOrSym) = eigvals!(A, 1, size(A, 1))
+eigmax(A::HermOrSym) = eigvals(A, size(A, 1), size(A, 1))[1]
+eigmin(A::HermOrSym) = eigvals(A, 1, 1)[1]
+eigfact(A::HermOrSym) = eigfact!(copy(A))
 
-function eigfact!{T<:BlasReal}(A::Symmetric{T}, B::Symmetric{T})
-    vals, vecs, _ = LAPACK.sygvd!(1, 'V', A.uplo, A.S, B.uplo == A.uplo ? B.S : B.S')
-    GeneralizedEigen(vals, vecs)
+for (elty, ty) in ((:BlasFloat, :Hermitian), (:BlasReal, :Symmetric))
+    @eval begin
+        eigfact!{T<:$elty }(A::$ty{T}) = Eigen(LAPACK.syevr!('V', 'A', A.uplo, A.S, 0.0, 0.0, 0, 0, -1.0)...)
+        eigvals!{T<:$elty}(A::$ty{T}, il::Int, ih::Int) = LAPACK.syevr!('N', 'I', A.uplo, A.S, 0.0, 0.0, il, ih, -1.0)[1]
+        eigvals!{T<:$elty}(A::$ty{T}, vl::Real, vh::Real) = LAPACK.syevr!('N', 'V', A.uplo, A.S, vl, vh, 0, 0, -1.0)[1]
+        function eigfact!{T<:$elty}(A::$ty{T}, B::$ty{T})
+            vals, vecs, _ = LAPACK.sygvd!(1, 'V', A.uplo, A.S, B.uplo == A.uplo ? B.S : B.S')
+            GeneralizedEigen(vals, vecs)
+        end
+        eigfact(A::$ty, B::$ty) = eigfact!(copy(A), copy(B))
+        eigvals!{T<:$elty}(A::$ty{T}, B::$ty{T}) = LAPACK.sygvd!(1, 'N', A.uplo, A.S, B.uplo == A.uplo ? B.S : B.S')[1]
+    end
 end
-eigfact(A::Symmetric, B::Symmetric) = eigfact!(copy(A), copy(B))
-eigvals!{T<:BlasReal}(A::Symmetric{T}, B::Symmetric{T}) = LAPACK.sygvd!(1, 'N', A.uplo, A.S, B.uplo == A.uplo ? B.S : B.S')[1]
 
-function expm{T<:Real}(A::Symmetric{T})
-    F = eigfact(A)
-    scale(F[:vectors], exp(F[:values])) * F[:vectors]'
-end
+#Matrix-valued functions
+for (elty, ty) in ((:Any, :Hermitian), (:Real, :Symmetric))
+    @eval begin
+        function expm{T<:$elty}(A::$ty{T})
+            F = eigfact(A)
+            F.vectors * Diagonal(exp(F.values)) * F.vectors'
+        end
 
-function sqrtm{T<:Real}(A::Symmetric{T}, cond::Bool)
-    F = eigfact(A)
-    vsqrt = sqrt(complex(F[:values]))
-    if all(imag(vsqrt) .== 0)
-        retmat = copytri!(scale(F[:vectors], real(vsqrt)) * F[:vectors]', 'U')
-    else
-        zc = complex(F[:vectors])
-        retmat = copytri!(scale(zc, vsqrt) * zc', 'U')
+        function sqrtm{T<:$elty}(A::$ty{T}, cond::Bool=false)
+            F = eigfact(A)
+            retmat = F.vectors*Diagonal((isposdef(F) ? sqrt : x->sqrt(complex(x)))(F.values))*F.vectors'
+            return cond ? (retmat, norm(vsqrt, Inf)^2/norm(F[:values], Inf)) : retmat
+        end
     end
-    cond ? (retmat, norm(vsqrt, Inf)^2/norm(F[:values], Inf)) : retmat
-end
+end