Make warnings optional and test non-hermitian case.

src/spectralanalysis.jl

@@ -5,6 +5,8 @@ using ..bases, ..states, ..operators, ..operators_dense, ..operators_sparse
 export eigenstates, eigenenergies, simdiag
 
 
+warn_nonhermitian() = warn("The given operator is not hermitian. If this is due to a numerical error make the operator hermitian first by calculating (x+dagger(x))/2 first.")
+
 """
     eigenstates(op::Operator[, n::Int])
 
@@ -16,13 +18,14 @@ about the way the calculation is done is needed, use the functions directly.
 More details can be found at
 [http://docs.julialang.org/en/stable/stdlib/linalg/].
 """
-function eigenstates(op::DenseOperator, n::Int=length(basis(op)))
+function eigenstates(op::DenseOperator, n::Int=length(basis(op)); warning=true)
     b = basis(op)
     if ishermitian(op)
         D, V = eig(Hermitian(op.data), 1:n)
         states = [Ket(b, V[:, k]) for k=1:length(D)]
         return D, states
     else
+        warning && warn_nonhermitian()
         D, V = eig(op.data)
         states = [Ket(b, V[:, k]) for k=1:length(D)]
         perm = sortperm(D, by=real)
@@ -35,11 +38,12 @@ end
 """
 For sparse operators by default it only returns the 6 lowest eigenvalues.
 """
-function eigenstates(op::SparseOperator, n::Int=length(basis(op)))
+function eigenstates(op::SparseOperator, n::Int=length(basis(op)); warning=true)
     b = basis(op)
     if ishermitian(op)
         data = Hermitian(op.data)
     else
+        warning && warn_nonhermitian()
         data = op.data
     end
     D, V = eigs(data; nev=n, which=:SR)
@@ -58,12 +62,13 @@ about the way the calculation is done is needed, use the function directly.
 More details can be found at
 [http://docs.julialang.org/en/stable/stdlib/linalg/].
 """
-function eigenenergies(op::DenseOperator, n::Int=length(basis(op)))
+function eigenenergies(op::DenseOperator, n::Int=length(basis(op)); warning=true)
     b = basis(op)
     if ishermitian(op)
         D = eigvals(Hermitian(op.data), 1:n)
         return D
     else
+        warning && warn_nonhermitian()
         D = eigvals(op.data)
         sort!(D, by=real)
         return D[1:n]
@@ -73,7 +78,7 @@ end
 """
 For sparse operators by default it only returns the 6 lowest eigenvalues.
 """
-eigenenergies(op::SparseOperator, n::Int=6) = eigenstates(op, n)[1]
+eigenenergies(op::SparseOperator, n::Int=6; warning=true) = eigenstates(op, n; warning=warning)[1]
 
 
 arithmetic_unary_error = operators.arithmetic_unary_error

test/test_spectralanalysis.jl

@@ -16,6 +16,8 @@ sprandop(b) = sparse(DenseOperator(b, rand(Complex128, length(b), length(b))))
 @test_throws ArgumentError eigenenergies(SpectralanalysisTestOperator())
 @test_throws bases.IncompatibleBases eigenenergies(DenseOperator(GenericBasis(3), GenericBasis(4)))
 
+
+# Test hermitian diagonalization
 b = GenericBasis(5)
 a = randoperator(b)
 H = (a+dagger(a))/2
@@ -30,8 +32,6 @@ R = U*D*dagger(U)
 Rsp = sparse(R)
 @test eigenenergies((R+dagger(R))/2) ≈ d
 @test eigenenergies((Rsp+dagger(Rsp))/2, 3) ≈ d[1:3]
-@test eigenenergies(R) ≈ d
-@test eigenenergies(Rsp, 3) ≈ d[1:3]
 
 E, states = eigenstates((R+dagger(R))/2, 3)
 Esp, states_sp = eigenstates((Rsp+dagger(Rsp))/2, 3)
@@ -44,6 +44,35 @@ for i=1:3
     @test states_sp[i].data*v ≈ U.data[:,i]
 end
 
+
+# Test nonhermitian diagonalization
+b = GenericBasis(5)
+a = randoperator(b)
+H = (a+dagger(a))/2
+U = expm(1im*H)
+d = [-3+0.2im, -2.6-0.1im, -0.1+0.5im, 0.0, 0.6+0.3im]
+D = DenseOperator(b, diagm(d))
+Dsp = sparse(D)
+@test eigenenergies(D; warning=false) ≈ d
+@test eigenenergies(Dsp, 3; warning=false) ≈ d[1:3]
+
+R = U*D*dagger(U)
+Rsp = sparse(R)
+@test eigenenergies(R; warning=false) ≈ d
+@test eigenenergies(Rsp, 3; warning=false) ≈ d[1:3]
+
+E, states = eigenstates(R, 3; warning=false)
+Esp, states_sp = eigenstates(Rsp, 3; warning=false)
+for i=1:3
+    @test E[i] ≈ d[i]
+    @test Esp[i] ≈ d[i]
+    v = U.data[1,i]/states[i].data[1]
+    @test states[i].data*v ≈ U.data[:,i]
+    v = U.data[1,i]/states_sp[i].data[1]
+    @test states_sp[i].data*v ≈ U.data[:,i]
+end
+
+
 # Test simdiag
 spinbasis = SpinBasis(1//2)
 sx = sigmax(spinbasis)