The code is well-typed – try @code_warntype ground() and you can see that all the local variables have concrete types. The inner loop allocates a lot of new arrays where it could update the same array in place, so you could try avoiding allocation in the inner loop.
On Mon, Dec 14, 2015 at 8:23 AM, Patrick Kofod Mogensen < [email protected]> wrote: > I get it, thanks for pointing that out. > > > On Monday, December 14, 2015 at 11:38:14 AM UTC+1, FQ wrote: >> >> The code does exit at some point, returning a 2049x2049 complex matrix. >> note that delta phi is reassigned inside the loop at some point. >> >> Am 14.12.2015 um 10:31 schrieb Patrick Kofod Mogensen: >> > Are you sure this version works? I get Δϕ == 1.0 in all iterations, so >> I >> > can't see how it would ever exit the while loop... >> > >> > On Monday, December 14, 2015 at 4:35:43 AM UTC+1, 博陈 wrote: >> > >> > This is my first julia code, I am happy it did the right thing, but >> > compared with the Matlab code that did the same thing, it runs so >> > slowly. The Matlab code takes about 90s, but the julia below takes >> 130s. >> > >> > Potential(x::Float64, y::Float64) = -2./sqrt(x^2+1.) >> > -2./sqrt(y^2+1.) + 1./sqrt((x-y)^2+1.) >> > Ko(x::Float64,y::Float64) = x^2+y^2 >> > >> > function ground() >> > R = 320. >> > Ks = 2^11 >> > x = linspace(-R, R, Ks+1) >> > dx = 2R/Ks >> > dt = 0.1 >> > pm = 2π/2dx >> > px = circshift( linspace(-pm, pm, Ks+1), round(Int64, Ks/2) ) >> > FFTW.set_num_threads(8) >> > >> > V = Float64[Potential(ix, iy) for ix in x, iy in x] >> > ko = Float64[e^(-dt/4* (Ko(ipx,ipy))) for ipx in px, ipy in px] >> > ko2 = ko.^2 >> > vo = Float64[e^(-dt* (Potential(ix,iy))) for ix in x, iy in x] >> > ϕ = Array(Complex128,(Ks+1,Ks+1)) >> > Pl = plan_fft!(ϕ; flags=FFTW.MEASURE) >> > ϕ = V.*complex(1.) >> > >> > normphi = sqrt(sumabs2(ϕ))*dx >> > ϕ /= normphi >> > >> > ϕ = Pl*ϕ >> > ϕ .*= ko >> > ϕ = Pl\ϕ >> > >> > Δϕ = 1. >> > nstep = 0 >> > print("start loop") >> > while Δϕ > 1.e-15 >> > ϕ′ = ϕ >> > ϕ .*= vo >> > >> > normphi = sqrt(sumabs2(ϕ))*dx >> > ϕ /= normphi >> > >> > ϕ = Pl*ϕ >> > ϕ .*= ko2 >> > ϕ = Pl\ϕ >> > # if check Δϕ for every step, 35s is needed. >> > if nstep>500 >> > Δϕ = maxabs(ϕ-ϕ′) >> > end >> > nstep += 1 >> > if mod(nstep,200)==0 >> > print(nstep," ") >> > print("Δϕ=",Δϕ," ") >> > end >> > >> > end >> > ϕ .*= vo >> > >> > ϕ = Pl*ϕ >> > ϕ .*= ko >> > ϕ = Pl\ϕ >> > >> > normphi = sqrt(sumabs2(ϕ))*dx >> > ϕ /= normphi >> > end >> > >> >>
