Advanced Computing Platform for Theoretical Physics

Skip to content

GitLab

Commits (23)
Hide whitespace changes
Inline Side-by-side
src/CMPSCompress.jl
View file @ 2f2350ff
......@@ -4,11 +4,12 @@
function compress_cmps(ψ0::AbstractCMPS{T, S, U}, χ::Integer, β::Real;
options::CompressOptions=CompressOptions()) where {T,S,U}
@unpack (init, show_trace, mera_update_options,
optim_options, trace_estimator) = options
optim_options, trace_estimator, processor) = options
if show_trace
println("----------------------------Compress CMPS-----------------------------")
end
solver = solver_function(processor)
χ0 = size(ψ0.Q, 1)
U <: AbstractMatrix ? vir_dim = 1 : vir_dim = size(ψ0.R, 3)
optim_result = nothing
......
src/CMPSOperations.jl
View file @ 2f2350ff
......@@ -105,6 +105,14 @@ function ==(a::T, b::T) where T<:AbstractCMPS
end
"""
==(a::T, b::T) where T<:CMPSMatrix
"""
function ==(a::T, b::T) where T<:CMPSMatrix
return ==(a.ψl, b.ψl) && ==(a.ψr, b.ψr)
end
"""
≈(a::T, b::T) where T<:AbstractCTensor
"""
......@@ -116,7 +124,16 @@ function ≈(a::T, b::T; kwargs...) where T<:AbstractCMPO
end
function ≈(a::T, b::T; kwargs...) where T<:AbstractCMPS
return ≈(a.Q, b.Q; kwargs...) && ≈(a.R, b.R)
return ≈(a.Q, b.Q; kwargs...) &&
≈(a.R, b.R; kwargs...)
end
"""
≈(a::T, b::T) where T<:CMPSMatrix
"""
function ≈(a::T, b::T; kwargs...) where T<:CMPSMatrix
return ≈(a.ψl, b.ψl; kwargs...) && ≈(a.ψr, b.ψr; kwargs...)
end
......
src/JuliaCMPO.jl
View file @ 2f2350ff
......@@ -12,9 +12,13 @@ using PhysModels, FiniteTLanczos
@reexport import Base: kron, *, eltype, size, length, getindex, iterate
@reexport import Base: ==, ≈, transpose, adjoint, cat
@reexport import LinearAlgebra: ishermitian, norm, normalize, diag, diagm
@reexport import LinearAlgebra: ishermitian, norm, normalize, diag
@reexport import KrylovKit.eigsolve
@reexport import FiniteTLanczos: eigensolver, symmetrize
@reexport import FiniteTLanczos: eigensolver, symmetrize,
Processor, CPU, GPU,
solver_function,
cpu_solver, gpu_solver,
FTLMOptions, TraceEstimator
# structs/CMPSandCMPO.jl
......@@ -35,11 +39,9 @@ export MeraUpdateOptions,
export CompressOptions, CompressResult
#utilities
export ⊗, diagm
export ⊗
#solver.jl
export cpu_solver, gpu_solver
# OptimFunctions.jl
export veclength, optim_functions
# SaveLoad.jl
......@@ -71,10 +73,11 @@ export make_operator
#structs/EvaluateOptions.jl
export EstimatorType, EvaluateOptions
#evaluate.jl
export evaluate,
hermitian_evaluate,
non_hermitian_evaluate
#export evaluate,
# hermitian_evaluate,
# non_hermitian_evaluate
include("./structs/CMPSandCMPO.jl")
......@@ -84,7 +87,6 @@ include("./structs/CMPSCompress.jl")
include("./utilities/otimes.jl")
include("./utilities/diagm.jl")
include("solver.jl")
include("OptimFunctions.jl")
......
src/eigensolver.jl
View file @ 2f2350ff
@reexport import FiniteTLanczos.eigensolver
@reexport import KrylovKit.eigsolve
import LinearAlgebra.Eigen
for elty in (:Float32, :Float64)
"""
eigensolver for CuMatrix
For CUDA dense matrix eigen solver `CUSOLVER.syevd!` and `CUSOLVER.heevd!`:
'N'/'V': return eigenvalues/both eigenvalues and eigenvectors
'U'/'L': Upper/Lower triangle of `M` is stored.
"""
@eval eigensolver(M::CuMatrix{$elty}) = CUSOLVER.syevd!('V', 'U', M)
function eigensolver(M::CMPSMatrix{Ts,T,S,U};
which::Symbol=:SR,
ishermitian = true,
kwargs...) where {Ts,T,S,U}
N = size(M,1)
Ts <: CMPS ? x0 = rand(T,N) : x0 = CUDA.rand(T, N)
vals, vecs, _ = eigsolve(M, x0, N, which, ishermitian = ishermitian, krylovdim = N, kwargs...)
vals = convert(typeof(x0), vals)
vecs = hcat(vecs[1:N]...)
return Eigen(vals, vecs)
end
for elty in (:ComplexF32, :ComplexF64)
"""
eigensolver for CuMatrix
For CUDA dense matrix eigen solver `CUSOLVER.syevd!` and `CUSOLVER.heevd!`:
'N'/'V': return eigenvalues/both eigenvalues and eigenvectors
'U'/'L': Upper/Lower triangle of `M` is stored.
"""
@eval eigensolver(M::CuArray{$elty}) = CUSOLVER.heevd!('V', 'U', M)
end
eigsolve(M::CMPSMatrix{T}, howmany, which; args...) where {T} =
eigsolve(M, size(M,1), howmany, which, T, ishermitian = true, args...)
#function eigsolve(M::CMPSMatrix{Ts,T,S,U}, x0, howmany, which; args...) where {Ts,T,S,U}
# Ts <: CMPS ? x0 = rand(size(M,1)) : x0 = CUDA.rand(size(M,1))
# return eigsolve(M, x0, howmany, which, T, ishermitian = true, args...)
#end
\ No newline at end of file
src/evaluate.jl
View file @ 2f2350ff
......@@ -20,8 +20,15 @@ end
"""
function hermitian_evaluate(m::PhysModel, bondD::Integer, β::Real, ResultFolder::String;
options::EvaluateOptions=EvaluateOptions())
@unpack (init, solver, trace_estimator,
@unpack (init, processor, trace_estimator,
compress_options, optim_options, tag) = options
if trace_estimator !== nothing
new_options = FTLMOptions(trace_estimator.options, processor=processor)
trace_estimator = TraceEstimator(trace_estimator.estimator, new_options)
compress_options = CompressOptions(compress_options, trace_estimator = trace_estimator)
end
solver = solver_function(processor)
"""
ResultFolder: classified by Model parameters(interaction, width), store CMPS and Obsv files
CMPSResultFolder: CMPS information, classified by bond dimension
......@@ -82,6 +89,13 @@ function non_hermitian_evaluate(m::PhysModel, bondD::Integer, β::Real, ResultFo
@unpack (init, solver, trace_estimator,
compress_options, optim_options,
Continue, max_pow_step, group, tag) = options
if trace_estimator !== nothing
new_options = FTLMOptions(trace_estimator.options, processor=processor)
trace_estimator = TraceEstimator(trace_estimator.estimator, new_options)
compress_options = CompressOptions(compress_options, trace_estimator = trace_estimator)
end
solver = solver_function(processor)
"""
ResultFolder: classified by Model parameters(interaction, width), store CMPS and Obsv files
CMPSResultFolder: CMPS information, classified by bond dimension
......
src/logtrexp.jl
View file @ 2f2350ff
......@@ -9,14 +9,19 @@ function logtrexp(t::Real, M::AbstractMatrix)
return logsumexp(t*vals)
end
logtrexp(t::Real, M::CMPSMatrix) = logtrexp(t, Matrix(M))
logtrexp(t, M, esitmator::Nothing) = logtrexp(t, M)
logtrexp(t, M, trace_esitmator::Nothing) = logtrexp(t, M)
function logtrexp(t::Real, M::CMPSMatrix, esitmator::TraceEstimator)
function logtrexp(t::Real, M, trace_estimator::TraceEstimator)
sign(t) == 1 ? which = :LR : which=:SR
e0, _, _ = eigsolve(M, size(M,1), 1, which, ishermitian = true)
e0 = e0[1]
expr = e -> exp(t * (e - e0))
res = FiniteTLanczos.evaluate(esitmator, M, expr)[1]
return log(res) + t*e0
@unpack estimator, options = trace_estimator
@unpack processor = options
processor == CPU ? x0 = rand(size(M,1)) : x0 = CUDA.rand(Float64, size(M,1))
e0, _, _ = eigsolve(M, x0, 1, which, ishermitian = true)
e1 = e0[1]
expr = e -> exp(t * (e - e1))
options = FTLMOptions(options, which = which)
res = FiniteTLanczos.evaluate(TraceEstimator(estimator, options), M, expr)[1]
return log(res) + t*e1
end
src/rrule/logtrexp.jl
View file @ 2f2350ff
Zygote.@adjoint Array(x::CuArray) = Array(x), dy->(CuArray(dy),)
Zygote.@adjoint CMPSMatrix(ψl, ψr) = CMPSMatrix(ψl, ψr), ȳ->(ȳ.ψl, ȳ.ψr)
Zygote.@adjoint CTensor(x::T) where T<:CMPSMatrix = CTensor(x), ȳ->(CTensor(CMPSMatrix(ȳ.ψl, ȳ.ψr)), )
Zygote.@adjoint CuCTensor(x::T) where T<:CMPSMatrix = CuCTensor(x), ȳ->(CuCTensor(CMPSMatrix(ȳ.ψl, ȳ.ψr)),)
include("./logtrexp_Full_ED.jl")
include("./logtrexp_simple_FTLM.jl")
include("./logtrexp_orthogonalized_FTLM.jl")
include("./logtrexp_replaced_FTLM.jl")
include("./logtrexp_FullSampling_FTLM.jl")
\ No newline at end of file
include("./logtrexp_FullSampling_FTLM.jl")
src/rrule/logtrexp_FullSampling_FTLM.jl
View file @ 2f2350ff
"""
rrule for logtrexp(tM) function where `typeof(M)=CMPSMatrix`,
`estimator = FullSampling_FTLM`
`trace_estimator = FullSampling_FTLM`
"""
function ChainRules.rrule(::typeof(logtrexp),
t::Real, M::CMPSMatrix{Ts,T,S,U},
estimator::TraceEstimator{Tf, To}
trace_estimator::TraceEstimator{Tf, To}
) where {Ts,T,S,U,Tf<:typeof(FullSampling_FTLM),To}
@unpack options = estimator
@unpack Nk = options
@unpack options = trace_estimator
@unpack Nk, processor = options
@unpack ψl, ψr = M
Ns = size(M, 1)
Nk = min(Nk, Ns)
χl, χr = size(ψl.Q, 1), size(ψr.Q, 1)
solver = solver_function(processor)
sign(t) == 1 ? which = :LR : which=:SR
e0, _, _ = eigsolve(M, size(M,1), 1, which, ishermitian = true)
e0 = e0[1]
expr_Λ = e -> exp(t*(e-e0))
expr_∂y_∂t = e -> e * exp(t*(e-e0))
processor == CPU ? x0 = rand(T,Ns) : x0 = CUDA.rand(T, Ns)
e0, _, _ = eigsolve(M, x0, 1, which, ishermitian = true)
e1 = e0[1]
expr_Λ = e -> exp(t*(e-e1))
expr_∂y_∂t = e -> e * exp(t*(e-e1))
expr = (expr_Λ, expr_∂y_∂t)
ortho_basis = basis_generate(Ns)
res = zeros(2)
∂y_∂Ql, ∂y_∂Qr, ∂y_∂Rl, ∂y_∂Rr = zeros(4)
if processor == CPU
∂y_∂Ql, ∂y_∂Qr = zeros(T, size(ψl.Q)), zeros(T, size(ψr.Q))
∂y_∂Rl, ∂y_∂Rr = zeros(T, size(ψl.R)), zeros(T, size(ψr.R))
else
∂y_∂Ql, ∂y_∂Qr = CUDA.zeros(T, size(ψl.Q)), CUDA.zeros(T, size(ψr.Q))
∂y_∂Rl, ∂y_∂Rr = CUDA.zeros(T, size(ψl.R)), CUDA.zeros(T, size(ψr.R))
end
for r = 1: Ns
v0 = ortho_basis[:, r]
v0 = solver(x->x, v0)
@unpack init_vector, weight, values, vectors = itFOLM(M, init_vector = v0, Nk = Nk) |> eigensolver
func = f -> map((e,w)->f(e)* w * w', values, weight) |> sum
#Nk = size(values,1)
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))
Λ = map(expr_Λ, values)
vecs = reshape(vectors, χr, χl, Nk)
vectors = reshape(vectors, χr, χl, Nk)
v0 = reshape(init_vector, χr, χl)
Onel = ones(χl, χl); Onel = convert(S, Onel)
Oner = ones(χr, χr); Oner = convert(S, Oner)
∂y_∂Ql += -t * ein"n,n,kbn,ka,kk -> ab"(Λ, weight, conj(vecs), v0, Oner)
∂y_∂Qr += -t * ein"n,n,bkn,ak,kk -> ab"(Λ, weight, conj(vecs), v0, Onel)
if processor == CPU
Onel = ones(T, χl, χl)
Oner = ones(T, χr, χr)
else
Onel = CUDA.ones(T, χl, χl)
Oner = CUDA.ones(T, χr, χr)
end
∂y_∂Ql_temp = -t * ein"n,n,kbn,ka,kk -> ab"(Λ, weight, conj(vectors), v0, Oner)
∂y_∂Qr_temp = -t * ein"n,n,bkn,ak,kk -> ab"(Λ, weight, conj(vectors), v0, Onel)
∂y_∂Ql = map(+, ∂y_∂Ql, ∂y_∂Ql_temp)
∂y_∂Qr = map(+, ∂y_∂Qr, ∂y_∂Qr_temp)
if U <: AbstractMatrix
∂y_∂Rl += -t * ein"n,n,kl,lbn,ka,kl -> ab"(Λ, weight, ψr.R, conj(vecs), v0, Oner)
∂y_∂Rr += -t * ein"n,n,kl,bln,ak,kl -> ab"(Λ, weight, ψl.R, conj(vecs), v0, Onel)
∂y_∂Rl_temp = -t * ein"n,n,kl,lbn,ka,kl -> ab"(Λ, weight, ψr.R, conj(vectors), v0, Oner)
∂y_∂Rr_temp = -t * ein"n,n,kl,bln,ak,kl -> ab"(Λ, weight, ψl.R, conj(vectors), v0, Onel)
else
Onel = ones(χl, χl, size(ψ.R,3)); Onel = convert(U, Onel)
Oner = ones(χr, χr, size(ψ.R,3)); Oner = convert(U, Oner)
∂y_∂Rl += -t * ein"n,n,klm,lbn,ka,klm -> ab"(Λ, weight, ψr.R, conj(vecs), v0, Oner)
∂y_∂Rr += -t * ein"n,n,klm,bln,ak,klm -> ab"(Λ, weight, ψl.R, conj(vecs), v0, Onel)
if processor == CPU
Onel = ones(T, χl, χl, size(ψl.R,3))
Oner = ones(T, χr, χr, size(ψr.R,3))
else
Onel = CUDA.ones(T, χl, χl, size(ψl.R,3))
Oner = CUDA.ones(T, χr, χr, size(ψr.R,3))
end
∂y_∂Rl_temp = -t * ein"n,n,klm,lbn,ka,klm -> abm"(Λ, weight, ψr.R, conj(vectors), v0, Oner)
∂y_∂Rr_temp = -t * ein"n,n,klm,bln,ak,klm -> abm"(Λ, weight, ψl.R, conj(vectors), v0, Onel)
end
∂y_∂Rl = map(+, ∂y_∂Rl, ∂y_∂Rl_temp)
∂y_∂Rr = map(+, ∂y_∂Rr, ∂y_∂Rr_temp)
end
y, ∂y_∂t = map(x -> x * Ns/Nr, res)
∂y_∂Ql, ∂y_∂Qr, ∂y_∂Rl, ∂y_∂Rr = map(x -> x * Ns/Nr, (∂y_∂Ql, ∂y_∂Qr, ∂y_∂Rl, ∂y_∂Rr))
y, ∂y_∂t = res
function logtrexp_pullback(ȳ)
∂y_∂t, ∂y_∂Ql, ∂y_∂Qr, ∂y_∂Rl, ∂y_∂Rr = map(x -> x/y, (∂y_∂t, ∂y_∂Ql, ∂y_∂Qr, ∂y_∂Rl, ∂y_∂Rr))
......
src/rrule/logtrexp_Full_ED.jl
View file @ 2f2350ff
......@@ -64,15 +64,18 @@ end
"""
rrule for logtrexp(tM) function where `typeof(M)=CMPSMatrix`,
`∂y_∂M` also return a `CMPSMatrix`. When `estimator = Full_ED`,
`∂y_∂M` also return a `CMPSMatrix`. When `trace_estimator = Full_ED`,
`exp(tM)` is not explicitly constructed.
(Note there is no memory saving in this case because one has to
save all eigen vectors of the CMPSMatrix.)
"""
function ChainRules.rrule(::typeof(logtrexp),
t::Real, M::CMPSMatrix{Ts,T,S,U},
estimator::typeof(Full_ED)) where {Ts,T,S,U}
trace_estimator::TraceEstimator{Tf, To}
) where {Ts,T,S,U,Tf<:typeof(Full_ED),To}
@unpack ψl, ψr = M
@unpack processor = trace_estimator.options
χl, χr = size(ψl.Q, 1), size(ψr.Q, 1)
vals, vecs = eigensolver(M)
y = logsumexp(t*vals)
......@@ -85,8 +88,13 @@ function ChainRules.rrule(::typeof(logtrexp),
vecs = reshape(vecs, χr, χl, Ns)
χl, χr = size(ψl.Q, 1), size(ψr.Q, 1)
Onel = ones(χl, χl); Onel = convert(S, Onel)
Oner = ones(χr, χr); Oner = convert(S, Oner)
if processor == CPU
Onel = ones(T, χl, χl)
Oner = ones(T, χr, χr)
else
Onel = CUDA.ones(T, χl, χl)
Oner = CUDA.ones(T, χr, χr)
end
∂y_∂Ql = -t * ein"n,kbn,kan,kk -> ab"(Λ, conj(vecs), vecs, Oner)
∂y_∂Qr = -t * ein"n,bkn,akn,kk -> ab"(Λ, conj(vecs), vecs, Onel)
......@@ -94,10 +102,15 @@ function ChainRules.rrule(::typeof(logtrexp),
∂y_∂Rl = -t * ein"n,kl,lbn,kan,kl -> ab"(Λ, ψr.R, conj(vecs), vecs, Oner)
∂y_∂Rr = -t * ein"n,kl,bln,akn,kl -> ab"(Λ, ψl.R, conj(vecs), vecs, Onel)
else
Onel = ones(χl, χl, size(ψ.R,3)); Onel = convert(U, Onel)
Oner = ones(χr, χr, size(ψ.R,3)); Oner = convert(U, Oner)
∂y_∂Rl = -t * ein"n,klm,lbn,kan,kl -> ab"(Λ, ψr.R, conj(vecs), vecs, Oner)
∂y_∂Rr = -t * ein"n,klm,bln,akn,kl -> ab"(Λ, ψl.R, conj(vecs), vecs, Onel)
if processor == CPU
Onel = ones(T, χl, χl, size(ψl.R,3))
Oner = ones(T, χr, χr, size(ψr.R,3))
else
Onel = CUDA.ones(T, χl, χl, size(ψl.R,3))
Oner = CUDA.ones(T, χr, χr, size(ψr.R,3))
end
∂y_∂Rl = -t * ein"n,klm,lbn,kan,klm -> abm"(Λ, ψr.R, conj(vecs), vecs, Oner)
∂y_∂Rr = -t * ein"n,klm,bln,akn,klm -> abm"(Λ, ψl.R, conj(vecs), vecs, Onel)
end
Q̄l = ȳ * ∂y_∂Ql
R̄l = ȳ * ∂y_∂Rl
......
src/rrule/logtrexp_orthogonalized_FTLM.jl
View file @ 2f2350ff
"""
rrule for logtrexp(tM) function where `typeof(M)=CMPSMatrix`,
`estimator = orthogonalized_FTLM`
`trace_estimator = orthogonalized_FTLM`
"""
function ChainRules.rrule(::typeof(logtrexp),
t::Real, M::CMPSMatrix{Ts,T,S,U},
estimator::TraceEstimator{Tf, To}
trace_estimator::TraceEstimator{Tf, To}
) where {Ts,T,S,U,Tf<:typeof(orthogonalized_FTLM),To}
@unpack options = estimator
@unpack distr, Nr, Nk, Ne = options
@unpack options = trace_estimator
@unpack distr, Nr, Nk, Ne, processor = options
@unpack ψl, ψr = M
solver = solver_function(processor)
Ns = size(M, 1)
Ne = min(Ne, Ns-1)
Ne = min(Ne, Ns)
Nk = min(Nk, Ns - Ne)
χl, χr = size(ψl.Q, 1), size(ψr.Q, 1)
sign(t) == 1 ? which = :LR : which=:SR
e0, _, _ = eigsolve(M, size(M,1), 1, which, ishermitian = true)
e0 = e0[1]
expr_Λ = e -> exp(t*(e-e0))
expr_∂y_∂t = e -> e * exp(t*(e-e0))
processor == CPU ? x0 = rand(T, Ns) : x0 = CUDA.rand(T, Ns)
e0, _, _ = eigsolve(M, x0, 1, which, ishermitian = true)
e1 = e0[1]
expr_Λ = e -> exp(t*(e-e1))
expr_∂y_∂t = e -> e * exp(t*(e-e1))
expr = (expr_Λ, expr_∂y_∂t)
#exact eigen pair part
evals, evecs, _ = eigsolve(M, Ns, Ne, :SR, ishermitian = true)
evals = evals[1:Ne]
evecs = hcat(evecs[1:Ne]...)
krylovdim = max(30, Ne+1)
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]...)
func1 = f -> map(f, evals) |> sum
func1 = f -> reduce(+, map(f, eigen_vals))
res1 = map(func1, expr) #store y1, ∂y_∂t1
Λ1 = map(expr_Λ, eigen_vals)
eigen_vecs2 = reshape(eigen_vecs, χr, χl, Ne)
if processor == CPU
Onel = ones(T, χl, χl)
Oner = ones(T, χr, χr)
else
Onel = CUDA.ones(T, χl, χl)
Oner = CUDA.ones(T, χr, χr)
end
Λ = map(expr_Λ, evals)
evecs = reshape(evecs, χr, χl, Ne)
Onel = ones(χl, χl); Onel = convert(S, Onel)
Oner = ones(χr, χr); Oner = convert(S, Oner)
∂y_∂Ql1 = -t * ein"n,kbn,kan,kk -> ab"(Λ, conj(evecs), evecs, Oner)
∂y_∂Qr1 = -t * ein"n,bkn,akn,kk -> ab"(Λ, conj(evecs), evecs, Onel)
∂y_∂Ql1 = -t * ein"n,kbn,kan,kk -> ab"(Λ1, conj(eigen_vecs2), eigen_vecs2, Oner)
∂y_∂Qr1 = -t * ein"n,bkn,akn,kk -> ab"(Λ1, conj(eigen_vecs2), eigen_vecs2, Onel)
if U <: AbstractMatrix
∂y_∂Rl1 = -t * ein"n,kl,lbn,kan,kl -> ab"(Λ, ψr.R, conj(evecs), evecs, Oner)
∂y_∂Rr1 = -t * ein"n,kl,bln,akn,kl -> ab"(Λ, ψl.R, conj(evecs), evecs, Onel)
∂y_∂Rl1 = -t * ein"n,kl,lbn,kan,kl -> ab"(Λ1, ψr.R, conj(eigen_vecs2), eigen_vecs2, Oner)
∂y_∂Rr1 = -t * ein"n,kl,bln,akn,kl -> ab"(Λ1, ψl.R, conj(eigen_vecs2), eigen_vecs2, Onel)
else
Onel = ones(χl, χl, size(ψ.R,3)); Onel = convert(U, Onel)
Oner = ones(χr, χr, size(ψ.R,3)); Oner = convert(U, Oner)
∂y_∂Rl1 = -t * ein"n,klm,lbn,kan,kl -> ab"(Λ, ψr.R, conj(evecs), evecs, Oner)
∂y_∂Rr1 = -t * ein"n,klm,bln,akn,kl -> ab"(Λ, ψl.R, conj(evecs), evecs, Onel)
if processor == CPU
Onel = ones(T, χl, χl, size(ψl.R,3))
Oner = ones(T, χr, χr, size(ψr.R,3))
else
Onel = CUDA.ones(T, χl, χl, size(ψl.R,3))
Oner = CUDA.ones(T, χr, χr, size(ψr.R,3))
end
∂y_∂Rl1 = -t * ein"n,klm,lbn,kan,klm -> abm"(Λ1, ψr.R, conj(eigen_vecs2), eigen_vecs2, Oner)
∂y_∂Rr1 = -t * ein"n,klm,bln,akn,klm -> abm"(Λ1, ψl.R, conj(eigen_vecs2), eigen_vecs2, Onel)
end
#lanczos part
res2 = zeros(2)
∂y_∂Ql2, ∂y_∂Qr2, ∂y_∂Rl2, ∂y_∂Rr2 = zeros(4)
for r = 1: Nr
v0 = random_unit_vector(Ns, distr)
@unpack init_vector, weight, values, vectors =
itFOLM(M, evecs, init_vector = v0, Nk = Nk) |> eigensolver
if processor == CPU
∂y_∂Ql2, ∂y_∂Qr2 = zeros(T, size(ψl.Q)), zeros(T, size(ψr.Q))
∂y_∂Rl2, ∂y_∂Rr2 = zeros(T, size(ψl.R)), zeros(T, size(ψr.R))
else
∂y_∂Ql2, ∂y_∂Qr2 = CUDA.zeros(T, size(ψl.Q)), CUDA.zeros(T, size(ψr.Q))
∂y_∂Rl2, ∂y_∂Rr2 = CUDA.zeros(T, size(ψl.R)), CUDA.zeros(T, size(ψr.R))
end
if Nk > 0
for r = 1: Nr
v0 = random_unit_vector(Ns, distr)
v0 = solver(x->x, v0)
@unpack init_vector, weight, values, vectors =
itFOLM(M, eigen_vecs, init_vector = v0, Nk = Nk) |> eigensolver
#Nk = size(values,1)
func2 = f -> map((e,w)->f(e)* w * w', values, weight) |> sum
res2 = map(+, res2, map(func2, expr))
func2 = f -> begin
Λ = map(f, values)
Z = ein"i,i,i -> "(Λ, weight, conj(weight))
return Array(Z)[1]
end
res2 = map(+, res2, map(func2, expr))
Λ = map(expr_Λ, values)
vecs = reshape(vectors, χr, χl, Nk)
v0 = reshape(init_vector, χr, χl)
Λ2 = map(expr_Λ, values)
vectors = reshape(vectors, χr, χl, Nk)
v0 = reshape(init_vector, χr, χl)
Onel = ones(χl, χl); Onel = convert(S, Onel)
Oner = ones(χr, χr); Oner = convert(S, Oner)
∂y_∂Ql2 += -t * ein"n,n,kbn,ka,kk -> ab"(Λ, weight, conj(vecs), v0, Oner)
∂y_∂Qr2 += -t * ein"n,n,bkn,ak,kk -> ab"(Λ, weight, conj(vecs), v0, Onel)
if processor == CPU
Onel = ones(T, χl, χl)
Oner = ones(T, χr, χr)
else
Onel = CUDA.ones(T, χl, χl)
Oner = CUDA.ones(T, χr, χr)
end
∂y_∂Ql2_temp = -t * ein"n,n,kbn,ka,kk -> ab"(Λ2, weight, conj(vectors), v0, Oner)
∂y_∂Qr2_temp = -t * ein"n,n,bkn,ak,kk -> ab"(Λ2, weight, conj(vectors), v0, Onel)
∂y_∂Ql2 = map(+, ∂y_∂Ql2, ∂y_∂Ql2_temp)
∂y_∂Qr2 = map(+, ∂y_∂Qr2, ∂y_∂Qr2_temp)
if U <: AbstractMatrix
∂y_∂Rl2 += -t * ein"n,n,kl,lbn,ka,kl -> ab"(Λ, weight, ψr.R, conj(vecs), v0, Oner)
∂y_∂Rr2 += -t * ein"n,n,kl,bln,ak,kl -> ab"(Λ, weight, ψl.R, conj(vecs), v0, Onel)
else
Onel = ones(χl, χl, size(ψ.R,3)); Onel = convert(U, Onel)
Oner = ones(χr, χr, size(ψ.R,3)); Oner = convert(U, Oner)
∂y_∂Rl2 += -t * ein"n,n,klm,lbn,ka,klm -> ab"(Λ, weight, ψr.R, conj(vecs), v0, Oner)
∂y_∂Rr2 += -t * ein"n,n,klm,bln,ak,klm -> ab"(Λ, weight, ψl.R, conj(vecs), v0, Onel)
if U <: AbstractMatrix
∂y_∂Rl2_temp = -t * ein"n,n,kl,lbn,ka,kl -> ab"(Λ2, weight, ψr.R, conj(vectors), v0, Oner)
∂y_∂Rr2_temp = -t * ein"n,n,kl,bln,ak,kl -> ab"(Λ2, weight, ψl.R, conj(vectors), v0, Onel)
else
if processor == CPU
Onel = ones(T, χl, χl, size(ψl.R,3))
Oner = ones(T, χr, χr, size(ψr.R,3))
else
Onel = CUDA.ones(T, χl, χl, size(ψl.R,3))
Oner = CUDA.ones(T, χr, χr, size(ψr.R,3))
end
∂y_∂Rl2_temp = -t * ein"n,n,klm,lbn,ka,klm -> abm"(Λ2, weight, ψr.R, conj(vectors), v0, Oner)
∂y_∂Rr2_temp = -t * ein"n,n,klm,bln,ak,klm -> abm"(Λ2, weight, ψl.R, conj(vectors), v0, Onel)
end
∂y_∂Rl2 = map(+, ∂y_∂Rl2, ∂y_∂Rl2_temp)
∂y_∂Rr2 = map(+, ∂y_∂Rr2, ∂y_∂Rr2_temp)
end
end
res2 = map(x -> x * Ns/Nr, res2)
∂y_∂Ql2, ∂y_∂Qr2, ∂y_∂Rl2, ∂y_∂Rr2 = map(x -> x * Ns/Nr, (∂y_∂Ql, ∂y_∂Qr, ∂y_∂Rl, ∂y_∂Rr))
factor = (Ns-Ne)/Nr
res2 = map(x -> x * factor, res2)
∂y_∂Ql2, ∂y_∂Qr2, ∂y_∂Rl2, ∂y_∂Rr2 = map(x -> x * factor, (∂y_∂Ql2, ∂y_∂Qr2, ∂y_∂Rl2, ∂y_∂Rr2))
y, ∂y_∂t = map(+, res1, res2)
∂y_∂Ql, ∂y_∂Qr, ∂y_∂Rl, ∂y_∂Rr = map(+,(∂y_∂Ql1, ∂y_∂Qr1, ∂y_∂Rl1, ∂y_∂Rr1),
......
src/rrule/logtrexp_replaced_FTLM.jl
View file @ 2f2350ff
"""
rrule for logtrexp(tM) function where `typeof(M)=CMPSMatrix`,
`estimator = replaced_FTLM`
`trace_estimator = replaced_FTLM`
"""
function ChainRules.rrule(::typeof(logtrexp),
t::Real, M::CMPSMatrix{Ts,T,S,U},
estimator::TraceEstimator{Tf, To}
trace_estimator::TraceEstimator{Tf, To}
) where {Ts,T,S,U,Tf<:typeof(replaced_FTLM),To}
@unpack options = estimator
@unpack distr, Nr, Nk, Ne = options
@unpack options = trace_estimator
@unpack distr, Nr, Nk, Ne, processor = options
@unpack ψl, ψr = M
solver = solver_function(processor)
Ns = size(M, 1)
Ne = min(Ne, Ns-1)
Ne = min(Ne, Ns)
Nk = min(max(Nk, Ne), Ns)
χl, χr = size(ψl.Q, 1), size(ψr.Q, 1)
sign(t) == 1 ? which = :LR : which=:SR
e0, _, _ = eigsolve(M, size(M,1), 1, which, ishermitian = true)
e0 = e0[1]
expr_Λ = e -> exp(t*(e-e0))
expr_∂y_∂t = e -> e * exp(t*(e-e0))
processor == CPU ? x0 = rand(T,Ns) : x0 = CUDA.rand(T, Ns)
e0, _, _ = eigsolve(M, x0, 1, which, ishermitian = true)
e1 = e0[1]
expr_Λ = e -> exp(t*(e-e1))
expr_∂y_∂t = e -> e * exp(t*(e-e1))
expr = (expr_Λ, expr_∂y_∂t)
evals, evecs, _ = eigsolve(M, Ns, Ne, :SR, ishermitian = true)
evals = evals[1:Ne]
evecs = hcat(evecs[1:Ne]...)
krylovdim = max(30, Ne+1)
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]...)
res = zeros(2)
∂y_∂Ql, ∂y_∂Qr, ∂y_∂Rl, ∂y_∂Rr = zeros(4)
if processor == CPU
∂y_∂Ql, ∂y_∂Qr = zeros(T, size(ψl.Q)), zeros(T, size(ψr.Q))
∂y_∂Rl, ∂y_∂Rr = zeros(T, size(ψl.R)), zeros(T, size(ψr.R))
else
∂y_∂Ql, ∂y_∂Qr = CUDA.zeros(T, size(ψl.Q)), CUDA.zeros(T, size(ψr.Q))
∂y_∂Rl, ∂y_∂Rr = CUDA.zeros(T, size(ψl.R)), CUDA.zeros(T, size(ψr.R))
end
Nk = max(Nk, Ne)
for r = 1: Nr
v0 = random_unit_vector(Ns, distr)
v0 = solver(x->x, v0)
@unpack init_vector, weight, values, vectors =
itFOLM(M, init_vector = v0, Nk = Nk) |> eigensolver
#Nk = size(values,1)
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])
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])
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))
Λ = map(expr_Λ, values)
vecs = reshape(vectors, χr, χl, Nk)
v0 = reshape(init_vector, χr, χl)
Onel = ones(χl, χl); Onel = convert(S, Onel)
Oner = ones(χr, χr); Oner = convert(S, Oner)
∂y_∂Ql += -t * ein"n,n,kbn,ka,kk -> ab"(Λ, weight, conj(vecs), v0, Oner)
∂y_∂Qr += -t * ein"n,n,bkn,ak,kk -> ab"(Λ, weight, conj(vecs), v0, Onel)
if processor == CPU
Onel = ones(T, χl, χl)
Oner = ones(T, χr, χr)
else
Onel = CUDA.ones(T, χl, χl)
Oner = CUDA.ones(T, χr, χr)
end
∂y_∂Ql_temp = -t * ein"n,n,kbn,ka,kk -> ab"(Λ, weight, conj(vecs), v0, Oner)
∂y_∂Qr_temp = -t * ein"n,n,bkn,ak,kk -> ab"(Λ, weight, conj(vecs), v0, Onel)
∂y_∂Ql = map(+, ∂y_∂Ql, ∂y_∂Ql_temp)
∂y_∂Qr = map(+, ∂y_∂Qr, ∂y_∂Qr_temp)
if U <: AbstractMatrix
∂y_∂Rl += -t * ein"n,n,kl,lbn,ka,kl -> ab"(Λ, weight, ψr.R, conj(vecs), v0, Oner)
∂y_∂Rr += -t * ein"n,n,kl,bln,ak,kl -> ab"(Λ, weight, ψl.R, conj(vecs), v0, Onel)
∂y_∂Rl_temp = -t * ein"n,n,kl,lbn,ka,kl -> ab"(Λ, weight, ψr.R, conj(vecs), v0, Oner)
∂y_∂Rr_temp = -t * ein"n,n,kl,bln,ak,kl -> ab"(Λ, weight, ψl.R, conj(vecs), v0, Onel)
else
Onel = ones(χl, χl, size(ψ.R,3)); Onel = convert(U, Onel)
Oner = ones(χr, χr, size(ψ.R,3)); Oner = convert(U, Oner)
∂y_∂Rl += -t * ein"n,n,klm,lbn,ka,klm -> ab"(Λ, weight, ψr.R, conj(vecs), v0, Oner)
∂y_∂Rr += -t * ein"n,n,klm,bln,ak,klm -> ab"(Λ, weight, ψl.R, conj(vecs), v0, Onel)
if processor == CPU
Onel = ones(T, χl, χl, size(ψl.R,3))
Oner = ones(T, χr, χr, size(ψr.R,3))
else
Onel = CUDA.ones(T, χl, χl, size(ψl.R,3))
Oner = CUDA.ones(T, χr, χr, size(ψr.R,3))
end
∂y_∂Rl_temp = -t * ein"n,n,klm,lbn,ka,klm -> abm"(Λ, weight, ψr.R, conj(vecs), v0, Oner)
∂y_∂Rr_temp = -t * ein"n,n,klm,bln,ak,klm -> abm"(Λ, weight, ψl.R, conj(vecs), v0, Onel)
end
∂y_∂Rl = map(+, ∂y_∂Rl, ∂y_∂Rl_temp)
∂y_∂Rr = map(+, ∂y_∂Rr, ∂y_∂Rr_temp)
end
y, ∂y_∂t = map(x -> x * Ns/Nr, res)
∂y_∂Ql, ∂y_∂Qr, ∂y_∂Rl, ∂y_∂Rr = map(x -> x * Ns/Nr, (∂y_∂Ql, ∂y_∂Qr, ∂y_∂Rl, ∂y_∂Rr))
factor = Ns/Nr
y, ∂y_∂t = map(x -> x * factor, res)
∂y_∂Ql, ∂y_∂Qr, ∂y_∂Rl, ∂y_∂Rr = map(x -> x * factor, (∂y_∂Ql, ∂y_∂Qr, ∂y_∂Rl, ∂y_∂Rr))
function logtrexp_pullback(ȳ)
∂y_∂t, ∂y_∂Ql, ∂y_∂Qr, ∂y_∂Rl, ∂y_∂Rr = map(x -> x/y, (∂y_∂t, ∂y_∂Ql, ∂y_∂Qr, ∂y_∂Rl, ∂y_∂Rr))
......
src/rrule/logtrexp_simple_FTLM.jl
View file @ 2f2350ff
"""
rrule for logtrexp(tM) function where `typeof(M)=CMPSMatrix`,
`estimator = simple_FTLM`
`trace_estimator = simple_FTLM`
"""
function ChainRules.rrule(::typeof(logtrexp),
t::Real, M::CMPSMatrix{Ts,T,S,U},
estimator::TraceEstimator{Tf, To}
trace_estimator::TraceEstimator{Tf, To}
) where {Ts,T,S,U,Tf<:typeof(simple_FTLM),To}
@unpack options = estimator
@unpack distr, Nr, Nk = options
@unpack options = trace_estimator
@unpack distr, Nr, Nk, processor = options
@unpack ψl, ψr = M
Ns = size(M, 1)
Nk = min(Nk, Ns)
χl, χr = size(ψl.Q, 1), size(ψr.Q, 1)
solver = solver_function(processor)
sign(t) == 1 ? which = :LR : which=:SR
e0, _, _ = eigsolve(M, size(M,1), 1, which, ishermitian = true)
e0 = e0[1]
expr_Λ = e -> exp(t*(e-e0))
expr_∂y_∂t = e -> e * exp(t*(e-e0))
processor == CPU ? x0 = rand(T,Ns) : x0 = CUDA.rand(T, Ns)
e0, _, _ = eigsolve(M, x0, 1, which, ishermitian = true)
e1 = e0[1]
expr_Λ = e -> exp(t*(e-e1))
expr_∂y_∂t = e -> e * exp(t*(e-e1))
expr = (expr_Λ, expr_∂y_∂t)
res = zeros(2)
∂y_∂Ql, ∂y_∂Qr, ∂y_∂Rl, ∂y_∂Rr = zeros(4)
if processor == CPU
∂y_∂Ql, ∂y_∂Qr = zeros(T, size(ψl.Q)), zeros(T, size(ψr.Q))
∂y_∂Rl, ∂y_∂Rr = zeros(T, size(ψl.R)), zeros(T, size(ψr.R))
else
∂y_∂Ql, ∂y_∂Qr = CUDA.zeros(T, size(ψl.Q)), CUDA.zeros(T, size(ψr.Q))
∂y_∂Rl, ∂y_∂Rr = CUDA.zeros(T, size(ψl.R)), CUDA.zeros(T, size(ψr.R))
end
for r = 1: Nr
v0 = random_unit_vector(Ns, distr)
v0 = solver(x->x, v0)
@unpack init_vector, weight, values, vectors = itFOLM(M, init_vector = v0, Nk = Nk) |> eigensolver
func = f -> map((e,w)->f(e)* w * w', values, weight) |> sum
#Nk = size(values,1)
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))
Λ = map(expr_Λ, values)
vecs = reshape(vectors, χr, χl, Nk)
vectors = reshape(vectors, χr, χl, Nk)
v0 = reshape(init_vector, χr, χl)
Onel = ones(χl, χl); Onel = convert(S, Onel)
Oner = ones(χr, χr); Oner = convert(S, Oner)
∂y_∂Ql += -t * ein"n,n,kbn,ka,kk -> ab"(Λ, weight, conj(vecs), v0, Oner)
∂y_∂Qr += -t * ein"n,n,bkn,ak,kk -> ab"(Λ, weight, conj(vecs), v0, Onel)
if processor == CPU
Onel = ones(T, χl, χl)
Oner = ones(T, χr, χr)
else
Onel = CUDA.ones(T, χl, χl)
Oner = CUDA.ones(T, χr, χr)
end
∂y_∂Ql_temp = -t * ein"n,n,kbn,ka,kk -> ab"(Λ, weight, conj(vectors), v0, Oner)
∂y_∂Qr_temp = -t * ein"n,n,bkn,ak,kk -> ab"(Λ, weight, conj(vectors), v0, Onel)
∂y_∂Ql = map(+, ∂y_∂Ql, ∂y_∂Ql_temp)
∂y_∂Qr = map(+, ∂y_∂Qr, ∂y_∂Qr_temp)
if U <: AbstractMatrix
∂y_∂Rl += -t * ein"n,n,kl,lbn,ka,kl -> ab"(Λ, weight, ψr.R, conj(vecs), v0, Oner)
∂y_∂Rr += -t * ein"n,n,kl,bln,ak,kl -> ab"(Λ, weight, ψl.R, conj(vecs), v0, Onel)
∂y_∂Rl_temp = -t * ein"n,n,kl,lbn,ka,kl -> ab"(Λ, weight, ψr.R, conj(vectors), v0, Oner)
∂y_∂Rr_temp = -t * ein"n,n,kl,bln,ak,kl -> ab"(Λ, weight, ψl.R, conj(vectors), v0, Onel)
else
Onel = ones(χl, χl, size(ψ.R,3)); Onel = convert(U, Onel)
Oner = ones(χr, χr, size(ψ.R,3)); Oner = convert(U, Oner)
∂y_∂Rl += -t * ein"n,n,klm,lbn,ka,klm -> ab"(Λ, weight, ψr.R, conj(vecs), v0, Oner)
∂y_∂Rr += -t * ein"n,n,klm,bln,ak,klm -> ab"(Λ, weight, ψl.R, conj(vecs), v0, Onel)
if processor == CPU
Onel = ones(T, χl, χl, size(ψl.R,3))
Oner = ones(T, χr, χr, size(ψr.R,3))
else
Onel = CUDA.ones(T, χl, χl, size(ψl.R,3))
Oner = CUDA.ones(T, χr, χr, size(ψr.R,3))
end
∂y_∂Rl_temp = -t * ein"n,n,klm,lbn,ka,klm -> abm"(Λ, weight, ψr.R, conj(vectors), v0, Oner)
∂y_∂Rr_temp = -t * ein"n,n,klm,bln,ak,klm -> abm"(Λ, weight, ψl.R, conj(vectors), v0, Onel)
end
∂y_∂Rl = map(+, ∂y_∂Rl, ∂y_∂Rl_temp)
∂y_∂Rr = map(+, ∂y_∂Rr, ∂y_∂Rr_temp)
end
y, ∂y_∂t = map(x -> x * Ns/Nr, res)
∂y_∂Ql, ∂y_∂Qr, ∂y_∂Rl, ∂y_∂Rr = map(x -> x * Ns/Nr, (∂y_∂Ql, ∂y_∂Qr, ∂y_∂Rl, ∂y_∂Rr))
factor = Ns/Nr
y, ∂y_∂t = map(x -> x * factor, res)
∂y_∂Ql, ∂y_∂Qr, ∂y_∂Rl, ∂y_∂Rr = map(x -> x * factor, (∂y_∂Ql, ∂y_∂Qr, ∂y_∂Rl, ∂y_∂Rr))
function logtrexp_pullback(ȳ)
∂y_∂t, ∂y_∂Ql, ∂y_∂Qr, ∂y_∂Rl, ∂y_∂Rr = map(x -> x/y, (∂y_∂t, ∂y_∂Ql, ∂y_∂Qr, ∂y_∂Rl, ∂y_∂Rr))
......
src/solver.jl
View file @ 2f2350ff
function cpu_solver(f::Function, M::AbstractArray...)
M = map(x->convert(Array, x), M)
return f(M...)
end
function gpu_solver(f::Function, M::AbstractArray...)
M = map(x->convert(CuArray, x), M)
return f(M...)
end
TensorUnion = Union{AbstractCTensor, CMPSMatrix}
"""
`cpu_solver` for CTensors and CMPSMatrix
"""
function cpu_solver(f::Function, M::T...) where T<:TensorUnion
#M = Tuple(CTensor(x) for x in M)
M = map(x->CTensor(x), M)
return f(M...)
end
"""
`gpu_solver` for CTensors and CMPSMatrix
"""
function gpu_solver(f::Function, M::T...) where T<:TensorUnion
#M = Tuple(CuCTensor(x) for x in M)
M = map(x->CuCTensor(x), M)
......
src/structs/CMPSCompress.jl
View file @ 2f2350ff
@with_kw struct CompressOptions{Ti<:Union{AbstractCMPS, Nothing},
Tm<:MeraUpdateOptions,
To<:Optim.Options,
Ttrace<:EstimatorType}
Ttrace<:EstimatorType,
Tprocessor<:Processor}
init::Ti = nothing
show_trace::Bool = false
store_trace::Bool = true
trace_estimator::Ttrace = nothing
processor::Tprocessor = CPU
mera_update_options::Tm = MeraUpdateOptions(MeraUpdateOptions(),
show_trace = show_trace,
store_trace = store_trace,
......
src/structs/CMPSMatrix.jl
View file @ 2f2350ff
......@@ -99,4 +99,5 @@ iterate(A::CMPSMatrix, n) = (getindex(A,n), n+1)
convert tensors in CMPSMatrix to CTensor/CuCTensor
"""
CTensor(x::T) where T<:CMPSMatrix = CMPSMatrix(ψl=CTensor(x.ψl), ψr=CTensor(x.ψr))
CuCTensor(x::T) where T<:CMPSMatrix = CMPSMatrix(ψl=CuCTensor(x.ψl), ψr=CuCTensor(x.ψr))
\ No newline at end of file
CuCTensor(x::T) where T<:CMPSMatrix = CMPSMatrix(ψl=CuCTensor(x.ψl), ψr=CuCTensor(x.ψr))
src/structs/EvaluateOptions.jl
View file @ 2f2350ff
@with_kw struct EvaluateOptions{Tinit<:Union{AbstractCMPS, Nothing},
Tsolver<:Function,
Tprocessor<:Processor,
Ttrace<:EstimatorType,
Thermitian<:Union{Bool, Nothing},
Tcontinue<:Union{Bool, Integer},
Tcompress<:CompressOptions,
Toptim<:Optim.Options}
init::Tinit = nothing
solver::Tsolver = cpu_solver
processor::Tprocessor = CPU
trace_estimator::Ttrace = nothing
hermitian::Thermitian = nothing
Continue::Tcontinue = false
......@@ -18,7 +18,8 @@
compress_options::Tcompress = CompressOptions(CompressOptions(),
show_trace = show_trace,
store_trace = store_trace,
trace_estimator = trace_estimator)
trace_estimator = trace_estimator,
processor = processor)
#optim option for hermitian_evaluate
optim_options::Toptim = Optim.Options(f_tol = eps(),
g_tol = 1.e-8,
......
src/utilities/diagm.jl deleted 100644 → 0
View file @ dd9b0a15
#module utilities
"""
diagm(v) function for CuVector
"""
LinearAlgebra.diagm(v::CuVector) = ein"i->ii"(v)
#end # module utilities
test/CMPSMatrix.jl
View file @ 2f2350ff
using LinearAlgebra, Parameters, FiniteTLanczos
using Random; Random.seed!()
Nl, Nr = 10, 15
Nl, Nr = 10, 11
for vD = 1:2
ψl, ψr = init_cmps(Nl, vD), init_cmps(Nr, vD)
......
test/PhysicalModels.jl
View file @ 2f2350ff
using Test, JuliaCMPO
using LinearAlgebra
using LinearAlgebra, PhysModels
@testset "Ising_CMPO" begin
J = 1.0
......
test/cpu_test.jl
View file @ 2f2350ff
using JuliaCMPO, Test
solver = cpu_solver
processor = CPU
solver = solver_function(processor)
@show solver
@testset "PhysicalModels.jl" begin
include("PhysicalModels.jl")
end
@testset "otimes.jl" begin
include("otimes.jl")
......@@ -11,4 +17,12 @@ end
@testset "CMPSMatrix.jl" begin
include("CMPSMatrix.jl")
end
@testset "logtrexp.jl" begin
include("./logtrexp.jl")
end
@testset "gradient/logtrexp.jl" begin
include("./gradient/logtrexp.jl")
end
\ No newline at end of file

Provided by Institute of Theoretical Physics, Chinese Academy of Sciences. Please contact support@itp.ac.cn for more information.