Have you tried: (a) calling @code_typewarn on your function (b) using the built-in profiler?
On Tue, Mar 29, 2016 at 9:23 AM, 博陈 <[email protected]> wrote: > First of all, have a look at the result. > > > <https://lh3.googleusercontent.com/-anNt-E4P1vM/Vvp-TybegZI/AAAAAAAAABE/ZvDO2xarndMSgKVOXy_hcPd5NTh-7QcEA/s1600/QQ%25E5%259B%25BE%25E7%2589%258720160329210732.png> > > > > > > > > > My code calculates the evolution of 1-d 2-electron system in the electric > field, some variables are calculated during the evolution. > According to the result of @time evolution, my code must have a > pre-allocation problem. Before you see the long code, i suggest that the > hotspot might be around the Arrays prop_e, \phio, pp. I have learnt that I > can use m = Array(Float64, 1) outside a "for" loop and empty!(m) and > push!(m, new_m) inside the loop to pre-allocate the variable m, but in my > situations, I don't know how to pre-allocate these arrays. > > Below is the script (precisely, the main function) itself. > > function evolution(ϕ::Array{Complex{Float64}, 2}, > ele::Array{Float64, 1}, dx::Float64, dt::Float64, > flags::Tuple{Int64, Int64, Int64, Int64}) > ϕg = copy(ϕ) > FFTW.set_num_threads(8) > ns = length( ϕ[:, 1] ) > x = get_x(ns, dx) > p = get_p(ns, dx) > if flags[4] == 1 > pp = similar(p) > A = -cumsum(ele) * dt > A² = A.*A > ##### splitting > r_sp = 150.0 > δ_sp = 5.0 > splitter = Array(Float64, ns, ns) > end > nt = length( ele ) > > # ##### Pre-allocate result and temporary arrays > #if flags[1] == 1 > σ = zeros(Complex128, nt) > #end > #if flags[2] == 1 > a = zeros(Float64, nt) > #end > #if flags[3] == 1 > r_ionization = 20.0 > n1 = round(Int, ns/2 - r_ionization/dx) > n2 = round(Int, ns/2 + r_ionization/dx) > ip = zeros(Float64, nt) > #end > > ##### FFT plan > p_fft! = plan_fft!( similar(ϕ), flags=FFTW.MEASURE ) > > prop_x = similar( ϕ ) > prop_p = similar( prop_x ) > prop_e = similar( prop_x ) > # this two versions just cost the same time > xplusy = Array(Float64, ns, ns) > #xplusy = Array( Float64, ns^2) > > ##### absorb boundary > r_a = ns * dx /2 - 50.0 > δ = 10.0 > absorb = Array(Float64, ns, ns) > > k0 = 2π / (ns * dx) > > @inbounds for j in 1:ns > @inbounds for i in 1:ns > prop_x[i, j] = exp( -im * get_potential(x[i], x[j]) * dt / 2 ) > prop_p[i, j] = exp( -im * (p[i]^2 + p[j]^2)/2 * dt ) > > xplusy[i, j] = x[i] + x[j] > > absorb[i, j] = (1.0 - get_out(x[i], r_a, δ ))* (1.0 - > get_out(x[j], > r_a, δ)) > end > end > > if flags[2] == 1 > pvpx = Array(Float64, ns, ns) > @inbounds for j in 1:ns > @inbounds for i in 1:ns > pvpx[i, j] = get_pvpx(x[i], x[j]) > end > end > end > > if flags[4] == 1 > ϕo = zeros(Complex128, ns, ns) > ϕp = zeros(Complex128, ns, ns) > @inbounds for j in 1:ns > @inbounds for i in 1:ns > splitter[i, j] = get_out(x[i], r_sp, δ_sp) * get_out(x[j], > r_sp, δ_sp) > end > end > end > > for i in 1:nt > for j in eachindex(ϕ) > prop_e[j] = exp( -im * ele[i] * xplusy[j] * dt/2.0) > end > > for j in eachindex(ϕ) > ϕ[j] *= prop_x[j] * prop_e[j] > end > p_fft! * ϕ > for j in eachindex(ϕ) > ϕ[j] *= prop_p[j] > end > p_fft! \ ϕ > for j in eachindex(ϕ) > ϕ[j] *= prop_x[j] * prop_e[j] > end > ########## autocorrelation function σ(t) > if flags[1] == 1 > for j in eachindex(ϕ) > σ[i] += conj(ϕg[j]) * ϕ[j] > end > end > ########## dipole acceleration a(t) > if flags[2] == 1 > for j in eachindex(ϕ) > a[i] += abs(ϕ[j])^2 * (pvpx[j] + 2ele[i]) > end > end > ########## ionization probability ip(t) > if flags[3] == 1 > for j1 in n1:n2 > for j2 in 1:ns > ip[i] += abs( ϕ[j2+ns*(j1-1)] )^2 > end > end > for j1 in [1:n1-1; n2+1:ns] > for j2 in n1:n2 > ip[i] += abs( ϕ[j2+ns*(j1-1)] )^2 > end > end > end > ########## get momentum > if flags[4] == 1 > for j in eachindex(ϕo) > ϕo[j] = ϕ[j] * splitter[j] * exp( -im * A[i]*xplusy[j] ) > end > for j in eachindex(p) > pp[j] = p[j]^2 /2 * (nt-i) - p[j] *sum( A[i:nt] ) + sum( > A²[1:nt] ) /2 > end > for j2 in 1:ns > for j1 in 1:ns > ϕo[j1, j2] = ϕo[j1, j2] * exp( -im * (pp[j1] + pp[j2]) > * dt) > end > end > p_fft! * ϕo > for j in eachindex(ϕp) > ϕp[j] += ϕo[j] > end > end > > ########## absorb boundary > if mod(i, 300) == 0 > for j in eachindex(ϕ) > ϕ[j] *= absorb[j] > end > end > > if (mod(i, 500) == 0) > println("i = $i") > flush(STDOUT) > end > end > σ *= dx^2 > a *= dx^2 > ip *= dx^2 > > save("data/fs.jld", "ϕ", ϕ) > if flags[1] == 1 > save("data/sigma.jld", "σ", σ) > end > if flags[2] == 1 > save("data/a.jld", "a", a) > end > if flags[3] == 1 > save("data/ip.jld", "ip", ip) > end > if flags[4] == 1 > save("data/pf.jld", "ϕp", ϕp) > end > > #return σ, a, ip, ϕ > nothing > end > > >
