From 84c41e06468ef041bd6097ec64e9967f8024568a Mon Sep 17 00:00:00 2001 From: ysyapa Date: Sun, 19 Mar 2023 14:59:37 -0400 Subject: [PATCH 1/2] Removed testing files in favor of examples --- example.ipynb | 47 ++++++++++++++++++++++++++++++++++ tester.ipynb | 70 --------------------------------------------------- tester.jl | 43 ------------------------------- 3 files changed, 47 insertions(+), 113 deletions(-) create mode 100644 example.ipynb delete mode 100644 tester.ipynb delete mode 100644 tester.jl diff --git a/example.ipynb b/example.ipynb new file mode 100644 index 0000000..dc235e1 --- /dev/null +++ b/example.ipynb @@ -0,0 +1,47 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# prerequisites: KrylovKit, TensorOperations, LinearAlgebra, CUDA#tb/cutensor (for GPU mode)\n", + "\n", + "include(\"CPU.jl\") # using CPU mode\n", + "T = Float32\n", + "\n", + "V_test(r2::T)::T =\n", + " -4 * exp(-r2 / 4)\n", + "\n", + "d = 3\n", + "n = 3\n", + "N = 6\n", + "L::T = 12.0\n", + "mu::T = 0.5\n", + "n_imag = 1\n", + "\n", + "H = HOperator{T}(V_test, 3, 3, N, L, 0.5f0, 1)\n", + "@time evals, evecs, info = eig(H, 5)\n", + "print(info.numops, \" operations : \")\n", + "println(evals)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Julia 1.8.5", + "language": "julia", + "name": "julia-1.8" + }, + "language_info": { + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "1.8.5" + }, + "orig_nbformat": 4 + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/tester.ipynb b/tester.ipynb deleted file mode 100644 index 2878b6c..0000000 --- a/tester.ipynb +++ /dev/null @@ -1,70 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# ./En.run -d 3 -n 3 -e 5 -c eps=0 -c pot=v_gauss,v0=-4,r=2 -N 6 -L 5:14 -c n_imag=1\n", - "\n", - "include(\"CPU.jl\")\n", - "\n", - "T=Float32\n", - "\n", - "function V_test(r2::T)::T\n", - " return -4*exp(-r2/4)\n", - "end\n", - "\n", - "N=6\n", - "for L::T in 5.0:14.0\n", - " H=HOperator{T}(V_test,3,3,N,L,0.5f0,1)\n", - " @time evals,evecs,info=eig(H,5)\n", - " print(info.numops,\" operations : \")\n", - " println(evals)\n", - "end" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# ./En.run -d 3 -n 2 -e 5 -c eps=0 -c pot=v_gauss,v0=-4,r=2 -N 32 -L 5:14 -c n_imag=0\n", - "\n", - "include(\"CPU.jl\")\n", - "\n", - "T=Float64\n", - "\n", - "function V_test(r2::T)::T\n", - " return -4*exp(-r2/4)\n", - "end\n", - "\n", - "N=32\n", - "for L::T in 5.0:14.0\n", - " H=HOperator{T}(V_test,3,2,N,L)\n", - " @time evals,evecs,info=eig(H,5)\n", - " print(info.numops,\" operations : \")\n", - " println(evals)\n", - "end" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Julia 1.8.5", - "language": "julia", - "name": "julia-1.8" - }, - "language_info": { - "file_extension": ".jl", - "mimetype": "application/julia", - "name": "julia", - "version": "1.8.5" - }, - "orig_nbformat": 4 - }, - "nbformat": 4, - "nbformat_minor": 2 -} diff --git a/tester.jl b/tester.jl deleted file mode 100644 index 95f2069..0000000 --- a/tester.jl +++ /dev/null @@ -1,43 +0,0 @@ -# ./En.run -d 3 -n 3 -e 5 -c eps=0 -c pot=v_gauss,v0=-4,r=2 -N 6 -L 5:14 -c n_imag=1 - -using CUDA - -GPU_mode = !("CPU" in ARGS) && CUDA.functional() && CUDA.has_cuda() && CUDA.has_cuda_gpu() - -println("Running with ",Threads.nthreads()," thread(s)") - -if GPU_mode - include("GPU.jl") - println("Available GPUs:") - print(" ") - println.(name.(devices())) -else - include("CPU.jl") -end - -T=Float32 - -function V_zero(r2::T)::T - return 0.0 -end - -function V_test(r2::T)::T - return -4*exp(-r2/4) -end - -N=6 -n_image=1 -μ=0.5 - -levels=5 - -for L::T in 5.0:14.0 - H=HOperator{T}(V_test,3,3,N,L,convert(T,μ),n_image) - if GPU_mode - CUDA.@time evals,evecs,info=eig(H,levels) - else - @time evals,evecs,info=eig(H,levels) - end - print(info.numops," operations : ") - println(evals) -end \ No newline at end of file From abcf18fec26bd025dbfa78b884c69ce58c6b71bb Mon Sep 17 00:00:00 2001 From: ysyapa Date: Sun, 19 Mar 2023 22:28:05 -0400 Subject: [PATCH 2/2] Fixed V algorithm to pick closest pairs --- common.jl | 54 ++++++++++++++++++++++++++++++++---------------------- 1 file changed, 32 insertions(+), 22 deletions(-) diff --git a/common.jl b/common.jl index d187f99..6a63689 100644 --- a/common.jl +++ b/common.jl @@ -1,5 +1,7 @@ Float = Union{Float32,Float64} +norm_square(x::Array{Int})::Int = sum(x .* x) + "Eq (46): Partial derivative matrix element for 1 degree of freedom" function ∂_1DOF(L::T, N::Int, k::Int, l::Int)::Complex{T} where {T<:Float} if k == l @@ -9,35 +11,43 @@ function ∂_1DOF(L::T, N::Int, k::Int, l::Int)::Complex{T} where {T<:Float} end end -"Which index (dimension of the multidimensional array) corresponds to this dimension and coordinate?" -which_index(n::Int, dim::Int, coord::Int)::Int = (dim - 1) * (n - 1) + coord +"Which index (dimension of the multidimensional array) corresponds to spatial dimension 'dim' and particle 'p'?" +which_index(n::Int, dim::Int, p::Int)::Int = (dim - 1) * (n - 1) + p -"k value of the given degree of freedom at the corresponding index, with coord=0 always returning 0" -get_k(n::Int, N::Int, i::CartesianIndex, dim::Int, coord::Int)::Int = - coord == 0 ? 0 : i[which_index(n, dim, coord)] - N ÷ 2 - 1 - -"k value of the DOF at the specified cubic image" -get_shifted_k(n::Int, N::Int, i::CartesianIndex, dim::Int, coord::Int, image::Vector{Int})::Int = - get_k(n, N, i, dim, coord) + N * image[dim] - -"Difference of k values between two particles" -get_Δk(n::Int, N::Int, i::CartesianIndex, dim::Int, coord1::Int, coord2::Int, image::Vector{Int})::Int = - get_k(n, N, i, dim, coord1) - get_shifted_k(n, N, i, dim, coord2, image) +"Δk (distance in terms of lattice paramter) between two particles along the given dimension" +function get_Δk(n::Int, N::Int, i::CartesianIndex, dim::Int, p1::Int, p2::Int)::Int + if p1 == p2 + return 0 + elseif p1 == n + return -(i[which_index(n, dim, p2)] - N ÷ 2 - 1) + elseif p2 == n + return i[which_index(n, dim, p1)] - N ÷ 2 - 1 + else + return i[which_index(n, dim, p1)] - i[which_index(n, dim, p2)] + end +end "Calculate diagonal elements of the V matrix" function calculate_Vs(V_twobody::Function, d::Int, n::Int, N::Int, L::T, n_image::Int)::Array{T} where {T<:Float} L²_over_N² = (L / N)^2 - coords = n - 1 images = collect.(Iterators.product(fill(-n_image:n_image, d)...)) # TODO: Learn how to use tuples instead of vectors - Vs = zeros(T, fill(N, d * coords)...) - for image in images - Threads.@threads for i in CartesianIndices(Vs) - for coord1 in 1:coords - for coord2 in 0:coord1-1 - Δk² = 0 - for dim in 1:d - Δk² += get_Δk(n, N, i, dim, coord1, coord2, image)^2 + Vs = zeros(T, fill(N, d * (n - 1))...) + Threads.@threads for i in CartesianIndices(Vs) + for p1 in 1:n + for p2 in (p1 + 1):n + min_Δk = Array{Int}(undef, d) + for dim in 1:d + Δk = get_Δk(n, N, i, dim, p1, p2) + if Δk > N ÷ 2 + min_Δk[dim] = Δk - N + elseif Δk < -N ÷ 2 + min_Δk[dim] = Δk + N + else + min_Δk[dim] = Δk end + end + for image in images + Δk² = norm_square(min_Δk .- (N .* image)) Vs[i] += V_twobody(Δk² * L²_over_N²) end end