Skip to content
Permalink
Browse files

Fix diagonalization and conjugate issues

  • Loading branch information...
david-pl committed May 29, 2018
1 parent 44f6326 commit 5d3c7ae62c5bfd4129de7a5fe981b6659b44c220
@@ -4,7 +4,7 @@ export Operator, length, basis, dagger, ishermitian, tensor, embed,
trace, ptrace, normalize, normalize!, expect, variance,
expm, permutesystems, identityoperator

import Base: ==, +, -, *, /, length, trace, one, ishermitian, expm
import Base: ==, +, -, *, /, length, trace, one, ishermitian, expm, conj, conj!
import ..bases: basis, tensor, ptrace, permutesystems,
samebases, check_samebases, multiplicable
import ..states: dagger, normalize, normalize!
@@ -50,6 +50,9 @@ basis(a::Operator) = (check_samebases(a); a.basis_l)

dagger(a::Operator) = arithmetic_unary_error("Hermitian conjugate", a)

conj(a::Operator) = arithmetic_unary_error("Complex conjugate", a)
conj!(a::Operator) = conj(a::Operator)

"""
ishermitian(op::Operator)
@@ -85,6 +85,10 @@ operators.dagger(x::DenseOperator) = DenseOperator(x.basis_r, x.basis_l, x.data'
operators.ishermitian(A::DenseOperator) = ishermitian(A.data)

operators.tensor(a::DenseOperator, b::DenseOperator) = DenseOperator(tensor(a.basis_l, b.basis_l), tensor(a.basis_r, b.basis_r), kron(b.data, a.data))

operators.conj(a::DenseOperator) = DenseOperator(a.basis_l, a.basis_r, conj(a.data))
operators.conj!(a::DenseOperator) = conj!(a.data)

"""
tensor(x::Ket, y::Bra)
@@ -75,6 +75,9 @@ operators.tensor(a::SparseOperator, b::DenseOperator) = SparseOperator(tensor(a.

operators.trace(op::SparseOperator) = (check_samebases(op); trace(op.data))

operators.conj(op::SparseOperator) = SparseOperator(op.basis_l, op.basis_r, conj(op.data))
operators.conj!(op::SparseOperator) = conj!(op.data)

function operators.ptrace(op::SparseOperator, indices::Vector{Int})
operators.check_ptrace_arguments(op, indices)
shape = [op.basis_l.shape; op.basis_r.shape]
@@ -18,6 +18,10 @@ 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/].
NOTE: Especially for small systems full diagonalization with Julia's `eig`
function is often more desirable. You can convert a sparse operator `A` to a
dense one using `full(A)`.
If the given operator is non-hermitian a warning is given. This behavior
can be turned off using the keyword `warning=false`.
"""
@@ -41,15 +45,13 @@ end
"""
For sparse operators by default it only returns the 6 lowest eigenvalues.
"""
function eigenstates(op::SparseOperator, n::Int=length(basis(op)); warning=true)
function eigenstates(op::SparseOperator; warning=true, kwargs...)
b = basis(op)
if ishermitian(op)
data = Hermitian(op.data)
else
warning && warn(nonhermitian_warning)
data = op.data
# TODO: Change to sparese-Hermitian specific algorithm if more efficient
if !ishermitian(op) && warning
warn(nonhermitian_warning)
end
D, V = eigs(data; nev=n, which=:SR)
D, V = eigs(op.data; which=:SR, kwargs...)
states = [Ket(b, V[:, k]) for k=1:length(D)]
D, states
end
@@ -72,6 +72,10 @@ xbra2 = Bra(b_l, rand(Complex128, length(b_l)))
@test 1e-12 > D((op1 + op2)*dagger(0.3*op3), 0.3*op1*dagger(op3) + 0.3*op2*dagger(op3))
@test 1e-12 > D(0.3*(op1*dagger(op2)), op1*(0.3*dagger(op2)))

tmp = copy(op1)
conj!(tmp)
@test tmp == conj(op1) && conj(tmp.data) == op1.data

# Internal layout
b1 = GenericBasis(2)
b2 = GenericBasis(3)
@@ -87,6 +87,11 @@ xbra2 = dagger(randstate(b_l))
# Test division
@test 1e-14 > D(op1/7, op1_/7)

# Conjugation
tmp = copy(op1)
conj!(tmp)
@test tmp == conj(op1) && conj(tmp.data) == op1.data

# Test identityoperator
Idense = identityoperator(DenseOperator, b_r)
I = identityoperator(SparseOperator, b_r)

0 comments on commit 5d3c7ae

Please sign in to comment.
You can’t perform that action at this time.