Skip to content

Commit

Permalink
Browse files Browse the repository at this point in the history
New lazy tensor implementation.
Use two arrays "indices" and "operators" instead of a dictionary
i=>op_i.
  • Loading branch information
bastikr committed Mar 6, 2017
1 parent 22aeb8c commit 4bc626d
Show file tree
Hide file tree
Showing 12 changed files with 225 additions and 192 deletions.
34 changes: 21 additions & 13 deletions src/correlationexpansion.jl
Expand Up @@ -6,6 +6,7 @@ import ..operators: dagger, identityoperator,
gemv!, gemm!

using Combinatorics, Iterators
using ..sortedindices
using ..bases
using ..operators
using ..operators_dense
Expand Down Expand Up @@ -33,7 +34,8 @@ end

mask2indices(mask::Mask) = find(mask)

complement(N::Int, indices::Vector{Int}) = Int[i for i=1:N if i indices]
complement = sortedindices.complement
# complement(N::Int, indices::Vector{Int}) = Int[i for i=1:N if i ∉ indices]

as_mask(N::Int, m::Mask) = (@assert length(m)==N; m)
function as_mask(N::Int, m)
Expand Down Expand Up @@ -281,23 +283,26 @@ function gemm!(alpha, a::LazyTensor, b::CorrelationExpansion, beta, result::Corr
result.factor = a.factor*alpha*b.factor
@assert beta == complex(0.)
alpha = complex(1.)
n = 1
Na = length(a.indices)
for i=1:N
if i keys(a.operators)
gemm!(alpha, a.operators[i], b.operators[i], beta, result.operators[i])
if n <= Na && i == a.indices[n]
gemm!(alpha, a.operators[n], b.operators[i], beta, result.operators[i])
n += 1
else
copy!(result.operators[i].data, b.operators[i].data)
end
end
for mask in keys(b.correlations)
I = mask2indices(mask)
I_ = I keys(a.operators)
I_ = sortedindices.intersect(I, a.indices)
op = result.correlations[mask]
if isempty(I_)
copy!(op.data, b.correlations[mask].data)
else
indices = reduced_indices(I, I_)
operators = [a.operators[i] for i in I_]
a_ = embed(op.basis_l, op.basis_r, indices, operators)
operators = operators_lazy.suboperators(a, I_)
sortedindices.reducedindices!(I_, I)
a_ = LazyTensor(op.basis_l, op.basis_r, I_, operators)
gemm!(alpha, a_, b.correlations[mask], beta, op)
end
end
Expand All @@ -308,23 +313,26 @@ function gemm!(alpha, a::CorrelationExpansion, b::LazyTensor, beta, result::Corr
result.factor = a.factor*alpha*b.factor
@assert beta == complex(0.)
alpha = complex(1.)
n = 1
Nb = length(b.indices)
for i=1:N
if i keys(b.operators)
gemm!(alpha, a.operators[i], b.operators[i], beta, result.operators[i])
if n <= Nb && i == b.indices[n]
gemm!(alpha, a.operators[i], b.operators[n], beta, result.operators[i])
n += 1
else
copy!(result.operators[i].data, a.operators[i].data)
end
end
for mask in keys(a.correlations)
I = mask2indices(mask)
I_ = I keys(b.operators)
I_ = sortedindices.intersect(I, b.indices)
op = result.correlations[mask]
if isempty(I_)
copy!(op.data, a.correlations[mask].data)
else
indices = reduced_indices(I, I_)
operators = [b.operators[i] for i in I_]
b_ = embed(op.basis_l, op.basis_r, indices, operators)
operators = operators_lazy.suboperators(b, I_)
sortedindices.reducedindices!(I_, I)
b_ = LazyTensor(op.basis_l, op.basis_r, I_, operators)
gemm!(alpha, a.correlations[mask], b_, beta, op)
end
end
Expand Down
88 changes: 47 additions & 41 deletions src/cumulantexpansion.jl
Expand Up @@ -26,36 +26,36 @@ ProductDensityOperator() = error("ProductDensityOperator needs at least one oper
ProductDensityOperator(states::Ket...) = ProductDensityOperator(DenseOperator[tensor(state, dagger(state)) for state in states]...)

*(a::ProductDensityOperator, b::ProductDensityOperator) = (check_multiplicable(a.basis_r, b.basis_l); ProductDensityOperator(a.basis_l, b.basis_r, DenseOperator[a_i*b_i for (a_i, b_i) in zip(a.operators, b.operators)]))
function *(a::LazyTensor, b::ProductDensityOperator)
check_multiplicable(a.basis_r, b.basis_l);
operators = DenseOperator[]
for (alpha, b_alpha) in enumerate(b.operators)
if alpha in keys(a.operators)
push!(operators, a.operators[alpha]*b_alpha)
else
push!(operators, deepcopy(b_alpha))
end
end
if a.factor != 1.
operators[0] *= a.factor
end
ProductDensityOperator(operators...)
end
function *(a::ProductDensityOperator, b::LazyTensor)
check_multiplicable(a.basis_r, b.basis_l);
operators = DenseOperator[]
for (alpha, a_alpha) in enumerate(a.operators)
if alpha in keys(b.operators)
push!(operators, a_alpha*b.operators[alpha])
else
push!(operators, deepcopy(a_alpha))
end
end
if b.factor != 1.
operators[0] *= b.factor
end
ProductDensityOperator(operators...)
end
# function *(a::LazyTensor, b::ProductDensityOperator)
# check_multiplicable(a.basis_r, b.basis_l);
# operators = DenseOperator[]
# for (alpha, b_alpha) in enumerate(b.operators)
# if alpha in keys(a.operators)
# push!(operators, a.operators[alpha]*b_alpha)
# else
# push!(operators, deepcopy(b_alpha))
# end
# end
# if a.factor != 1.
# operators[0] *= a.factor
# end
# ProductDensityOperator(operators...)
# end
# function *(a::ProductDensityOperator, b::LazyTensor)
# check_multiplicable(a.basis_r, b.basis_l);
# operators = DenseOperator[]
# for (alpha, a_alpha) in enumerate(a.operators)
# if alpha in keys(b.operators)
# push!(operators, a_alpha*b.operators[alpha])
# else
# push!(operators, deepcopy(a_alpha))
# end
# end
# if b.factor != 1.
# operators[0] *= b.factor
# end
# ProductDensityOperator(operators...)
# end

dims(bl::CompositeBasis, br::CompositeBasis) = [length(bl_i)*length(br_i) for (bl_i, br_i) in zip(bl.bases, br.bases)]
dims(x::ProductDensityOperator) = dims(x.basis_l, x.basis_r)
Expand All @@ -79,7 +79,9 @@ end
# Ignores the factor in the LazyTensor
function operators.gemm!(alpha, a::LazyTensor, b::ProductDensityOperator, beta, result::ProductDensityOperator)
@assert abs(beta)==0.
for (k, a_k) in a.operators
for n in 1:length(a.indices)
k = a.indices[n]
a_k = a.operators[n]
operators.gemm!(complex(1.), a_k, b.operators[k], complex(0.), result.operators[k])
end
end
Expand All @@ -91,9 +93,11 @@ function dmaster(rho0::ProductDensityOperator, H::LazySum,
for h_k in H.operators
operators.gemm!(1., h_k, rho0, 0., tmp)
subtraces = traces(tmp)
for (alpha, h_k_alpha) in h_k.operators
for n = 1:length(h_k.indices)
alpha = h_k.indices[n]
h_k_alpha = h_k.operators[n]
factor = h_k.factor
for gamma in keys(h_k.operators)
for gamma in h_k.indices
if alpha!=gamma
factor *= subtraces[gamma]
end
Expand All @@ -106,21 +110,23 @@ function dmaster(rho0::ProductDensityOperator, H::LazySum,
operators.gemm!(1., J[k], rho0, 0., tmp)
operators.gemm!(1., Jdagger[k], tmp, 0., tmp2)
subtraces = traces(tmp)
subindices = keys(J[k].operators)
for alpha in subindices
# subindices = keys(J[k].operators)
for n in 1:length(J[k].indices)
# for alpha in subindices
alpha = J[k].indices[n]
factor = Jdagger[k].factor * J[k].factor
for gamma in subindices
for gamma in J[k].indices
if alpha!=gamma
factor *= subtraces[gamma]
end
end
operators.gemm!(complex(factor), J[k].operators[alpha], rho0.operators[alpha], complex(0.), tmp.operators[alpha])
operators.gemm!(complex(1.), tmp.operators[alpha], Jdagger[k].operators[alpha], complex(1.), drho.operators[alpha])
operators.gemm!(complex(factor), J[k].operators[n], rho0.operators[alpha], complex(0.), tmp.operators[alpha])
operators.gemm!(complex(1.), tmp.operators[alpha], Jdagger[k].operators[n], complex(1.), drho.operators[alpha])

operators.gemm!(complex(-0.5), Jdagger[k].operators[alpha], tmp.operators[alpha], complex(1.), drho.operators[alpha])
operators.gemm!(complex(-0.5), Jdagger[k].operators[n], tmp.operators[alpha], complex(1.), drho.operators[alpha])

operators.gemm!(complex(-0.5*factor), rho0.operators[alpha], Jdagger[k].operators[alpha], complex(0.), tmp.operators[alpha])
operators.gemm!(complex(1.), tmp.operators[alpha], J[k].operators[alpha], complex(1.), drho.operators[alpha])
operators.gemm!(complex(-0.5*factor), rho0.operators[alpha], Jdagger[k].operators[n], complex(0.), tmp.operators[alpha])
operators.gemm!(complex(1.), tmp.operators[alpha], J[k].operators[n], complex(1.), drho.operators[alpha])
end
end
end
Expand Down
2 changes: 1 addition & 1 deletion src/operators_lazy.jl
Expand Up @@ -5,7 +5,7 @@ import ..operators: dagger, identityoperator,
trace, ptrace, normalize!, tensor, permutesystems

using Base.Cartesian
using ..bases, ..states, ..operators, ..operators_dense, ..operators_sparse
using ..bases, ..states, ..operators, ..operators_dense, ..operators_sparse, ..sortedindices

export lazy, LazyWrapper, LazyTensor, LazySum, LazyProduct

Expand Down

0 comments on commit 4bc626d

Please sign in to comment.