Following your advice, I tried @code_warntype ground(), but I didn't see any concrete types, I don't know whether there is any other package you use as default.
在 2015年12月15日星期二 UTC+8上午1:45:41，Stefan Karpinski写道： > > 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] <javascript:>> 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 >>> > >>> >>> >
