yapanuwan
/
DVR-jl
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Rename HOperator to Hamiltonian

This commit is contained in:
Nuwan Yapa 2023-04-07 20:02:34 -04:00
parent 03605c060a
commit a8a7bdb44f
3 changed files with 16 additions and 16 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

22
HOperator.jl → Hamiltonian.jl
View File

@ -1,10 +1,10 @@
include("common.jl")
using TensorOperations, KrylovKit, LinearAlgebra, CUDA, CUDA.CUTENSOR
@enum HOperator_backend cpu_tensor gpu_cutensor
@enum Hamiltonian_backend cpu_tensor gpu_cutensor
"A Hamiltonian that can be applied to a vector"
struct HOperator{T}
struct Hamiltonian{T}
d::Int
n::Int
N::Int
@ -15,8 +15,8 @@ struct HOperator{T}
K_mixed # CuTensor{Complex{T}} or Nothing
Vs # Array{Complex{T}} or CuArray{Complex{T}}
hermitian::Bool
mode::HOperator_backend
function HOperator{T}(V_twobody::Function, d::Int, n::Int, N::Int, L::T, ϕ::T, μ::T, n_image::Int, mode::HOperator_backend) where {T<:Float}
mode::Hamiltonian_backend
function Hamiltonian{T}(V_twobody::Function, d::Int, n::Int, N::Int, L::T, ϕ::T, μ::T, n_image::Int, mode::Hamiltonian_backend) where {T<:Float}
@assert mode != gpu_cutensor || CUDA.functional() && CUDA.has_cuda() && CUDA.has_cuda_gpu() "CUDA not available"
k = -N÷2:N÷2-1
Vs = calculate_Vs(V_twobody, d, n, N, L, ϕ, n_image)
@ -33,14 +33,14 @@ struct HOperator{T}
end
end
Base.size(H::HOperator, i::Int)::Int = (i == 1 || i == 2) ? H.N^(H.d * (H.n - 1)) : throw(ArgumentError("HOperator only has 2 dimesions"))
Base.size(H::HOperator)::Dims{2} = (size(H, 1), size(H, 2))
Base.size(H::Hamiltonian, i::Int)::Int = (i == 1 || i == 2) ? H.N^(H.d * (H.n - 1)) : throw(ArgumentError("Hamiltonian only has 2 dimesions"))
Base.size(H::Hamiltonian)::Dims{2} = (size(H, 1), size(H, 2))
"Dimensions of a vector to which 'H' can be applied"
vectorDims(H::HOperator)::Dims = tuple(fill(H.N, H.d * (H.n - 1))...)
vectorDims(H::Hamiltonian)::Dims = tuple(fill(H.N, H.d * (H.n - 1))...)
"Apply 'H' on 'v' and store the result in 'out' using the 'cpu_tensor' backend"
function LinearAlgebra.mul!(out::Array{Complex{T}}, H::HOperator{T}, v::Array{Complex{T}})::Array{Complex{T}} where {T<:Float}
function LinearAlgebra.mul!(out::Array{Complex{T}}, H::Hamiltonian{T}, v::Array{Complex{T}})::Array{Complex{T}} where {T<:Float}
#LinearMaps.check_dim_mul(out,H,v) --- dimensions don't match
# apply V operator
@. out = H.Vs * v
@ -78,7 +78,7 @@ function contract_accumulate!(C::CuTensor, A::CuTensor, B::CuTensor)::CuTensor
end
"Apply 'H' on 'v' and store the result in 'out' using the 'gpu_cutensor' backend"
function LinearAlgebra.mul!(out::CuArray{Complex{T}}, H::HOperator{T}, v::CuArray{Complex{T}})::CuArray{Complex{T}} where {T<:Float}
function LinearAlgebra.mul!(out::CuArray{Complex{T}}, H::Hamiltonian{T}, v::CuArray{Complex{T}})::CuArray{Complex{T}} where {T<:Float}
#LinearMaps.check_dim_mul(out,H,v) --- dimensions don't match
ctx = context()
# apply V operator
@ -123,7 +123,7 @@ function LinearAlgebra.mul!(out::CuArray{Complex{T}}, H::HOperator{T}, v::CuArra
end
"Apply 'H' on 'v' and return the result"
function (H::HOperator)(v)
function (H::Hamiltonian)(v)
out = similar(v)
return mul!(out, H, v)
end
@ -131,7 +131,7 @@ end
tolerance = 1e-6
"Wrapper for KrylovKit.eigsolve"
function eig(H::HOperator{T}, levels::Int; resonances = !H.hermitian)::Tuple{Vector{Complex{T}},Any,Any} where {T<:Float}
function eig(H::Hamiltonian{T}, levels::Int; resonances = !H.hermitian)::Tuple{Vector{Complex{T}},Any,Any} where {T<:Float}
if H.mode == cpu_tensor
x₀ = rand(Complex{T}, vectorDims(H)...)
elseif H.mode == gpu_cutensor

4
benchmark.jl
View File

@ -1,4 +1,4 @@
include("HOperator.jl")
include("Hamiltonian.jl")
GPU_mode = !("CPU" in ARGS) && CUDA.functional() && CUDA.has_cuda() && CUDA.has_cuda_gpu()
@ -31,7 +31,7 @@ n_image=1
for L::T in 5.0:14.0
println("Constructing H operator...")
@time H=HOperator{T}(V_test,3,3,N,L,convert(T,0),convert(T,μ),n_image,mode)
@time H=Hamiltonian{T}(V_test,3,3,N,L,convert(T,0),convert(T,μ),n_image,mode)
println("Applying H 1000 times...")
if GPU_mode
v=CUDA.rand(Complex{T},vectorDims(H)...)

6
example.ipynb
View File

@ -7,7 +7,7 @@
"outputs": [],
"source": [
"# prerequisite packages: KrylovKit, TensorOperations, LinearAlgebra, CUDA#tb/cutensor, Plots\n",
"include(\"HOperator.jl\")\n",
"include(\"Hamiltonian.jl\")\n",
"mode = cpu_tensor # using CPU mode\n",
"T = Float32 # single-precision mode"
]
@ -29,7 +29,7 @@
"μ::T = 0.5\n",
"n_imag = 1\n",
"\n",
"H = HOperator{T}(V_gauss, d, n, N, L, ϕ, μ, n_imag, mode)\n",
"H = Hamiltonian{T}(V_gauss, d, n, N, L, ϕ, μ, n_imag, mode)\n",
"@time evals, evecs, info = eig(H, 5)\n",
"print(info.numops, \" operations : \")\n",
"println(evals)"
@ -54,7 +54,7 @@
"μ::T = 0.5\n",
"n_imag = 0\n",
"\n",
"H = HOperator{T}(V_gauss, d, n, N, L, ϕ, μ, n_imag, mode)\n",
"H = Hamiltonian{T}(V_gauss, d, n, N, L, ϕ, μ, n_imag, mode)\n",
"@time evals, evecs, info = eig(H, 20)\n",
"print(info.numops, \" operations : \")\n",
"print(evals)\n",