问题描述
我想从矩阵值函数H（Lx，Ly，Lz）获得Array {Float64,3}，其中Lx，Ly，Lz是参数，H是（Lx×Ly×Lz）×（Lx ×Ly×Lz）矩阵。
I want to get Array{Float64,3} from a matrix-valued function H(Lx,Ly,Lz) ,where Lx,Ly,Lz are parameters and H is a (Lx×Ly×Lz)×(Lx×Ly×Lz) matrix.
示例代码为
using LinearAlgebra
eye(T::Type,n) = Diagonal{T}(I, n)
eye(n) = eye(Float64,n)
function H(Lx,Ly,Lz) #def of H
N = Lx*Ly*Lz
mat_Htb = zeros(Complex{Float64},N,N)
for iz = 1:Lz
for ix = 1:Lx
for iy=1:Ly
for dz in -1:1
jz = iz + dz
for dx in -1:1
jx = ix + dx
for dy in -1:1
jy = iy + dy
ii = (iz-1)*Lx*Ly + (ix-1)*Ly + (iy-1) + 1
jj = (jz-1)*Lx*Ly + (jx-1)*Ly + (jy-1) + 1
if 1 <= jx <= Lx && 1 <= jy <= Ly && 1 <= jz <= Lz
if dx == +1 && dy == 0 && dz == 0
mat_Htb[ii,jj] += im
end
if dx == -1 && dy == 0 && dz == 0
mat_Htb[ii,jj] += im/4
end
if dx == 0 && dy == +1 && dz == 0
mat_Htb[ii,jj] += im/2
end
if dx == 0 && dy == -1 && dz == 0
mat_Htb[ii,jj] += im
end
if dx == 0 && dy == 0 && dz == +1
mat_Htb[ii,jj] += -im
end
if dx == 0 && dy == 0 && dz == -1
mat_Htb[ii,jj] += im*(3/7)
end
if dx == 0 && dy == 0 && dz == 0
mat_Htb[ii,jj] += im
end
end
end
end
end
end
end
end
return mat_Htb
end
Lx = 10 #systemsize-parameters
Ly = 10
Lz = 10
ψ0 = Complex{Float64}[] #def of \psi0 ,(Lx×Ly×Lz)×1 vector
for iz = 1:Lz
for ix = 1:Lx
for iy=1:Ly
gauss = exp(-((ix-5)^2 + (iy-5)^2 + (iz-5)^2))
push!(ψ0,gauss)
end
end
end
ψ(t) = exp((-im*t).*H(Lx,Ly,Lz))*ψ0 #time-evolution
abs2ψ(t) = abs2.(ψ(t)./norm(ψ(t))) #normalized density
然后，我尝试制作一个数组{Float64,3}
Then, I tried to make an Array{Float64,3} like this.
x = 1:Lx # our value range
y = 1:Ly
z = 1:Lz
t = 15 #time
ρ(ix,iy,iz) = abs2ψ(t)[(iz-1)*Lx*Ly + (ix-1)*Ly + (iy-1) + 1]
density = Float64[ρ(ix,iy,iz) for ix in x, iy in y,iz in z]
H（Lx，Ly，Lz），ψ（t），abs2ψ（t），ρ（ix， iy，iz）计算顺利。
但是密度大约需要30分钟。
H(Lx,Ly,Lz),ψ(t),abs2ψ(t),ρ(ix,iy,iz) are calculated smoothly.But the density takes about 30 minutes.
最后，我将对t进行循环计算。
所以我想减少计算时间。
您能告诉我如何解决此问题吗？
Ultimately, I will do loop-calculations for t.So I want to reduce the calculation time.Could you tell me how to solve this problem?
推荐答案
仍然有很多事情可能需要进行了改进，但以下版本应该已经比您的版本快得多。
There are still numerous things that would probably need to be improved, but the following version should already be much faster than yours.
要记住的关键是尝试不要多次重复计算同一件事（特别是如果花一些时间来计算）
The key thing to remember is to try and not recompute the same thing several times (especially if it takes some time to compute it, and you're going to re-use the result a large number of times).
在您的示例中，这适用于：
In your example, this applies to :
-
H，仅取决于Lx，Ly和Lz，因此可以一劳永逸地计算 -
ψ和abs2ψ，它们取决于H和t，因此应在每个时间步进行更新-但可以重新用于所有（ix，iy，iz）三元组。
H, which only depends onLx,LyandLz, and as such can be computed once and for allψandabs2ψ, which depend onHandt, and should therefore get updated at each time step -- but can be re-used for all(ix, iy, iz)triplets.
using LinearAlgebra
function H(Lx,Ly,Lz)
N = Lx*Ly*Lz
mat_Htb = zeros(Complex{Float64},N,N)
for iz = 1:Lz
for ix = 1:Lx
for iy=1:Ly
for dz in -1:1
jz = iz + dz
for dx in -1:1
jx = ix + dx
for dy in -1:1
jy = iy + dy
ii = (iz-1)*Lx*Ly + (ix-1)*Ly + (iy-1) + 1
jj = (jz-1)*Lx*Ly + (jx-1)*Ly + (jy-1) + 1
if 1 <= jx <= Lx && 1 <= jy <= Ly && 1 <= jz <= Lz
if dx == +1 && dy == 0 && dz == 0
mat_Htb[ii,jj] += im
end
if dx == -1 && dy == 0 && dz == 0
mat_Htb[ii,jj] += im/4
end
if dx == 0 && dy == +1 && dz == 0
mat_Htb[ii,jj] += im/2
end
if dx == 0 && dy == -1 && dz == 0
mat_Htb[ii,jj] += im
end
if dx == 0 && dy == 0 && dz == +1
mat_Htb[ii,jj] += -im
end
if dx == 0 && dy == 0 && dz == -1
mat_Htb[ii,jj] += im*(3/7)
end
if dx == 0 && dy == 0 && dz == 0
mat_Htb[ii,jj] += im
end
end
end
end
end
end
end
end
return mat_Htb
end
function run(Lx, Ly, Lz)
ψ0 = Complex{Float64}[] #def of \psi0 ,(Lx×Ly×Lz)×1 vector
for iz = 1:Lz
for ix = 1:Lx
for iy=1:Ly
gauss = exp(-((ix-5)^2 + (iy-5)^2 + (iz-5)^2))
push!(ψ0,gauss)
end
end
end
x = 1:Lx # our value range
y = 1:Ly
z = 1:Lz
t = 15 #time
H_ = H(Lx,Ly,Lz)
ψ = exp((-im*t).*H_)*ψ0 #time-evolution
abs2ψ = abs2.(ψ./norm(ψ)) #normalized density
ρ(ix,iy,iz) = abs2ψ[(iz-1)*Lx*Ly + (ix-1)*Ly + (iy-1) + 1]
density = Float64[ρ(ix,iy,iz) for ix in x, iy in y,iz in z]
end
Lx = 10 #systemsize-parameters
Ly = 10
Lz = 10
run(Lx, Ly, Lz)
对于像这样的问题，这些问题非常特定于您要优化的代码，我倾向于认为在会更合适。
这篇关于Julia-在短时间内制作数组{Float64,3}的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！