Support in-the-wild Julia example ported from Matlab

[lkuper]

A couple weeks ago, someone on julia-users complained about their code running more slowly when porting from Matlab to Julia. Their code uses some Unicode characters which can look weird in some editors, so here's an ASCIIfied version:

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]
    phi = Array(Complex128,(Ks+1,Ks+1))
    Pl = plan_fft!(phi; flags=FFTW.MEASURE)
    phi = V.*complex(1.)

    normphi = sqrt(sumabs2(phi))*dx
    phi /= normphi

    phi = Pl*phi
    phi .*= ko
    phi = Pl\phi

    deltaphi = 1.
    nstep = 0
    print("start loop")
     while deltaphi > 1.e-15
        phi′ = phi
        phi .*= vo

        normphi = sqrt(sumabs2(phi))*dx
        phi /= normphi

        phi = Pl*phi
        phi .*= ko2
        phi = Pl\phi
# if check  deltaphi for every step, 35s is needed.
        if nstep>500
            deltaphi = maxabs(phi-phi′)
        end
        nstep += 1
        if mod(nstep,200)==0
            print(nstep,"  ")
            print("deltaphi=",deltaphi,"   ")
        end

    end
    phi .*= vo

    phi = Pl*phi
    phi .*= ko
    phi = Pl\phi

    normphi = sqrt(sumabs2(phi))*dx
    phi /= normphi
end

ground()

I think it's instructive to try to figure out what it would take for us to be able to handle this. Right now we are missing support for:

• exponentiation
• the updating multiplication operator, *=

There might be more stuff I'm not aware of.

[ninegua]

things like *= are already supported. I just added .*= and a few others.

There is a ticket to support left division A \ B. See #33

There are a few things out side the loop that CGen cannot support, such as the use of keyword argument in plan_fft!(phi; flags=FFTW.MEASURE). But usually these are not considered as part of kernel computation.

[lkuper]

Oops, yeah, I meant .*=, not just *=.

What's the story with exponentiation? How hard would it be to support?

To me, keyword arguments seem less essential.