Advanced Computing Platform for Theoretical Physics

Commits (6)
......@@ -8,14 +8,15 @@ using Random; Random.seed!()
using CUDA; CUDA.allowscalar(false)
using OMEinsum
@reexport import LinearAlgebra:Eigen, diagm
@reexport import LinearAlgebra: Eigen, diagm
using PhysModels
#Operator.jl
export Operator
#solver.jl
export Device, CPU, GPU,
export Processor, CPU, GPU,
solver_function,
cpu_solver, gpu_solver
......@@ -36,24 +37,13 @@ export icgs,
basis_generate
#GeneralFTLM.jl
export Full_ED,
FTLMOptions,
export FTLMOptions, TraceEstimator,
Full_ED,
simple_FTLM,
replaced_FTLM,
orthogonalized_FTLM,
FullSampling_FTLM
#TraceEstimator.jl
export TraceEstimator #, evaluate
#Thermaldynamic.jl
#export partitian,
# free_energy,
# energy,
# specific_heat,
# entropy,
# thermal_average
include("Operator.jl")
......@@ -63,7 +53,9 @@ include("eigensolver.jl")
include("FullOrthoLanczos.jl")
include("GeneralFTLM.jl")
include("TraceEstimator.jl")
include("GeneralFTLMForAverage.jl")
include("evaluate.jl")
#PhysModels
include("Thermaldynamic.jl")
......
......@@ -3,8 +3,7 @@
GeneralFTLM Methods(alse konwn as trace estimators)
For a given method, denoted as `f`,
`f(M, expr)` is to estimate `tr(expr(M))`,
`f(M, expr, Op)` is to evaluate `tr(Op expr(M))`,
where `M` and `Op` are matrices or matrix-like instances,
where `M` is a matrix or matrix-like instances,
and `M` is required to be hermitian.
Note:
......@@ -30,14 +29,23 @@ end
where `:SR` means smallest real part and `:LR` means largest real part.
"""
@with_kw struct FTLMOptions{Tdistr<:Function,
Tdevice<:Device,
Tprocessor<:Processor,
Tn<:Integer}
distr::Tdistr = Rademacher
Nr::Tn = 100
Nk::Tn = 50
Ne::Tn = 10
which::Symbol = :SR
device::Tdevice = CPU
processor::Tprocessor = CPU
end
"""
TraceEstimator struct
"""
@with_kw struct TraceEstimator{Tf<:Function, To<:FTLMOptions}
estimator::Tf
options::To
end
......@@ -47,28 +55,14 @@ end
function Full_ED(M, expr::Vararg{T,N};
options::FTLMOptions = FTLMOptions()
) where {T<:Function, N}
@unpack device = options
solver = solver_function(device)
@unpack processor = options
solver = solver_function(processor)
M = solver(x->x, M)
val, _ = eigensolver(M)
func = f -> map(e->f(e), val) |> sum
vals, _ = eigensolver(M)
func = f -> reduce(+, map(f , vals))
return map(func, expr)
end
function Full_ED(M, Op::Vararg{T,N};
expr::Te,
options::FTLMOptions=FTLMOptions()
) where {T<:Operator, N, Te<:Function}
@unpack device = options
solver = solver_function(device)
M = solver(x->x, M)
val, vec = eigensolver(M)
Z = map(expr, val) |> sum
f = x -> map((i -> expr(val[i]) * vec[:,i]' * x.O * vec[:,i]), 1:length(val)) |> sum
res = map(f, Op)
return map(x -> x/Z, res)
end
"""
Simple FTLM algorithom, where `M` is required to be hermitian.
......@@ -76,8 +70,8 @@ end
function simple_FTLM(M, expr::Vararg{T,N};
options::FTLMOptions = FTLMOptions()
) where {T<:Function, N}
@unpack distr, Nr, Nk, device = options
solver = solver_function(device)
@unpack distr, Nr, Nk, processor = options
solver = solver_function(processor)
M = solver(x->x, M)
Ns = size(M, 1)
......@@ -86,38 +80,16 @@ function simple_FTLM(M, expr::Vararg{T,N};
v0 = random_unit_vector(Ns, distr)
v0 = solver(x->x, v0)
@unpack values, weight = itFOLM(M, init_vector = v0, Nk = Nk) |> eigensolver
func = f -> map((e,w)->f(e)* w * w', values, weight) |> sum
func = f -> begin
Λ = map(f, values)
Z = ein"i,i,i -> "(Λ, weight, conj(weight))
return Array(Z)[1]
end
res = map(+, res, map(func, expr))
end
return map(x -> x * Ns/Nr, res)
end
function simple_FTLM(M, Op::Vararg{T,N};
expr::Te,
options::FTLMOptions = FTLMOptions()
) where {T<:Operator,N,Te<:Function}
@unpack distr, Nr, Nk, device = options
solver = solver_function(device)
M = solver(x->x, M)
Ns = size(M,1)
res = zeros(N)
Z = 0.0
for r = 1: Nr
v0 = random_unit_vector(Ns, distr)
v0 = solver(x->x, v0)
@unpack weight, values, vectors = itFOLM(M, init_vector = v0, Nk = Nk) |> eigensolver
Z += map((e,w)-> expr(e) * w * w', values, weight) |> sum
operator_average = x -> begin
weight2 = map(i -> v0' * x.O * vectors[:,i], 1:Nk)
return map((e, w1, w2)-> expr(e) * w1 * w2, values, weight, weight2) |> sum
end
res = map(+, res, map(operator_average, Op))
end
Z = Z * Ns/Nr
return map(x -> x/Z * Ns/Nr, res)
end
"""
Replaced FTLM algorithom
......@@ -125,141 +97,77 @@ end
function replaced_FTLM(M, expr::Vararg{T,N};
options::FTLMOptions = FTLMOptions()
) where {T<:Function,N}
@unpack distr, Nr, Nk, Ne, device, which = options
solver = solver_function(device)
@unpack distr, Nr, Nk, Ne, processor, which = options
solver = solver_function(processor)
M = solver(x->x, M)
Ns = size(M, 1)
res = zeros(N)
Ne = min(Ne, Ns)
krylovdim = max(30, Ne+1)
evals, evecs, _ = eigsolve(M, Ns, Ne, which, ishermitian = true, krylovdim=krylovdim)
evals = evals[1:Ne]
evecs = hcat(evecs[1:Ne]...)
processor == CPU ? x0 = rand(Ns) : x0 = CUDA.rand(Ns)
vals, vecs, _ = eigsolve(M, x0, Ne, which, ishermitian = true, krylovdim=krylovdim)
eigen_vals = solver(x->x, vals[1:Ne])
eigen_vecs = hcat(vecs[1:Ne]...)
Nk = max(Nk, Ne)
for r = 1: Nr
v0 = random_unit_vector(Ns, distr)
v0 = solver(x->x, v0)
@unpack weight, values = itFOLM(M, init_vector = v0, Nk = Nk) |> eigensolver
values = vcat(evals, values[Ne+1:end])
eweight = map(i-> v0' * evecs[:,i], 1:Ne)
weight = vcat(eweight, weight[Ne+1:end])
func = f -> map((e,w)->f(e)* w * w', values, weight) |> sum
res = map(+, res, map(func, expr))
end
return map(x -> x * Ns/Nr, res)
end
function replaced_FTLM(M, Op::Vararg{T,N};
expr::Te,
options::FTLMOptions=FTLMOptions()
) where {T<:Operator,N,Te<:Function}
@unpack distr, Nr, Nk, Ne, device, which = options
solver = solver_function(device)
M = solver(x->x, M)
eigen_weight = ein"i,ij->j"(conj(v0), eigen_vecs)
weight = vcat(eigen_weight, weight[Ne+1:end])
values = vcat(eigen_vals, values[Ne+1:end])
Ns = size(M, 1)
res = zeros(N)
Ne = min(Ne, Ns)
krylovdim = max(30, Ne+1)
evals, evecs, _ = eigsolve(M, Ns, Ne, which, ishermitian = true, krylovdim=krylovdim)
evals = evals[1:Ne]
evecs = hcat(evecs[1:Ne]...)
Z = 0.0
Nk = max(Nk, Ne)
for r = 1: Nr
v0 = random_unit_vector(Ns, distr)
v0 = solver(x->x, v0)
@unpack weight, values, vectors = itFOLM(M, init_vector = v0, Nk = Nk) |> eigensolver
values = vcat(evals, values[Ne+1:end])
vectors = hcat(evecs, vectors[:, Ne+1:end])
eweight = map(i-> v0' * evecs[:,i], 1:Ne)
weight = vcat(eweight, weight[Ne+1:end])
Z += map((e,w)-> expr(e) * w * w', values, weight) |> sum
operator_average = x -> begin
weight2 = map(i -> v0' * x.O * vectors[:,i], 1:Nk)
return map((e, w1, w2)-> expr(e)* w1 * w2, val, weight, weight2) |> sum
func = f -> begin
Λ = map(f, values)
Z = ein"i,i,i -> "(Λ, weight, conj(weight))
return Array(Z)[1]
end
res = map(+, res, map(operator_average, Op))
res = map(+, res, map(func, expr))
end
Z = Z * Ns/Nr
return map(x -> x/Z * Ns/Nr, res)
factor = Ns/Nr
return map(x -> x * factor, res)
end
"""
Orthogonalized FTLM algorithom
`M`: Symmetric matrix
"""
function orthogonalized_FTLM(M, expr::Vararg{T,N};
options::FTLMOptions=FTLMOptions()
) where {T<:Function,N}
@unpack distr, Nr, Nk, Ne, device, which = options
solver = solver_function(device)
@unpack distr, Nr, Nk, Ne, processor, which = options
solver = solver_function(processor)
M = solver(x->x, M)
Ns = size(M, 1)
Ne = min(Ne, Ns)
krylovdim = max(30, Ne+1)
evals, evecs, _ = eigsolve(M, Ns, Ne, which, ishermitian = true, krylovdim=krylovdim)
evals = evals[1:Ne]
evecs = hcat(evecs[1:Ne]...)
processor == CPU ? x0 = rand(Ns) : x0 = CUDA.rand(Ns)
vals, vecs, _ = eigsolve(M, x0, Ne, which, ishermitian = true, krylovdim=krylovdim)
eigen_vals = vals[1:Ne] #always be Vector type, not CuVector
eigen_vecs = hcat(vecs[1:Ne]...)
func1 = f -> map(e->f(e), evals) |> sum
func1 = f -> reduce(+, map(f, eigen_vals))
res1 = map(func1, expr)
res2 = zeros(N)
for r = 1: Nr
v0 = random_unit_vector(Ns, distr)
v0 = solver(x->x, v0)
@unpack weight, values = itFOLM(M, evecs, init_vector = v0, Nk = Nk) |> eigensolver
func2 = f -> map((e,w)->f(e)* w * w', values, weight) |> sum
res2 = map(+, res2, map(func2, expr))
end
return map((x1, x2) -> x1 + x2 * (Ns-Ne)/Nr, res1, res2)
end
function orthogonalized_FTLM(M, Op::Vararg{T,N};
expr::Te,
options::FTLMOptions=FTLMOptions()
) where {T<:Operator,N,Te<:Function}
@unpack distr, Nr, Nk, Ne, device, which = options
solver = solver_function(device)
M = solver(x->x, M)
Ns = size(M, 1)
Ne = min(Ne, Ns)
krylovdim = max(30, Ne+1)
evals, evecs, _ = eigsolve(M, Ns, Ne, which, ishermitian = true, krylovdim=krylovdim)
evals = evals[1:Ne]
evecs = hcat(evecs[1:Ne]...)
@unpack weight, values = itFOLM(M, eigen_vecs, init_vector = v0, Nk = Nk) |> eigensolver
Z1 = map(func, evals) |> sum
func1 = x -> map((i -> expr(evals[i]) * evecs[:,i]' * x.O * evecs[:,i]), 1:length(evals)) |> sum
res1 = map(func1, Op)
Z2 = 0.0
res2 = zeros(N)
for r = 1: Nr
v0 = random_unit_vector(Ns, distr)
v0 = solver(x->x, v0)
@unpack weight, values, vectors = itFOLM(M, evecs, init_vector = v0, Nk = Nk) |> eigensolver
Z2 += map((e,w)-> expr(e) * w * w', values, weight) |> sum
operator_average = x -> begin
weight2 = map(i -> v0' * x.O * vectors[:,i], 1:Nk)
return map((e, w1, w2)-> expr(e)* w1 * w2, val, weight, weight2) |> sum
func2 = f -> begin
Λ = map(f, values)
Z = ein"i,i,i -> "(Λ, weight, conj(weight))
return Array(Z)[1]
end
res2 = map(+, res2, map(operator_average, Op))
res2 = map(+, res2, map(func2, expr))
end
Z = Z1 + Z2 * (Ns-Ne)/Nr
return map((x1, x2) -> (x1 + x2*(Ns-Ne)/Nr)/Z, res1, res2)
factor = (Ns-Ne)/Nr
return map((x1, x2) -> x1 + x2 * factor, res1, res2)
end
......@@ -271,8 +179,8 @@ function FullSampling_FTLM(M, expr::Vararg{T,N};
) where {T<:Function,N}
Ns = size(M, 1)
ortho_basis = basis_generate(Ns)
@unpack Nk, which, device = options
solver = solver_function(device)
@unpack Nk, which, processor = options
solver = solver_function(processor)
M = solver(x->x, M)
#Full FullSampling_FTLM
......@@ -281,36 +189,15 @@ function FullSampling_FTLM(M, expr::Vararg{T,N};
v0 = ortho_basis[:, i]
v0 = solver(x->x, v0)
@unpack weight, values = itFOLM(M, init_vector = v0, Nk = Nk) |> eigensolver
func2 = f -> map((e,w)->f(e)* w * w', values, weight) |> sum
res = map(+, res, map(func2, expr))
end
return res
end
function FullSampling_FTLM(M, Op::Vararg{T,N};
expr::Te,
options::FTLMOptions = FTLMOptions()
) where {T<:Operator,N,Te<:Function}
Ns = size(M, 1)
ortho_basis = basis_generate(Ns)
@unpack Nk, which, device = options
solver = solver_function(device)
M = solver(x->x, M)
Z = 0.0
res = zeros(N)
for i = 1: Ns
v0 = ortho_basis[:, i]
v0 = solver(x->x, v0)
@unpack weight, values, vectors = itFOLM(M, init_vector = v0, Nk = Nk) |> eigensolver
Z += map((e,w)-> expr(e) * w * w', values, weight) |> sum
operator_average = x -> begin
weight2 = map(i -> v0' * x.O * vectors[:,i], 1:Nk)
return map((e, w1, w2)-> expr(e)* w1 * w2, val, weight, weight2) |> sum
func = f -> begin
Λ = map(f, values)
Z = ein"i,i,i -> "(Λ, weight, conj(weight))
return Array(Z)[1]
end
res = map(+, res, map(operator_average, Op))
res = map(+, res, map(func, expr))
end
return map(x -> x/Z, res)
return res
end
#end #module methods
#module methods
#TODO: use Einsum instead of elementwise operations
#CUARRAY 会有 allowscalor 的问题
"""
GeneralFTLM Methods(alse konwn as trace estimators)
For a given method, denoted as `f`,
`f(M, expr, Op)` is to evaluate `tr(Op expr(M))`,
where `M` and `Op` are matrices or matrix-like instances,
and `M` is required to be hermitian.
Note:
1. A matrix-like instance `M` behaves as a matrix,
whose `eltype`, `size`, `length`, `getindex`, `iterate`
and matrix-vector products are required.
2. `M` is not explicity symmetrized in the following methods.
Full exact diagonalization method
"""
function Full_ED(M, Op::Vararg{T,N};
expr::Te,
options::FTLMOptions=FTLMOptions()
) where {T<:Operator, N, Te<:Function}
@unpack processor = options
solver = solver_function(processor)
M = solver(x->x, M)
vals, vecs = eigensolver(M)
Λ = map(expr, vals)
Z = reduce(+, Λ)
func = x -> begin
ave = ein"i,ik,kl,li->"(Λ, conj(vecs), x.O, vecs)
return Array(ave)[1]
end
res = map(func, Op)
return map(x -> x/Z, res)
end
"""
Simple FTLM algorithom, where `M` is required to be hermitian.
"""
function simple_FTLM(M, Op::Vararg{T,N};
expr::Te,
options::FTLMOptions = FTLMOptions()
) where {T<:Operator,N,Te<:Function}
@unpack distr, Nr, Nk, processor = options
solver = solver_function(processor)
M = solver(x->x, M)
Ns = size(M,1)
res = zeros(N)
Z = 0.0
for r = 1: Nr
v0 = random_unit_vector(Ns, distr)
v0 = solver(x->x, v0)
@unpack weight, values, vectors = itFOLM(M, init_vector = v0, Nk = Nk) |> eigensolver
Λ = map(expr, vals)
Z += Array(ein"i,i,i -> "(Λ, weight, conj(weight)))[1]
func = x -> begin
weight2 = ein"i,ik,kj->j"(conj(v0), x.O, vectors)
ave = ein"i,i,i->"(Λ, weight, weight2)[1]
return Array(ave)[1]
end
res = map(+, res, map(func, Op))
end
Z = Z * Ns/Nr
return map(x -> x/Z * Ns/Nr, res)
end
"""
Replaced FTLM algorithom
"""
function replaced_FTLM(M, Op::Vararg{T,N};
expr::Te,
options::FTLMOptions=FTLMOptions()
) where {T<:Operator,N,Te<:Function}
@unpack distr, Nr, Nk, Ne, processor, which = options
solver = solver_function(processor)
M = solver(x->x, M)
Ns = size(M, 1)
res = zeros(N)
Ne = min(Ne, Ns)
krylovdim = max(30, Ne+1)
processor == CPU ? x0 = rand(Ns) : x0 = CUDA.rand(Ns)
vals, vecs, _ = eigsolve(M, x0, Ne, which, ishermitian = true, krylovdim=krylovdim)
eigen_vals = solver(x->x, vals[1:Ne])
eigen_vecs = hcat(vecs[1:Ne]...)
Z = 0.0
Nk = max(Nk, Ne)
for r = 1: Nr
v0 = random_unit_vector(Ns, distr)
v0 = solver(x->x, v0)
@unpack weight, values, vectors = itFOLM(M, init_vector = v0, Nk = Nk) |> eigensolver
eigen_weight = ein"i,ij->j"(conj(v0), eigen_vecs)
weight = vcat(eigen_weight, weight[Ne+1:end])
values = vcat(eigen_vals, values[Ne+1:end])
vectors = hcat(eigen_vecs, vectors[:, Ne+1:end])
Λ = map(expr, values)
Z += Array(ein"i,i,i -> "(Λ, weight, conj(weight)))[1]
func = x -> begin
weight2 = ein"i,ik,kj->j"(conj(v0), x.O, vectors)
ave = ein"i,i,i->"(Λ, weight, weight2)[1]
return Array(ave)[1]
end
res = map(+, res, map(func, Op))
end
factor = Ns/Nr
Z = Z * factor
return map(x -> x * factor/Z, res)
end
"""
Orthogonalized FTLM algorithom
"""
function orthogonalized_FTLM(M, Op::Vararg{T,N};
expr::Te,
options::FTLMOptions=FTLMOptions()
) where {T<:Operator,N,Te<:Function}
@unpack distr, Nr, Nk, Ne, processor, which = options
solver = solver_function(processor)
M = solver(x->x, M)
Ns = size(M, 1)
Ne = min(Ne, Ns)
krylovdim = max(30, Ne+1)
processor == CPU ? x0 = rand(Ns) : x0 = CUDA.rand(Ns)
vals, vecs, _ = eigsolve(M, x0, Ne, which, ishermitian = true, krylovdim=krylovdim)
eigen_vals = vals[1:Ne]
eigen_vecs = hcat(vecs[1:Ne]...)
Λ1 = map(expr, eigen_vals)
Z1 = reduce(+, Λ1)
func1 = x -> begin
ave = ein"i,ik,kl,li->"(Λ1, conj(eigen_vecs), x.O, eigen_vecs)
return Array(ave)[1]
end
res1 = map(func1, Op)
Z2 = 0.0
res2 = zeros(N)
for r = 1: Nr
v0 = random_unit_vector(Ns, distr)
v0 = solver(x->x, v0)
@unpack weight, values, vectors = itFOLM(M, eigen_vecs, init_vector = v0, Nk = Nk) |> eigensolver
Λ2 = map(expr, values)
Z2 += ein"i,i,i -> "(Λ2, weight, conj(weight))[1]
func2 = x -> begin
weight2 = ein"i,ik,kj->j"(conj(v0), x.O, vectors)
ave = ein"i,i,i->"(Λ, weight, weight2)[1]
return Array(ave)[1]
end
res2 = map(+, res2, map(func2, Op))
end
factor = (Ns-Ne)/Nr
Z = Z1 + Z2 * factor
return map((x1, x2) -> (x1 + x2 * factor)/Z, res1, res2)
end
"""
Full sampling simple_FTLM
"""
function FullSampling_FTLM(M, Op::Vararg{T,N};
expr::Te,
options::FTLMOptions = FTLMOptions()
) where {T<:Operator,N,Te<:Function}
Ns = size(M, 1)
ortho_basis = basis_generate(Ns)
@unpack Nk, which, processor = options
solver = solver_function(processor)
M = solver(x->x, M)
Z = 0.0
res = zeros(N)
for i = 1: Ns
v0 = ortho_basis[:, i]
v0 = solver(x->x, v0)
@unpack weight, values, vectors = itFOLM(M, init_vector = v0, Nk = Nk) |> eigensolver
Λ = map(expr, values)
Z += Array(ein"i,i,i -> "(Λ, weight, conj(weight)))[1]
func = x -> begin
weight2 = ein"i,ik,kj->j"(conj(v0), x.O, vectors)
ave = ein"i,i,i->"(Λ, weight, weight2)[1]
return Array(ave)[1]
end
res = map(+, res, map(func, Op))
end
return map(x -> x/Z, res)
end
#end #module methods
"""
partitian function `Z = ∑exp(-β(E-E0))`
"""
function partitian(M, β; estimator::TraceEstimator)
e0, _, _ = eigsolve(M, size(M,1), 1, :SR, ishermitian = true)
e0 = e0[1]
exprZ = e -> exp(-β*(e - e0))
function partitian(M, β; trace_estimator::TraceEstimator)
@unpack processor = trace_estimator.options
processor == CPU ? x0 = rand(size(M,1)) : x0 = CUDA.rand(size(M,1))
e0, _, _ = eigsolve(M, x0, 1, :SR, ishermitian = true)
e1 = e0[1]
exprZ = e -> exp(-β*(e - e1))
return evaluate(estimator, M, exprZ)[1]
end
......@@ -13,10 +15,12 @@ end
total free energy `F = -lnZ/β`,
use `F/L` to calculate free energy per site
"""
function free_energy(M, β; estimator::TraceEstimator)
e0, _, _ = eigsolve(M, size(M,1), 1, :SR, ishermitian = true)
e0 = e0[1]
exprZ = e -> exp(-β*(e - e0))
function free_energy(M, β; trace_estimator::TraceEstimator)
@unpack processor = trace_estimator.options
processor == CPU ? x0 = rand(size(M,1)) : x0 = CUDA.rand(size(M,1))
e0, _, _ = eigsolve(M, x0, 1, :SR, ishermitian = true)
e1 = e0[1]
exprZ = e -> exp(-β*(e - e1))
res = evaluate(estimator, M, exprZ)[1]
return -log(res)/β + e0
end
......@@ -26,11 +30,13 @@ end
total energy `E`,
use `E/L` to calculate energy per site
"""
function energy(M, β; estimator::TraceEstimator)
e0, _, _ = eigsolve(M, size(M,1), 1, :SR, ishermitian = true)
e0 = e0[1]
exprZ = e -> exp(-β*(e - e0))
exprE = e -> e * exp(-β*(e - e0))
function energy(M, β; trace_estimator::TraceEstimator)
@unpack processor = trace_estimator.options
processor == CPU ? x0 = rand(size(M,1)) : x0 = CUDA.rand(size(M,1))
e0, _, _ = eigsolve(M, x0, 1, :SR, ishermitian = true)
e1 = e0[1]
exprZ = e -> exp(-β*(e - e1))
exprE = e -> e * exp(-β*(e - e1))
Z, E = evaluate(estimator, M, exprZ, exprE)
return E/Z
end
......@@ -41,12 +47,14 @@ end
total specific heat `Cv`,
use `Cv/L` to calculate specific heat per site
"""
function specific_heat(M, β; estimator::TraceEstimator)
e0, _, _ = eigsolve(M, size(M,1), 1, :SR, ishermitian = true)
e0 = e0[1]
exprZ = e -> exp(-β*(e - e0))
exprE = e -> e * exp(-β*(e - e0))
exprCv = e -> e^2 * exp(-β *(e - e0))
function specific_heat(M, β; trace_estimator::TraceEstimator)
@unpack processor = trace_estimator.options
processor == CPU ? x0 = rand(size(M,1)) : x0 = CUDA.rand(size(M,1))
e0, _, _ = eigsolve(M, x0, 1, :SR, ishermitian = true)
e1 = e0[1]
exprZ = e -> exp(-β*(e - e1))
exprE = e -> e * exp(-β*(e - e1))
exprCv = e -> e^2 * exp(-β *(e - e1))
Z, E, Cv = evaluate(estimator, M, exprZ, exprE, exprCv)
E = E/Z
......@@ -58,11 +66,13 @@ end
total entropy `S = β(E-F)`,
use `S/L` to calculate entropy per site
"""
function entropy(M, β; estimator::TraceEstimator)
e0, _, _ = eigsolve(M, size(M,1), 1, :SR, ishermitian = true)
e0 = e0[1]
exprZ = e -> exp(-β*(e - e0))
exprE = e -> e * exp(-β*(e - e0))
function entropy(M, β; trace_estimator::TraceEstimator)
@unpack processor = trace_estimator.options
processor == CPU ? x0 = rand(size(M,1)) : x0 = CUDA.rand(size(M,1))
e0, _, _ = eigsolve(M, x0, 1, :SR, ishermitian = true)
e1 = e0[1]
exprZ = e -> exp(-β*(e - e1))
exprE = e -> e * exp(-β*(e - e1))
Z, E = evaluate(estimator, M, exprZ, exprE)
F = -log(Z)/β + e0
E = E/Z
......@@ -74,10 +84,12 @@ end
thermal average `⟨O⟩ = Tr(exp(-βH)O)/Z` of a operator `O`
"""
function thermal_average(M, β, Op::Vararg{T,N};
estimator::TraceEstimator) where {T<:Operator, N}
e0, _, _ = eigsolve(M, size(M,1), 1, :SR, ishermitian = true)
e0 = e0[1]
expr = e -> exp(-β*(e - e0))
trace_estimator::TraceEstimator) where {T<:Operator, N}
@unpack processor = trace_estimator.options
processor == CPU ? x0 = rand(size(M,1)) : x0 = CUDA.rand(size(M,1))
e0, _, _ = eigsolve(M, x0, 1, :SR, ishermitian = true)
e1 = e0[1]
expr = e -> exp(-β*(e - e1))
return evaluate(estimator, M, Op, expr = expr)
end
import LinearAlgebra: Eigen
"""
`eigensolver`: functions to solve all eigenpairs for a instance `M` which behaves as a matrix
"""
......
"""
TraceEstimator struct
evaluate function
"""
@with_kw struct TraceEstimator{Tf<:Function, To<:FTLMOptions}
estimator::Tf
options::To
end
function evaluate(method::TraceEstimator, M, expr::Vararg{T,N}) where {T<:Function,N}
@unpack estimator, options = method
return estimator(M, expr..., options=options)
......@@ -15,7 +9,4 @@ end
function evaluate(method::TraceEstimator, M, Op::Vararg{T,N}; expr::Te) where {T<:Operator,N,Te<:Function}
@unpack estimator, options = method
return estimator(M, Op..., expr = expr, options=method.options)
end
\ No newline at end of file
end
\ No newline at end of file
@enum Device CPU GPU
@enum Processor CPU GPU
function cpu_solver end
function gpu_solver end
function solver_function(d::Device)
d == CPU ? (return cpu_solver) : (return cpu_solver)
function solver_function(processor::Processor)
processor == CPU ? (return cpu_solver) : (return gpu_solver)
end
function cpu_solver(f::Function, M::AbstractArray...)
......
using Test, FiniteTLanczos
using LinearAlgebra, Parameters
using Random; Random.seed!()
solver = solver_function(device)
@testset "icgs algorithm" begin
N = 256
......
using Test, FiniteTLanczos
using LinearAlgebra, KrylovKit, LogExpFunctions, Parameters
using Random; Random.seed!()
solver = solver_function(device)
function logtrexp(M::AbstractMatrix)
vals, _ = eigensolver(M)
return logsumexp(vals)
end
function logtrexp(M, estimator::TraceEstimator)
e0, _, _ = eigsolve(M, size(M,1), 1, :LR, ishermitian = true)
e0 = e0[1]
expr = e -> exp(e .- e0)
@unpack estimator, options = estimator
function logtrexp(M, trace_estimator::TraceEstimator)
@unpack estimator, options = trace_estimator
@unpack processor = options
options = FTLMOptions(options, which = :LR)
processor == CPU ? x0 = rand(size(M,1)) : x0 = CUDA.rand(size(M,1))
e0, _, _ = eigsolve(M, x0, 1, :LR, ishermitian = true)
e1 = e0[1]
expr = e -> exp(e - e1)
res = FiniteTLanczos.evaluate(TraceEstimator(estimator, options), M, expr)[1]
return log(res) + e0
return log(res) + e1
end
N = 256
......@@ -25,21 +27,21 @@ M = solver(x->x, M)
log_trexp = logtrexp(M)
@testset "Full_ED method" begin
estimator = TraceEstimator(Full_ED, FTLMOptions(FTLMOptions(), device = device))
res = logtrexp(M, estimator)
trace_estimator = TraceEstimator(Full_ED, FTLMOptions(FTLMOptions(), processor = processor))
res = logtrexp(M, trace_estimator)
@test res log_trexp
end
@testset "simple_FTLM method" begin
Ns = size(M,1)
options = FTLMOptions(FTLMOptions(), Nk=Ns, device = device)
estimator = TraceEstimator(FullSampling_FTLM, options)
res = logtrexp(M, estimator)
options = FTLMOptions(FTLMOptions(), Nk=Ns, processor = processor)
trace_estimator = TraceEstimator(FullSampling_FTLM, options)
res = logtrexp(M, trace_estimator)
@test res log_trexp
options = FTLMOptions(FTLMOptions(), device = device)
estimator = TraceEstimator(simple_FTLM, options)
res = logtrexp(M, estimator)
options = FTLMOptions(FTLMOptions(), processor = processor)
trace_estimator = TraceEstimator(simple_FTLM, options)
res = logtrexp(M, trace_estimator)
@unpack Nr = options
rtol = 1/sqrt(min(Nr, Ns))
......@@ -48,15 +50,15 @@ end
@testset "orthogonalized_FTLM method" begin
Ns = size(M,1)
options = FTLMOptions(FTLMOptions(), Ne = Ns , device = device)
estimator = TraceEstimator(orthogonalized_FTLM, options)
res = logtrexp(M, estimator)
@test res log_trexp
options = FTLMOptions(FTLMOptions(), Ne = Ns , processor = processor)
trace_estimator = TraceEstimator(orthogonalized_FTLM, options)
res = logtrexp(M, trace_estimator)
@test (res, log_trexp, rtol = 1.e-5)
for Ne in [1, div(Ns,3), div(Ns,2), Ns-1]
options = FTLMOptions(FTLMOptions(), Ne = Ne, device = device)
estimator = TraceEstimator(orthogonalized_FTLM, options)
res = logtrexp(M, estimator)
options = FTLMOptions(FTLMOptions(), Ne = Ne, processor = processor)
trace_estimator = TraceEstimator(orthogonalized_FTLM, options)
res = logtrexp(M, trace_estimator)
@unpack Nr = options
rtol = 1/Nr
......@@ -67,15 +69,15 @@ end
@testset "replaced_FTLM method" begin
Ns = size(M,1)
options = FTLMOptions(FTLMOptions(), Ne = Ns , device = device)
estimator = TraceEstimator(replaced_FTLM, options)
res = logtrexp(M, estimator)
options = FTLMOptions(FTLMOptions(), Ne = Ns , processor = processor)
trace_estimator = TraceEstimator(replaced_FTLM, options)
res = logtrexp(M, trace_estimator)
@test (res, log_trexp, rtol = 0.1)
for Ne in [1, div(Ns,3), div(Ns,2), Ns-1]
options = FTLMOptions(FTLMOptions(), Ne = Ne, device = device)
estimator = TraceEstimator(replaced_FTLM, options)
res = logtrexp(M, estimator)
options = FTLMOptions(FTLMOptions(), Ne = Ne, processor = processor)
trace_estimator = TraceEstimator(replaced_FTLM, options)
res = logtrexp(M, trace_estimator)
@unpack Nr = options
rtol = 1/sqrt(min(Nr, Ns))
......
using Test, FiniteTLanczos
device = CPU
processor = CPU
solver = solver_function(processor)
@show solver
@testset "FullOrthoLanczos.jl" begin
include("FullOrthoLanczos.jl")
......
using CUDA; CUDA.allowscalar(false)
using Test, FiniteTLanczos
device = GPU
processor = GPU
solver = solver_function(processor)
@show solver
@testset "FullOrthoLanczos.jl" begin
include("FullOrthoLanczos.jl")
......