Hermitian, ishermitian now check diagonals

base/linalg/symmetric.jl

@@ -8,7 +8,15 @@ immutable Hermitian{T,S<:AbstractMatrix} <: AbstractMatrix{T}
     data::S
     uplo::Char
 end
-Hermitian(A::AbstractMatrix, uplo::Symbol=:U) = (chksquare(A);Hermitian{eltype(A),typeof(A)}(A, char_uplo(uplo)))
+function Hermitian(A::AbstractMatrix, uplo::Symbol=:U)
+    n = chksquare(A)
+    for i=1:n
+        isreal(A[i, i]) || throw(ArgumentError(
+            "Cannot construct Hermitian from matrix with nonreal diagonals"))
+    end
+    Hermitian{eltype(A),typeof(A)}(A, char_uplo(uplo))
+end
+
 typealias HermOrSym{T,S} Union(Hermitian{T,S}, Symmetric{T,S})
 typealias RealHermSymComplexHerm{T<:Real,S} Union(Hermitian{T,S}, Symmetric{T,S}, Hermitian{Complex{T},S})
 
@@ -51,27 +59,27 @@ inv{T<:BlasFloat,S<:StridedMatrix}(A::Hermitian{T,S}) = Hermitian{T,S}(inv(bkfac
 inv{T<:BlasFloat,S<:StridedMatrix}(A::Symmetric{T,S}) = Symmetric{T,S}(inv(bkfact(A.data, symbol(A.uplo), true)), A.uplo)
 
 eigfact!{T<:BlasReal,S<:StridedMatrix}(A::RealHermSymComplexHerm{T,S}) = Eigen(LAPACK.syevr!('V', 'A', A.uplo, A.data, 0.0, 0.0, 0, 0, -1.0)...)
-# Because of #6721 it is necessay to specify the parameters explicitly here.
+# Because of #6721 it is necessary to specify the parameters explicitly here.
 eigfact{T1<:Real,T2}(A::RealHermSymComplexHerm{T1,T2}) = (T = eltype(A); S = promote_type(Float32, typeof(zero(T)/norm(one(T)))); eigfact!(S != T ? convert(AbstractMatrix{S}, A) : copy(A)))
 
 eigfact!{T<:BlasReal,S<:StridedMatrix}(A::RealHermSymComplexHerm{T,S}, irange::UnitRange) = Eigen(LAPACK.syevr!('V', 'I', A.uplo, A.data, 0.0, 0.0, irange.start, irange.stop, -1.0)...)
-# Because of #6721 it is necessay to specify the parameters explicitly here.
+# Because of #6721 it is necessary to specify the parameters explicitly here.
 eigfact{T1<:Real,T2}(A::RealHermSymComplexHerm{T1,T2}, irange::UnitRange) = (T = eltype(A); S = promote_type(Float32, typeof(zero(T)/norm(one(T)))); eigfact!(S != T ? convert(AbstractMatrix{S}, A) : copy(A), irange))
 
 eigfact!{T<:BlasReal,S<:StridedMatrix}(A::RealHermSymComplexHerm{T,S}, vl::Real, vh::Real) = Eigen(LAPACK.syevr!('V', 'V', A.uplo, A.data, convert(T, vl), convert(T, vh), 0, 0, -1.0)...)
-# Because of #6721 it is necessay to specify the parameters explicitly here.
+# Because of #6721 it is necessary to specify the parameters explicitly here.
 eigfact{T1<:Real,T2}(A::RealHermSymComplexHerm{T1,T2}, vl::Real, vh::Real) = (T = eltype(A); S = promote_type(Float32, typeof(zero(T)/norm(one(T)))); eigfact!(S != T ? convert(AbstractMatrix{S}, A) : copy(A), vl, vh))
 
 eigvals!{T<:BlasReal,S<:StridedMatrix}(A::RealHermSymComplexHerm{T,S}) = LAPACK.syevr!('N', 'A', A.uplo, A.data, 0.0, 0.0, 0, 0, -1.0)[1]
-# Because of #6721 it is necessay to specify the parameters explicitly here.
+# Because of #6721 it is necessary to specify the parameters explicitly here.
 eigvals{T1<:Real,T2}(A::RealHermSymComplexHerm{T1,T2}) = (T = eltype(A); S = promote_type(Float32, typeof(zero(T)/norm(one(T)))); eigvals!(S != T ? convert(AbstractMatrix{S}, A) : copy(A)))
 
 eigvals!{T<:BlasReal,S<:StridedMatrix}(A::RealHermSymComplexHerm{T,S}, irange::UnitRange) = LAPACK.syevr!('N', 'I', A.uplo, A.data, 0.0, 0.0, irange.start, irange.stop, -1.0)[1]
-# Because of #6721 it is necessay to specify the parameters explicitly here.
+# Because of #6721 it is necessary to specify the parameters explicitly here.
 eigvals{T1<:Real,T2}(A::RealHermSymComplexHerm{T1,T2}, irange::UnitRange) = (T = eltype(A); S = promote_type(Float32, typeof(zero(T)/norm(one(T)))); eigvals!(S != T ? convert(AbstractMatrix{S}, A) : copy(A), irange))
 
 eigvals!{T<:BlasReal,S<:StridedMatrix}(A::RealHermSymComplexHerm{T,S}, vl::Real, vh::Real) = LAPACK.syevr!('N', 'V', A.uplo, A.data, convert(T, vl), convert(T, vh), 0, 0, -1.0)[1]
-# Because of #6721 it is necessay to specify the parameters explicitly here.
+# Because of #6721 it is necessary to specify the parameters explicitly here.
 eigvals{T1<:Real,T2}(A::RealHermSymComplexHerm{T1,T2}, vl::Real, vh::Real) = (T = eltype(A); S = promote_type(Float32, typeof(zero(T)/norm(one(T)))); eigvals!(S != T ? convert(AbstractMatrix{S}, A) : copy(A), vl, vh))
 
 eigmax{T<:Real,S<:StridedMatrix}(A::RealHermSymComplexHerm{T,S}) = eigvals(A, size(A, 1):size(A, 1))[1]

test/choosetests.jl

@@ -47,7 +47,7 @@ function choosetests(choices = [])
         prepend!(tests, ["linalg1", "linalg2", "linalg3", "linalg4",
             "linalg/lapack", "linalg/triangular", "linalg/tridiag",
             "linalg/pinv", "linalg/givens", "linalg/cholesky", "linalg/lu",
-            "linalg/arnoldi"])
+            "linalg/arnoldi", "linalg/symmetric"])
         end
 
     net_required_for = ["socket", "parallel"]

test/linalg/symmetric.jl

@@ -0,0 +1,83 @@
+using Base.Test
+
+debug = false #Turn on for more debugging info
+
+#Pauli σ-matrices
+for σ in map(Hermitian, Any[ eye(2), [0 1; 1 0], [0 -im; im 0], [1 0; 0 -1] ])
+    @test ishermitian(σ)
+end
+
+# Hermitian matrix exponential
+let A1 = randn(4,4) + im*randn(4,4)
+    A2 = A1 + A1'
+    @test_approx_eq expm(A2) expm(Hermitian(A2))
+end
+
+let n=10
+    areal = randn(n,n)/2
+    aimg  = randn(n,n)/2
+    debug && println("symmetric eigendecomposition")
+    for eltya in (Float32, Float64, Complex64, Complex128, BigFloat, Int)
+        a = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex(areal, aimg) : areal)
+        asym = a'+a                 # symmetric indefinite
+        ε = εa = eps(abs(float(one(eltya))))
+
+        debug && println("\ntype of a: ", eltya, "\n")
+
+        eltya == BigFloat && continue # Revisit when implemented in julia
+        d, v = eig(asym)
+        @test_approx_eq asym*v[:,1] d[1]*v[:,1]
+        @test_approx_eq v*Diagonal(d)*v' asym
+        @test isequal(eigvals(asym[1]), eigvals(asym[1:1,1:1]))
+        @test_approx_eq abs(eigfact(Hermitian(asym), 1:2)[:vectors]'v[:,1:2]) eye(eltya, 2)
+        eig(Hermitian(asym), 1:2) # same result, but checks that method works
+        @test_approx_eq abs(eigfact(Hermitian(asym), d[1]-10*eps(d[1]), d[2]+10*eps(d[2]))[:vectors]'v[:,1:2]) eye(eltya, 2)
+        eig(Hermitian(asym), d[1]-10*eps(d[1]), d[2]+10*eps(d[2])) # same result, but checks that method works
+        @test_approx_eq eigvals(Hermitian(asym), 1:2) d[1:2]
+        @test_approx_eq eigvals(Hermitian(asym), d[1]-10*eps(d[1]), d[2]+10*eps(d[2])) d[1:2]
+
+        # relation to svdvals
+        @test sum(sort(abs(eigvals(Hermitian(asym))))) == sum(sort(svdvals(Hermitian(asym))))
+
+        # cond
+        @test_approx_eq cond(Hermitian(asym)) cond(asym)
+
+        # rank
+        let A = a[:,1:5]*a[:,1:5]'
+            @test rank(A) == rank(Hermitian(A))
+        end
+    end
+end
+
+#Issue #7647: test xsyevr, xheevr, xstevr drivers
+for Mi7647 in (Symmetric(diagm(1.0:3.0)),
+               Hermitian(diagm(1.0:3.0)),
+               Hermitian(diagm(complex(1.0:3.0))),
+               SymTridiagonal([1.0:3.0;], zeros(2)))
+    debug && println("Eigenvalues in interval for $(typeof(Mi7647))")
+    @test eigmin(Mi7647)  == eigvals(Mi7647, 0.5, 1.5)[1] == 1.0
+    @test eigmax(Mi7647)  == eigvals(Mi7647, 2.5, 3.5)[1] == 3.0
+    @test eigvals(Mi7647) == eigvals(Mi7647, 0.5, 3.5) == [1.0:3.0;]
+end
+
+#Issue #7933
+let A7933 = [1 2; 3 4]
+    B7933 = copy(A7933)
+    C7933 = full(Symmetric(A7933))
+    @test A7933 == B7933
+end
+
+# Issues #8057 and #8058
+for f in (eigfact, eigvals, eig)
+    for A in (Symmetric([0 1; 1 0]), Hermitian([0 im; -im 0]))
+        @test_throws ArgumentError f(A, 3, 2)
+        @test_throws ArgumentError f(A, 1:4)
+    end
+end
+
+#Issue 10671
+let A = [1.0+im 2.0; 2.0 0.0]
+    @test !ishermitian(A)
+    @test_throws ArgumentError Hermitian(A)
+end
+

test/linalg1.jl

@@ -97,32 +97,6 @@ debug && println("(Automatic) Thin (pivoted) QR decomposition") # Pivoting is on
     @test_approx_eq q*r isa(qrpa, QRPivoted) ? a[:,p] : a[:,1:n1]
     @test_approx_eq isa(qrpa, QRPivoted) ? q*r[:,invperm(p)] : q*r a[:,1:n1]
 
-debug && println("symmetric eigen-decomposition")
-    if eltya != BigFloat && eltyb != BigFloat # Revisit when implemented in julia
-        d,v   = eig(asym)
-        @test_approx_eq asym*v[:,1] d[1]*v[:,1]
-        @test_approx_eq v*Diagonal(d)*v' asym
-        @test isequal(eigvals(asym[1]), eigvals(asym[1:1,1:1]))
-        @test_approx_eq abs(eigfact(Hermitian(asym), 1:2)[:vectors]'v[:,1:2]) eye(eltya, 2)
-        eig(Hermitian(asym), 1:2) # same result, but checks that method works
-        @test_approx_eq abs(eigfact(Hermitian(asym), d[1]-10*eps(d[1]), d[2]+10*eps(d[2]))[:vectors]'v[:,1:2]) eye(eltya, 2)
-        eig(Hermitian(asym), d[1]-10*eps(d[1]), d[2]+10*eps(d[2])) # same result, but checks that method works
-        @test_approx_eq eigvals(Hermitian(asym), 1:2) d[1:2]
-        @test_approx_eq eigvals(Hermitian(asym), d[1]-10*eps(d[1]), d[2]+10*eps(d[2])) d[1:2]
-
-        # relation to svdvals
-        @test sum(sort(abs(eigvals(Hermitian(asym))))) == sum(sort(svdvals(Hermitian(asym))))
-
-        # cond
-        @test_approx_eq cond(Hermitian(asym)) cond(asym)
-
-        # rank
-        let
-            A = a[:,1:5]*a[:,1:5]'
-            @test rank(A) == rank(Hermitian(A))
-        end
-    end
-
 debug && println("non-symmetric eigen decomposition")
     if eltya != BigFloat && eltyb != BigFloat # Revisit when implemented in julia
         d,v   = eig(a)

test/linalg3.jl

@@ -175,11 +175,6 @@ for elty in (Float32, Float64, Complex64, Complex128)
                                               0  -0.000000000000002   3.000000000000000])
 end
 
-# Hermitian matrix exponential
-A1 = randn(4,4) + im*randn(4,4)
-A2 = A1 + A1'
-@test_approx_eq expm(A2) expm(Hermitian(A2))
-
 # matmul for types w/o sizeof (issue #1282)
 A = Array(Complex{Int},10,10)
 A[:] = complex(1,1)

test/linalg4.jl

@@ -25,17 +25,6 @@ end
 
 n=12 #Size of matrix problem to test
 
-#Issue #7647: test xsyevr, xheevr, xstevr drivers
-for Mi7647 in (Symmetric(diagm(1.0:3.0)),
-               Hermitian(diagm(1.0:3.0)),
-               Hermitian(diagm(complex(1.0:3.0))),
-               SymTridiagonal([1.0:3.0;], zeros(2)))
-    debug && println("Eigenvalues in interval for $(typeof(Mi7647))")
-    @test eigmin(Mi7647)  == eigvals(Mi7647, 0.5, 1.5)[1] == 1.0
-    @test eigmax(Mi7647)  == eigvals(Mi7647, 2.5, 3.5)[1] == 3.0
-    @test eigvals(Mi7647) == eigvals(Mi7647, 0.5, 3.5) == [1.0:3.0;]
-end
-
 debug && println("Bidiagonal matrices")
 for relty in (Float32, Float64, BigFloat), elty in (relty, Complex{relty})
     debug && println("elty is $(elty), relty is $(relty)")
@@ -228,20 +217,6 @@ end
 @test_throws ErrorException transpose(sub(sprandn(10, 10, 0.3), 1:4, 1:4))
 @test_throws ErrorException ctranspose(sub(sprandn(10, 10, 0.3), 1:4, 1:4))
 
-# Issue #7933
-A7933 = [1 2; 3 4]
-B7933 = copy(A7933)
-C7933 = full(Symmetric(A7933))
-@test A7933 == B7933
-
-# Issues #8057 and #8058
-for f in (eigfact, eigvals, eig)
-    for A in (Symmetric(randn(2,2)), Hermitian(complex(randn(2,2), randn(2,2))))
-        @test_throws ArgumentError f(A, 3, 2)
-        @test_throws ArgumentError f(A, 1:4)
-    end
-end
-
 # test diag
 A = eye(4)
 @test diag(A) == ones(4)

test/sparsedir/cholmod.jl

@@ -180,7 +180,7 @@ ACSC = sprandn(10, 10, 0.3) + I
 @test ishermitian(Sparse(Hermitian(complex(ACSC), :U)))
 
 # test Sparse constructor for c_Sparse{Tv,Ti} input( and Sparse*Sparse)
-B = CHOLMOD.Sparse(SparseMatrixCSC{Float64,Int32}(sprandn(48, 48, 0.1))) # A has Int32 indeces
+B = CHOLMOD.Sparse(SparseMatrixCSC{Float64,Int32}(sprandn(48, 48, 0.1))) # A has Int32 indices
 @test_approx_eq A*B sparse(A)*sparse(B)
 
 # test Sparse constructor for c_SparseVoid (and read_sparse)
@@ -215,7 +215,7 @@ end
 @test_throws ArgumentError CHOLMOD.Sparse(convert(Ptr{CHOLMOD.C_SparseVoid}, C_NULL))
 
 ## The struct pointer must be constructed by the library constructor and then modified afterwards to checks that the method throws
-### illegal dtype (for now but should be supprted at some point)
+### illegal dtype (for now but should be supported at some point)
 p = ccall((:cholmod_l_allocate_sparse, :libcholmod), Ptr{CHOLMOD.C_SparseVoid},
     (Csize_t, Csize_t, Csize_t, Cint, Cint, Cint, Cint, Ptr{Void}),
     1, 1, 1, true, true, 0, CHOLMOD.REAL, CHOLMOD.common(CHOLMOD.SuiteSparse_long))
@@ -265,7 +265,7 @@ end
 
 # Test Dense wrappers (only Float64 supported a present)
 
-## High level interfact
+## High level interface
 for elty in (Float64, Complex{Float64})
     if elty == Float64
         A = randn(5, 5)
@@ -328,6 +328,7 @@ for elty in (Float64, Complex{Float64})
     @test_throws BoundsError A1Sparse[6, 1]
     @test_throws BoundsError A1Sparse[1, 6]
     @test sparse(A1Sparse) == A1
+    for i=1:size(A1, 1) A1[i, i] = real(A1[i, i]) end #Construct Hermitian matrix properly
     @test CHOLMOD.sparse(CHOLMOD.Sparse(Hermitian(A1, :L))) == Hermitian(A1, :L)
     @test CHOLMOD.sparse(CHOLMOD.Sparse(Hermitian(A1, :U))) == Hermitian(A1, :U)
     @test_throws ArgumentError convert(SparseMatrixCSC{elty,Int}, A1pdSparse)
@@ -338,7 +339,7 @@ for elty in (Float64, Complex{Float64})
     end
     @test copy(A1Sparse) == A1Sparse
     @test size(A1Sparse, 3) == 1
-    if elty <: Real # multiplcation only defined for real matrices in CHOLMOD
+    if elty <: Real # multiplication only defined for real matrices in CHOLMOD
         @test_approx_eq A1Sparse*A2Sparse A1*A2
         @test_throws DimensionMismatch CHOLMOD.Sparse(A1[:,1:4])*A2Sparse
         @test_approx_eq A1Sparse'A2Sparse A1'A2
@@ -359,8 +360,8 @@ for elty in (Float64, Complex{Float64})
     @test_throws ArgumentError cholfact(A1)
     @test_throws Base.LinAlg.PosDefException cholfact(A1 + A1' - 2eigmax(full(A1 + A1'))I)
     @test_throws Base.LinAlg.PosDefException cholfact(A1 + A1', -2eigmax(full(A1 + A1')))
-    @test_throws Base.LinAlg.ArgumentError ldltfact(A1 + A1' - 2real(A1[1,1])I)
-    @test_throws Base.LinAlg.ArgumentError ldltfact(A1 + A1', -2real(A1[1,1]))
+    @test_throws ArgumentError ldltfact(A1 + A1' - 2real(A1[1,1])I)
+    @test_throws ArgumentError ldltfact(A1 + A1', -2real(A1[1,1]))
     @test_throws ArgumentError cholfact(A1)
     @test_throws ArgumentError cholfact(A1, 1.0)
     @test_throws ArgumentError ldltfact(A1)