# Saving values from GPU during Euler stepping in ODE

**URL:** <https://discourse.julialang.org/t/saving-values-from-gpu-during-euler-stepping-in-ode/7406>\
**Category:** Performance\
**Tags:** arrayfire, gpu, gpuarrays\
**Created:** [November 30, 2017, 7:58am UTC](https://discourse.julialang.org/t/saving-values-from-gpu-during-euler-stepping-in-ode/7406 "2017-11-30T07:58:21Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [November 30, 2017, 7:58am UTC](https://discourse.julialang.org/t/saving-values-from-gpu-during-euler-stepping-in-ode/7406/1 "2017-11-30T07:58:21Z")

</div>

Dear all,

I have written a simple Euler stepping algorithm for a particular problem which runs on GPU. The vector `y_current` is big (~200\_000) and it is updated like `y_current .= y_current .+ dt * RHS`. A working example can be found [here](https://github.com/gaika/ArrayFire.jl/issues/35). It was done in ArrayFire. I have a working example in CLArrays but it is 10x slower.

I can’t save every step in a matrix because I need a lot of steps. Instead, I want to record `sum(y_current)` at every step. However, including a save mechanism slows down the whole thing 50 times and it is not useful then.

Hence, I am doing something like

```julia
for ii=1:10000
     updateEuler!(y_current,dt)
     #save data, a is own by the GPU
     a[ii] = sum(y_current)
end

```

Not being written in an array fashion, this is slow indeed.

Does anyone ones a trick to save the sum value by staying on the GPU? May be one has to write a specific kernel…

Thank you for your help and suggestions,

Best regards,

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [November 30, 2017, 2:12pm UTC](https://discourse.julialang.org/t/saving-values-from-gpu-during-euler-stepping-in-ode/7406/2 "2017-11-30T14:12:47Z")

</div>

> [@rveltz](#):
>
> Does anyone ones a trick to save the sum value by staying on the GPU? May be one has to write a specific kernel…

I don’t think that’s the issue here. Returning a scalar from a GPU reduction should be cheap. Why is `a` on the GPU you though?

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [November 30, 2017, 4:29pm UTC](https://discourse.julialang.org/t/saving-values-from-gpu-during-euler-stepping-in-ode/7406/3 "2017-11-30T16:29:51Z")

</div>

The Euler step is something like

`y_current .= y_current .+ dt .* (F.(y_current) .+ sum(y_current) )`

Removing `.+ sum(y_current)` runs much faster going from 1.1s to 0.16s. It don’t expect computing a sum to be so slow.

> Why is a on the GPU you though?

I thought `sum(y_current)` would stay on the gpu, so I initialised `a` on the gpu as well.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [November 30, 2017, 4:35pm UTC](https://discourse.julialang.org/t/saving-values-from-gpu-during-euler-stepping-in-ode/7406/4 "2017-11-30T16:35:13Z")

</div>

> [@rveltz](#):
>
> I thought sum(y\_current) would stay on the gpu, so I initialised a on the gpu as well.

No, it’ll reduce to a scalar on the CPU. That’s why the indexing is slow.

> [@rveltz](#):
>
> y\_current .= y\_current .+ dt .\* (F.(y\_current) .+ sum(y\_current) )

I believe you don’t want to broadcast ArrayFire commands and let its internal memory management handle it. I think if you do that, the scalar add will be much faster?

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [November 30, 2017, 4:56pm UTC](https://discourse.julialang.org/t/saving-values-from-gpu-during-euler-stepping-in-ode/7406/5 "2017-11-30T16:56:24Z")

</div>

> Blockquote I think if you do that, the scalar add will be much faster?

Indeed, it improves a bit.

Now removing the sum changes the timing from 0.85s to 0.13s.

Still, I did not expect the sum to decrease the speed like that.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [November 30, 2017, 5:00pm UTC](https://discourse.julialang.org/t/saving-values-from-gpu-during-euler-stepping-in-ode/7406/6 "2017-11-30T17:00:10Z")

</div>

> [@rveltz](#):
>
> Still, I did not expect the sum to decrease the speed like that.

Do `allowslow(AFArray, false)` and see if your code still runs.

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [November 30, 2017, 5:15pm UTC](https://discourse.julialang.org/t/saving-values-from-gpu-during-euler-stepping-in-ode/7406/7 "2017-11-30T17:15:10Z")

</div>

I use @anon43133060 [version](https://github.com/gaika/ArrayFire.jl/issues) which does not seem to have that functionality. Thank you for the suggestion.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [November 30, 2017, 5:17pm UTC](https://discourse.julialang.org/t/saving-values-from-gpu-during-euler-stepping-in-ode/7406/8 "2017-11-30T17:17:03Z")

</div>

> [@rveltz](#):
>
> I use @anon43133060 version which does not seem to have that functionality. Thank you for the suggestion.

But @anon43133060 implemented that function? Try `ArrayFire.allowslow(AFArray, false)`

---

<div class="post-metadata">

**Author:** ![anon43133060](https://avatars.discourse-cdn.com/v4/letter/a/57b2e6/32.png) [@anon43133060](https://discourse.julialang.org/u/anon43133060)\
**Post date:** [November 30, 2017, 5:27pm UTC](https://discourse.julialang.org/t/saving-values-from-gpu-during-euler-stepping-in-ode/7406/9 "2017-11-30T17:27:06Z")

</div>

The code was merged with official version, my branch is stale now.

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [November 30, 2017, 6:19pm UTC](https://discourse.julialang.org/t/saving-values-from-gpu-during-euler-stepping-in-ode/7406/10 "2017-11-30T18:19:23Z")

</div>

@anon43133060 You mean I can use the official ArrayFire to get your code?

---

<div class="post-metadata">

**Author:** ![anon43133060](https://avatars.discourse-cdn.com/v4/letter/a/57b2e6/32.png) [@anon43133060](https://discourse.julialang.org/u/anon43133060)\
**Post date:** [November 30, 2017, 6:21pm UTC](https://discourse.julialang.org/t/saving-values-from-gpu-during-euler-stepping-in-ode/7406/11 "2017-11-30T18:21:50Z")

</div>

Main reason the `sum` is slow is that it stops JIT and forces evaluation of all `AFArrays`, even temporary.

---

<div class="post-metadata">

**Author:** ![anon43133060](https://avatars.discourse-cdn.com/v4/letter/a/57b2e6/32.png) [@anon43133060](https://discourse.julialang.org/u/anon43133060)\
**Post date:** [November 30, 2017, 6:23pm UTC](https://discourse.julialang.org/t/saving-values-from-gpu-during-euler-stepping-in-ode/7406/12 "2017-11-30T18:23:27Z")

</div>

> [@rveltz](#):
>
> You mean I can use the official ArrayFire to get your code?

Yes, updated description of my fork to that effect.

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [November 30, 2017, 6:26pm UTC](https://discourse.julialang.org/t/saving-values-from-gpu-during-euler-stepping-in-ode/7406/13 "2017-11-30T18:26:56Z")

</div>

Oh, can we “avoid” this behaviour?

The dot product behaves the same?

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [November 30, 2017, 6:34pm UTC](https://discourse.julialang.org/t/saving-values-from-gpu-during-euler-stepping-in-ode/7406/14 "2017-11-30T18:34:40Z")

</div>

Then using u\*y\_current with u=ones(AFArray,1,N) is faster than sum (I tried)… 😂

---

<div class="post-metadata">

**Author:** ![sdanisch](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sdanisch/32/1406_2.png) [@sdanisch](https://discourse.julialang.org/u/sdanisch)\
**Post date:** [December 1, 2017, 11:59am UTC](https://discourse.julialang.org/t/saving-values-from-gpu-during-euler-stepping-in-ode/7406/15 "2017-12-01T11:59:46Z")

</div>

Can you post the full CLArray version? I can see if I can make it faster 🙂

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [December 2, 2017, 12:59pm UTC](https://discourse.julialang.org/t/saving-values-from-gpu-during-euler-stepping-in-ode/7406/16 "2017-12-02T12:59:10Z")

</div>

I can already tell you that I made a mistake. It is not that slow compared to ArrayFire…

```julia
using CLArrays, GPUArrays
CLArrays.init(CLArrays.devices()[2])
const TY = Float32
# using ProgressMeter, StatsBase, Distributions, JLD

const N = 200_000
const w = TY(1.0)
const H0 = TY(2.)
const Iext = TY(1.0)

F(v::TY) = Iext - v*0.1f0
f(v::TY) = max(v-0.2f0,0f0).^20

# F(v) = Iext - v*0.1
# f(v) = max(v-0.2,0).^20

function clstep!(y_current,rates,kappa,dt,a,ii,rnd_cl)
	kappa .= sqrt.(2 * H0 * sum(rates,1)/N) .* w
	y_current .= y_current .+ dt .* (F.(y_current) .+ kappa)
	rates .= f.(y_current)
	GPUArrays.rand!(rnd_cl)
	y_current .= y_current .* (rates .< -log.(rnd_cl) ./dt )
	if mod(ii,500) == 0
		push!(a,kappa)
end

function euler_cl(T::TY,dt::TY,n::Int32,y0, plot_ = false)
    println("\n--> Euler algo, up to T = ",n*dt)
    n = min(n,Int(ceil(T/dt)))
	println("OK")	
    y_current = copy(y0)
    rndcl = copy(y0) # to hold randomness
    rates = f.(y_current)
	println("OK")
    kappa = sqrt.(2 * H0 * sum(rates,1)/N) .* w
    nb_jumps = 0
	y_current .= y_current .+ dt .* (F.(y_current) .+ kappa)
	println("OK for Euler step")	

    # variables for saving data
    tout = TY[]
    a = TY[]
	
    for ii=1:n
        # @show ii
        clstep!(y_current,rates,kappa,dt,a,ii,rndcl)
        # save data
        push!(tout,ii*dt)

    end
    return tout,a
end

#######################################################################################
xc0 = rand(TY,N) |> CLArray;
@show typeof(xc0)

parms = TY.([0., 0.])
dt = TY(1./(5.4N))
t,a_tot = euler_cl(TY(5.5),dt*2,Int32(10),xc0);
t,a_tot = @time euler_cl(TY(5.5),dt,Int32(1000),xc0,false);

```

---

<div class="post-metadata">

**Author:** ![sdanisch](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sdanisch/32/1406_2.png) [@sdanisch](https://discourse.julialang.org/u/sdanisch)\
**Post date:** [December 2, 2017, 2:08pm UTC](https://discourse.julialang.org/t/saving-values-from-gpu-during-euler-stepping-in-ode/7406/17 "2017-12-02T14:08:46Z")

</div>

Okay I got it a bit faster:

Initial version: 2.363006s

1.0s → Insert forgotten points (meaning only partially fused kernels, and tmp allocations)  
0.7s → using sum(x) instead of sum(x, 1), possibly doesn’t generalize anymore? But the current version was 1D anyways  
0.4s → moving everything into one gpu kernel

final version:

```julia

function clstep!(state, randstate, ratesum, kappa, y_current, rates, dt)
    i = @linearidx(y_current, state)
    yi = y_current[i] # create tmp variables to not needlessly access array multiple times
    ytmp = yi + dt * (F(yi) + kappa)
    rate = (rates[i] = f(ytmp))
    rnd = GPUArrays.gpu_rand(Float32, state, randstate)
    y_current[i] = ytmp * (rate < -log(rnd) / dt)
    return
end

function euler_cl(T::TY,dt::TY,n::Int32,y0, plot_ = false)
    println("\n--> Euler algo, up to T = ",n*dt)
    n = min(n,Int(ceil(T/dt)))
	println("OK")
    y_current = copy(y0)
    rndcl = copy(y0) # to hold randomness
    rates = f.(y_current)
	println("OK")
    kappa = sqrt.(2 * H0 * sum(rates)/N) .* w
    nb_jumps = 0
	y_current .= y_current .+ dt .* (F.(y_current) .+ kappa)
	println("OK for Euler step")

    # variables for saving data
    tout = TY[]
    a = TY[]
    randstate = GPUArrays.cached_state(y_current)
    for ii=1:n
        ratesum = sum(rates)
        kappa = sqrt(2f0 * H0 * ratesum ./ N) .* w
        gpu_call(clstep!, y_current, (randstate, ratesum, kappa, y_current, rates, dt))
        # save data
        push!(tout, ii*dt)
    end
    return tout, a
end

#######################################################################################
xc0 = rand(TY,N) |> CLArray;

parms = TY.([0., 0.])
dt = TY(1. / (5.4N))
t,a_tot = euler_cl(TY(5.5),dt*2,Int32(10),xc0);

t,a_tot = @time euler_cl(TY(5.5), dt,Int32(1000),xc0);

```

I removed the `push!(a, kappa)` since there was a type mismatch… I guess that shouldn’t change anything for the benchmarks…

Best,  
Simon

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [December 2, 2017, 2:45pm UTC](https://discourse.julialang.org/t/saving-values-from-gpu-during-euler-stepping-in-ode/7406/18 "2017-12-02T14:45:53Z")

</div>

> @sdanisch a bit?

😍

---

<div class="post-metadata">

**Author:** ![sdwfrost](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sdwfrost/32/2831_2.png) [@sdwfrost](https://discourse.julialang.org/u/sdwfrost)\
**Post date:** [December 4, 2017, 8:19pm UTC](https://discourse.julialang.org/t/saving-values-from-gpu-during-euler-stepping-in-ode/7406/19 "2017-12-04T20:19:51Z")

</div>

I thought it would be straightforward to port this example to CuArrays from CLArrays, but I’m getting some errors that a GPU newbie like myself can’t easily resolve. Is there a guide to go between the two?

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [December 4, 2017, 8:21pm UTC](https://discourse.julialang.org/t/saving-values-from-gpu-during-euler-stepping-in-ode/7406/20 "2017-12-04T20:21:54Z")

</div>

> [@sdanisch](#):
>
> rate = (rates[i] = f(ytmp))

Do you really mean it or `rates = (rates[i] = f(ytmp))`?

[Next page](https://discourse.julialang.org/t/saving-values-from-gpu-during-euler-stepping-in-ode/7406.md?page=2)