# Moving from \`Float64\` to \`Float32\` not improving performance

**URL:** <https://discourse.julialang.org/t/moving-from-float64-to-float32-not-improving-performance/88516>\
**Category:** Performance\
**Created:** [October 10, 2022, 9:43am UTC](https://discourse.julialang.org/t/moving-from-float64-to-float32-not-improving-performance/88516 "2022-10-10T09:43:55Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![JordiBolibar](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jordibolibar/32/24307_2.png) [@JordiBolibar](https://discourse.julialang.org/u/JordiBolibar)\
**Post date:** [October 10, 2022, 9:43am UTC](https://discourse.julialang.org/t/moving-from-float64-to-float32-not-improving-performance/88516/1 "2022-10-10T09:43:55Z")

</div>

Hi, I’m trying to optimize a script. I first did the usual stuff to avoid memory allocation, but when moving all data structures from `Float64` to `Float32`, which should result in a reduced memory usage, the code has become super slow.

I’m still trying to get my head around this, but I was wondering if anyone could help me spot where my mistakes are.

- [Here](https://github.com/ODINN-SciML/iceflow_sandbox/blob/main/1D_SIA_raw.jl)’s the `Float64` version. Benchmarked time is `0.392s`.  
[Here](https://github.com/ODINN-SciML/iceflow_sandbox/blob/main/1D_SIA_raw_Float32.jl)’s the `Float32` version. The simulation takes forever, hinting at an error during the code conversion.

Here’s an overview of the `Float64` script to be optimized:

```julia
using Revise, BenchmarkTools
using ProgressMeter, Infiltrator
using Plots

const sec_in_day = 24.0 * 60.0 * 60.0
const sec_in_year = sec_in_day * 365.0
const glen_n = 3.0
const ρ = 900.0
const g = 9.81

@views diff1(A) = (A[2:end] .- A[1:end - 1])
@views diff2(A) = (A[3:end] .- A[1:end - 2])
@views avg(A) = (A[2:end] .+ A[1:end - 1])./2.0

function glacier_evolution_optim(;
    dx=100.0, # grid resolution in m
    nx=200, # grid size
    width=600.0, # glacier width in m
    top_h=3000.0, # bed top altitude
    bottom_h=1200.0, # bed bottom altitude
    glen_a=2.4e-24, # ice stiffness
    ela_h=2600.0, # mass balance model Equilibrium Line Altitude
    mb_grad=3.0, # linear mass balance gradient (unit: [mm w.e. yr-1 m-1])
    n_years=700 # simulation time in years
)
    function get_mb(heights)
        mb = (heights .- ela_h) .* mb_grad
        return mb ./ sec_in_year ./ ρ
    end

    let 
    bed_h = collect(LinRange(top_h, bottom_h, nx))
    surface_h = copy(bed_h)
    thick = bed_h .* 0.0

    t = 0.0
    dt = sec_in_day * 10.0

    years = collect(0:(n_years+1))
    volume = zeros(size(years))
    length = zeros(size(years))

    new_thick = zeros(nx)
    diffusivity = zeros(nx)
    diffusivity_s = zeros(nx)
    surface_gradient = zeros(nx)
    surface_gradient_s = zeros(nx)
    grad_x_diff = zeros(nx)
    flux_div = zeros(nx-1)
    mb = zeros(nx-1)

    for (i, y) in enumerate(years)
        let end_t = y * sec_in_year
        # Time integration
        while t < end_t
            # This is to guarantee a precise arrival on a specific date if asked
            remaining_t = end_t - t
            if remaining_t < dt
                dt = remaining_t
            end

            # Surface gradient
            surface_gradient[2:end-1] .= diff2(surface_h) ./ (2.0*dx)

            # Diffusivity
            diffusivity .= width * ((ρ*g)^3.0) .* (thick.^3.0) .* surface_gradient.^2.0
            diffusivity .*= 2.0/(glen_a+2.0) * glen_a .* thick.^2.0

            # Ice flux in a staggered grid
            diffusivity_s[2:end] .= avg(diffusivity)

            surface_gradient_s[2:end] .= diff1(surface_h) ./ dx

            grad_x_diff .= surface_gradient_s .* diffusivity_s
            flux_div .= diff1(grad_x_diff) ./ dx

            # Mass balance
            mb .= get_mb(surface_h[begin:end-1])

            # Ice thickness update: old + flux div + mb
            new_thick[begin:end-1] .= thick[begin:end-1] .+ (dt/width) .* flux_div .+ dt.*mb

            # We can have negative thickness because of MB - correct here
            thick .= ifelse.(new_thick.<0.0, 0.0, new_thick)
            
            @assert thick[end] == 0.0 "Glacier exceeding boundaries! at time $(t/sec_in_year)"

            # Prepare for next step 
            surface_h .= bed_h .+ thick
            t += dt
        end
        end # let

        volume[i] = sum(thick .* width .* dx)
        length[i] = sum(thick .> 0.0) .* dx
    end

    # xcoordinates
    xc = collect(0:nx-1) .* dx

    return xc, bed_h, surface_h, years, volume, length
    end # let
    
end

function wrapper_grad(grad)
    return glacier_evolution_optim(mb_grad=grad)
end

####### MAIN ########

# Let's test the performance
@btime xc, bed_h, surface_h, years, volume, length = glacier_evolution_optim()

```

Thanks in advance!

---

<div class="post-metadata">

**Author:** ![oheil](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oheil/32/220745_2.png) [@oheil](https://discourse.julialang.org/u/oheil)\
**Post date:** [October 10, 2022, 10:19am UTC](https://discourse.julialang.org/t/moving-from-float64-to-float32-not-improving-performance/88516/2 "2022-10-10T10:19:44Z")

</div>

If you do not run this on a GPU you will not gain ~~anything~~ much from Float32. See e.g. these discussion about it [Why no Float type alias?](https://discourse.julialang.org/t/why-no-float-type-alias/87947)

I just took a glance into your code and stumbled over:

- Why is function `get_mb(heights)` local in function `function glacier_evolution_optim`. Why not just a global function with additional parameters `ela_h` and `mb_grad` ?

- And I see those collect`s like

```julia
years = collect(0:(n_years+1))

```

Later you do

```julia
for (i, y) in enumerate(years)

```

but you don’t use `i` (perhaps I missed it?), so perhaps just a

```julia
years = 0:(n_years+1)
for y in years
...

```

will do.

- The purpose of those `let` 's is not clear to me. Perhaps you developed those parts in the REPL where you needed them for scope issues. You don’t need them in the function body.

I don’t know if those issues help for performance, I didn’t try to run your code. But cleaning it up from not needed structures will help a lot to get people into the maximize-performance-game.

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [October 10, 2022, 10:51am UTC](https://discourse.julialang.org/t/moving-from-float64-to-float32-not-improving-performance/88516/3 "2022-10-10T10:51:19Z")

</div>

It’s not clear where you are using `Float32`, but all your literal floats are `Float64`, everywhere you write `24.0` or `2.0` it’s a `Float64`, so you may get lots of type promotions or even instabilities.

Also, it’s better not to use floats for integer concepts. Replace

> [@JordiBolibar](#):
>
> `const sec_in_day = 24.0 * 60.0 * 60.0`

with

```julia
const sec_in_day = 24 * 60 * 60

```

You can see it here:

```julia
julia> x = 3f0 # Float32
3.0f0

julia> typeof(x * 2) # preserves type
Float32

julia> typeof(x * 2.0) # promotes to Float64
Float64

```

Multiplying with 2.0 promotes your variable to `Float64`, multiply with `2`, and the float type is preserved.

This:

```julia
diffusivity .= width * ((ρ*g)^3.0) .* (thick.^3.0) .* surface_gradient.^2.0
diffusivity .*= 2.0/(glen_a+2.0) * glen_a .* thick.^2.0

```

should be replaced with:

```julia
diffusivity .= width * ((ρ*g)^3) .* (thick.^3) .* surface_gradient.^2
diffusivity .*= 2/(glen_a+2) * glen_a .* thick.^2

```

Not just to avoid promotion, but because `x^2` is much faster than `x^2.0`. The former is just `x*x`, and `x^3` is `x*x*x`, while the latter needs to run code that calculates general floating point exponents.

Your loop over `years` seems to allocate a lot of temporary arrays. Can that be avoided? For example:

> [@JordiBolibar](#):
>
> ```julia
> volume[i] = sum(thick .* width .* dx)
> length[i] = sum(thick .> 0.0) .* dx
> 
> ```

could be changed to

```julia
volume[i] = dot(thick, width) * dx
length[i] = count(>(0), thick) * dx

```

---

<div class="post-metadata">

**Author:** ![JordiBolibar](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jordibolibar/32/24307_2.png) [@JordiBolibar](https://discourse.julialang.org/u/JordiBolibar)\
**Post date:** [October 10, 2022, 11:27am UTC](https://discourse.julialang.org/t/moving-from-float64-to-float32-not-improving-performance/88516/4 "2022-10-10T11:27:30Z")

</div>

Thanks for the feedback.

- Is there really a performance gain by moving a function to the global scope?
- In the for loop, I do use the `i` in `enumerate(years)`. So I can’t really remove that.
- The let block is to be able to access all those variables in the `for` and `while` loops scope. These introduce new local scopes, so without it I cannot access them.

To be honest, I’m not terribly worried about pushing the performance much further. I was mostly curious to understand why I wasn’t getting any performance gains by moving to `Float32`. After reading the discussion there it is still unclear to me why there shouldn’t be any benefits in using `Float32`. Shouldn’t it halve the memory allocation?

---

<div class="post-metadata">

**Author:** ![JordiBolibar](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jordibolibar/32/24307_2.png) [@JordiBolibar](https://discourse.julialang.org/u/JordiBolibar)\
**Post date:** [October 10, 2022, 11:31am UTC](https://discourse.julialang.org/t/moving-from-float64-to-float32-not-improving-performance/88516/5 "2022-10-10T11:31:01Z")

</div>

You’re looking at the `Float64` script. I gave direct links to the `Float32` and `Float64` scripts. In the `Float32` I’m using `Float32` types throughout the script as you mentioned.

I didn’t know about the `Int` exponents, thanks for that. I’ll aslo try the trick for `volume` and `length`.

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [October 10, 2022, 11:35am UTC](https://discourse.julialang.org/t/moving-from-float64-to-float32-not-improving-performance/88516/6 "2022-10-10T11:35:09Z")

</div>

What is important is the _number_ of allocations, not the amount of memory used. The slowdown comes from having to create blocks of memory to store data, and that does not change (much?) with the type of variable.

You can use `enumerate` without collecting the range.

If you want to accelerate the code, you need to get rid of temporary allocations. For example, this function:

> [@JordiBolibar](#):
>
> `@views diff2(A) = (A[3:end] .- A[1:end - 2])`

It returns a new array, and that makes the line where you use it allocating:

```julia

julia> @views diff2(A) = (A[3:end] .- A[1:end - 2])
diff2 (generic function with 3 methods)

julia> f(surface_gradient,surface_h,dx) = surface_gradient[2:end-1] .= diff2(surface_h) ./ (2.0*dx)
f (generic function with 2 methods)

julia> @btime f($surface_gradient, $surface_h, 0.1)
  36.695 ns (1 allocation: 128 bytes)

```

Now if you want that to be faster, you could update the `surface_gradient` differently, without creating the intermediate difference array:

```julia
julia> new_diff2!(surface_gradient, surface_h, dx) = @views @. surface_gradient[2:end-1] = (surface_h[3:end] - surface_h[1:end-2])/(2*dx)
new_diff2! (generic function with 1 method)

julia> @btime new_diff2!($surface_gradient, $surface_h, 0.1)
  18.282 ns (0 allocations: 0 bytes)

```

(I’m just using random arrays of length 10 here). If you solved these types of allocations you may obtain a much faster code.

`Float32`s can make the computations slightly faster at the end, but that will probably be very minor relative to solving these other issues.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [October 10, 2022, 11:46am UTC](https://discourse.julialang.org/t/moving-from-float64-to-float32-not-improving-performance/88516/7 "2022-10-10T11:46:56Z")

</div>

> [@JordiBolibar](#):
>
> `surface_gradient[2:end-1] .= diff2(surface_h) ./ (2.0*dx)`

Note this this allocates an extra temporary array for the result of `diff2`, because dot operations inside the function call are not merged with the dot operation outside the function call. Much better to do (with `@views` somewhere, e.g. on the whole `function`:

```julia
surface_gradient[2:end-1] .= (A[3:end] .- A[1:end-2]) ./ (2dx)

```

Note that it’s not necessary to write `2.0*dx` — `2dx` will automatically promote the result of the multiplication to the type of `dx`. Or even:

```julia
surface_gradient[2:end-1] .= (A[3:end] .- A[1:end-2]) .* eltype(A)(1/2dx)

```

so that the division is done outside the loop and is converted to the element type of `A` (e.g. `Float64`)

As another example:

> [@JordiBolibar](#):
>
> ```julia
> diffusivity .= width * ((ρ*g)^3.0) .* (thick.^3.0) .* surface_gradient.^2.0
> diffusivity .*= 2.0/(glen_a+2.0) * glen_a .* thick.^2.0
> 
> ```

has multiple problems.

1. Floating-point exponentials like `x^3.0` are _much_ slower to compute than integer exponents like `x^3`. The compiler doesn’t look at the _value_ of `3.0`, it just looks at the _type_ and sees `x^(somefloat)`, and calls a generic exponentiation routine that handles arbitrary non-integer exponents.
2. The `width *` is a non-dotted operator that allocates a new array for its result (when scaling an array)
3. By splitting this into two assignment statements, you have two loops instead of one.

Better to do something like:

```julia
diffusivity .= (width * (ρ*g)^3) .* thick.^3 .* surface_gradient.^2 .*
               (2glen_a/(glen_a+2)) .* thick.^2

```

Notice that I’m avoiding dots for the computation of purely scalar coefficients like `(2glen_a/(glen_a+2))` — this way, the scalar computation is done _outside_ the vectorized loop and the _result_ is dropped into the vectorized loop over the arrays.

There are many other similar problems. In general, it looks like you really need to understand (a) Julia’s type-promotion rules and (b) what the [dot vectorization syntax](https://julialang.org/blog/2017/01/moredots/) does and does not do.

---

<div class="post-metadata">

**Author:** ![JordiBolibar](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jordibolibar/32/24307_2.png) [@JordiBolibar](https://discourse.julialang.org/u/JordiBolibar)\
**Post date:** [October 10, 2022, 12:07pm UTC](https://discourse.julialang.org/t/moving-from-float64-to-float32-not-improving-performance/88516/8 "2022-10-10T12:07:59Z")

</div>

Thanks a lot for all this information. I’m aware I still have some important concepts to grasp correctly, but I also feel there are some grey areas which are not too clear to me.

@lmiq Thanks a lot for the suggestion. I updated the implementation of the functions and it helped boost the performance further.

@stevengj Thanks a lot for all the details. I’m a little bit confused with your comments on dot operations. I thought I was using those in the way you mention (i.e. no broadcasting for scalars, just when arrays are involved).

In your re-writing of this:

```julia
diffusivity .= (ρ*g)^3 .* thick.^3 .* surface_gradient.^2 .*
               (2glen_a/(glen_a+2)) .* thick^2

```

There are a few things I don’t understand. Why are you using `.^` in `thick` at the beginning but ommitting it in the second line? It is always an array. I understand that the way it was written before it could be improved in two ways: (1) using `Int` in the exponents, and (2) writing everything in a single line. But I’m was under the impression that I was already using broadcasting only when arrays were involved. Am I missing something?

Thanks again for all your feedback.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [October 10, 2022, 12:15pm UTC](https://discourse.julialang.org/t/moving-from-float64-to-float32-not-improving-performance/88516/9 "2022-10-10T12:15:55Z")

</div>

> [@JordiBolibar](#):
>
> Why are you using `.^` in `thick` at the beginning but ommitting it in the second line?

Typo, sorry. It’s now fixed.

> [@JordiBolibar](#):
>
> I was already using broadcasting only when arrays were involved

`width * (...)` in your code multiplies a scalar `width` times an array. While this is a valid operation, because it is not “dotted” it generates a new temporary array as a result, rather than “fusing” the scalar multiplication with the other dot operations into a single allocation-free loop.

`2.0/(glen_a+2.0) * glen_a .* thick.^2.0` is okay because `*` in Julia left-associative, so the scalar operations at left are done first, before the `.*`. But relying on associativity like this can be a bit fragile because it is vulnerable to slight code modifications. It’s safer to write `(2glen_a / (glen_a+2))` in parentheses to make it clear that it’s supposed to be calculated separately, in my opinion.

---

<div class="post-metadata">

**Author:** ![JordiBolibar](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jordibolibar/32/24307_2.png) [@JordiBolibar](https://discourse.julialang.org/u/JordiBolibar)\
**Post date:** [October 10, 2022, 12:21pm UTC](https://discourse.julialang.org/t/moving-from-float64-to-float32-not-improving-performance/88516/10 "2022-10-10T12:21:58Z")

</div>

Hmmm, in this line:

```julia
diffusivity .= width * ((ρ*g)^3) .* (thick.^3) .* surface_gradient.^2

```

I was following this logic: `width`, `ρ` and `g` are scalars, so all of them interact with non-dotted operators. However, the interface with that scalar block (i.e. `width * ((ρ*g)^3.0)`) with an array, uses dotted operators. Shouldn’t that be correct as well, following Julia’s left-association?

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [October 10, 2022, 12:26pm UTC](https://discourse.julialang.org/t/moving-from-float64-to-float32-not-improving-performance/88516/11 "2022-10-10T12:26:43Z")

</div>

That line is fine:

```julia
julia> f(diffusivity, width, ρ, g, thick, surface_gradient) = diffusivity .= width * ((ρ*g)^3) .* (thick.^3) .* surface_gradient.^2
f (generic function with 1 method)

julia> diffusivity = rand(10); thick=rand(10); surface_gradient=rand(10);

julia> @btime f($diffusivity, 0.1, 0.1, 0.1, $thick, $surface_gradient)
  9.385 ns (0 allocations: 0 bytes)

```

IMHO, independently on how much one internalizes the rules, having the habit of testing these code snippets, like I did there, is what will make you confident that you are doing the best on each line.

---

<div class="post-metadata">

**Author:** ![JordiBolibar](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jordibolibar/32/24307_2.png) [@JordiBolibar](https://discourse.julialang.org/u/JordiBolibar)\
**Post date:** [October 10, 2022, 12:36pm UTC](https://discourse.julialang.org/t/moving-from-float64-to-float32-not-improving-performance/88516/12 "2022-10-10T12:36:20Z")

</div>

I completely agree here. That’s what I normally do. But I must admit that I often come across counter intuitive behaviours. For example, some of you suggested merging these two lines in a single line in order to avoid having two loops:

```julia
diffusivity .= width * ((ρ*g)^3.0) .* (thick.^3.0) .* surface_gradient.^2.0
diffusivity .*= 2.0/(glen_a+2.0) * glen_a .* thick.^2.0

```

However, the benchmark indicates that it is actually faster to keep the code like this in two lines. I sometimes struggle to find a reasonable explanation to some practical performance results.

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [October 10, 2022, 12:42pm UTC](https://discourse.julialang.org/t/moving-from-float64-to-float32-not-improving-performance/88516/13 "2022-10-10T12:42:40Z")

</div>

note that as of Julia 1.8, x^y has a fast path for when y is an integer. it’s still better to use integer exponents directly, but it’s not as bad as it used to be.

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [October 10, 2022, 12:44pm UTC](https://discourse.julialang.org/t/moving-from-float64-to-float32-not-improving-performance/88516/14 "2022-10-10T12:44:30Z")

</div>

If you manage to write correctly the dots, the single line may be slightly faster, because you will be saving a loop. But that dot fusion syntax can become pretty hard to get right.

I tend to avoid these too complicated dot combinations, because of that, but it is a matter of taste. Sometimes writing the loop is just simpler:

```julia
julia> function f0(diffusivity, width, ρ, g, thick, surface_gradient, glen_a)
           diffusivity .= width * ((ρ*g)^3.0) .* (thick.^3.0) .* surface_gradient.^2.0
           diffusivity .*= 2.0/(glen_a+2.0) * glen_a .* thick.^2.0
           return diffusivity
       end
f0 (generic function with 1 method)

julia> diffusivity = rand(10); thick=rand(10); surface_gradient=rand(10);

julia> @btime f0($diffusivity, 1.0, 1.0, 1.0, $thick, $surface_gradient, 1.0);
  75.667 ns (0 allocations: 0 bytes)

julia> function f1(diffusivity, width, ρ, g, thick, surface_gradient, glen_a)
           for i in eachindex(diffusivity, thick, surface_gradient)
               diffusivity[i] = (width * ((ρ*g)^3.0) * (thick[i]^3.0) * surface_gradient[i]^2.0) *
                                (2.0/(glen_a+2.0) * glen_a * thick[i]^2.0)
           end
           return diffusivity
       end

f1 (generic function with 1 method)

julia> @btime f1($diffusivity, 1.0, 1.0, 1.0, $thick, $surface_gradient, 1.0);
  70.115 ns (0 allocations: 0 bytes)

```

ps: By fixing the exponents to integers, this is what we get:

```julia
julia> function f1(diffusivity, width, ρ, g, thick, surface_gradient, glen_a)
           for i in eachindex(diffusivity, thick, surface_gradient)
               diffusivity[i] = (width * ((ρ*g)^3) * (thick[i]^3) * surface_gradient[i]^2) *
                                          (2/(glen_a+2) * glen_a * thick[i]^2)
           end
           return diffusivity
       end
f1 (generic function with 1 method)

julia> @btime f1($diffusivity, 1.0, 1.0, 1.0, $thick, $surface_gradient, 1.0);
  9.594 ns (0 allocations: 0 bytes)

```

---

<div class="post-metadata">

**Author:** ![Palli](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/palli/32/3380_2.png) [@Palli](https://discourse.julialang.org/u/Palli)\
**Post date:** [October 10, 2022, 1:05pm UTC](https://discourse.julialang.org/t/moving-from-float64-to-float32-not-improving-performance/88516/15 "2022-10-10T13:05:37Z")

</div>

> [@oheil](#):
>
> If you do not run this on a GPU you will not gain anything from Float32.

I don’t think that’s true (in general). Yes, on GPUs Float32 is going to be a lot faster, and again Float16 (blfoat16 etc.), because you have more functional units for those (and better use of memory [bandwidth]).

> [@JordiBolibar](#):
>
> when moving all data structures from `Float64` to `Float32`, which should result in a reduced memory usage, the code has become super slow.

I believe that should not happen, on CPUs (or GPUs), it means you are doing something wrong, likely introduced a type-instability (noticeable by e.g. extra memory allocations). If you do it right I think at worst if should be as fast (though less precise, not always a problem, but if you need to work around that, then could get slower, but you haven’t tried that yet).

One exception I can think of where Float32 would make slower, is if you converge slower, might apply for you, but I would rule out other issues people have already pointed out first.

At best Float32 will give you 2x speedup on CPUs (and 4x for Float16, while often not achievable), because of memory bandwidth, right?

I see in your code, the first needless line:

```julia-auto
# -*- coding: utf-8 -*-

```

That isn’t needed/wanted in Julia, nor even in Python 3. Julia always uses UTF-8 for source code (and this is just a comment in Julia), and Python 3 defaults to it (do people ever use anything else by now? You could use such a special comment in Python for non-UTF-8,or for UTF-8 in outdated Python 2.

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [October 10, 2022, 1:10pm UTC](https://discourse.julialang.org/t/moving-from-float64-to-float32-not-improving-performance/88516/16 "2022-10-10T13:10:35Z")

</div>

> [@JordiBolibar](#):
>
> `diffusivity .= width * ((ρ*g)^3) .* (thick.^3) .* surface_gradient.^2`

The above works ok, but I suggest instead:

```julia
diffusivity .= (width * (ρ*g)^3) .* (thick.^3) .* surface_gradient.^2

```

It’s clearer that ` (width * (ρ*g)^3)` is a group.

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [October 10, 2022, 1:10pm UTC](https://discourse.julialang.org/t/moving-from-float64-to-float32-not-improving-performance/88516/17 "2022-10-10T13:10:46Z")

</div>

Concerning the original question, F32 can accelerate, if using LoopVectorization:

```julia
julia> using LoopVectorization

julia> function f1(diffusivity, width, ρ, g, thick, surface_gradient, glen_a)
           @turbo for i in eachindex(diffusivity, thick, surface_gradient)
               diffusivity[i] = (width * ((ρ*g)^3) * (thick[i]^3) * surface_gradient[i]^2) *
                                          (2/(glen_a+2) * glen_a * thick[i]^2)
           end
           return diffusivity
       end
f1 (generic function with 1 method)

julia> diffusivity = rand(10); thick=rand(10); surface_gradient=rand(10);

julia> @btime f1($diffusivity, 1.0, 1.0, 1.0, $thick, $surface_gradient, 1.0);
  11.353 ns (0 allocations: 0 bytes)

julia> diffusivity = rand(Float32,10); thick=rand(Float32,10); surface_gradient=rand(Float32,10);

julia> @btime f1($diffusivity, 1.0f0, 1.0f0, 1.0f0, $thick, $surface_gradient, 1.0f0);
  5.362 ns (0 allocations: 0 bytes)

```

But it is not about memory usage, but how many operations the processors can do simultaneously.

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [October 10, 2022, 1:11pm UTC](https://discourse.julialang.org/t/moving-from-float64-to-float32-not-improving-performance/88516/18 "2022-10-10T13:11:26Z")

</div>

> [@JordiBolibar](#):
>
> You’re looking at the `Float64` script.

I think it would have made more sense to share the `Float32` script in the OP instead.

Also, the link to the `Float32` is broken.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [October 10, 2022, 1:13pm UTC](https://discourse.julialang.org/t/moving-from-float64-to-float32-not-improving-performance/88516/19 "2022-10-10T13:13:48Z")

</div>

> [@JordiBolibar](#):
>
> Hmmm, in this line:
> 
> ```julia
> diffusivity .= width * ((ρ*g)^3) .* (thick.^3) .* surface_gradient.^2
> 
> ```
> 
> I was following this logic: `width`, `ρ` and `g` are scalars, so all of them interact with non-dotted operators. However, the interface with that scalar block (i.e. `width * ((ρ*g)^3.0)`) with an array, uses dotted operators. Shouldn’t that be correct as well, following Julia’s left-association?

Oh, I was mis-reading the parentheses. Yes, the `width * ((ρ*g)^3) ` multiplication will be done first as a scalar, and _then_ multiplied by the array. Sorry.

(But I still think it’s clearer to explicitly group scalar operations with parentheses rather than relying on associativity. e.g. if you merge the two lines then the associativity may get changed if you are not careful.)

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [October 10, 2022, 1:20pm UTC](https://discourse.julialang.org/t/moving-from-float64-to-float32-not-improving-performance/88516/20 "2022-10-10T13:20:57Z")

</div>

> [@JordiBolibar](#):
>
> However, the benchmark indicates that it is actually faster to keep the code like this in two lines

I can’t reproduce:

```julia
function f1(diffusivity, width, ρ, g, thick, surface_gradient, glen_a)
    diffusivity .= width * ((ρ*g)^3) .* (thick.^3) .* surface_gradient.^2
    diffusivity .*= 2.0/(glen_a+2.0) * glen_a .* thick.^2
end

function f2(diffusivity, width, ρ, g, thick, surface_gradient, glen_a)
    diffusivity .= (width * (ρ*g)^3) .* thick.^3 .* surface_gradient.^2 .*
                   (2glen_a/(glen_a+2)) .* thick.^2
end

using BenchmarkTools
for n in (10, 100, 1000, 10000)
    @show n
    @btime f1($(rand(n)), 1.0, 1.0, 1.0, $(rand(n)), $(rand(n)), 1.0)
    @btime f2($(rand(n)), 1.0, 1.0, 1.0, $(rand(n)), $(rand(n)), 1.0)
end

```

gives

```julia
n = 10
  12.929 ns (0 allocations: 0 bytes)
  11.053 ns (0 allocations: 0 bytes)
n = 100
  44.276 ns (0 allocations: 0 bytes)
  38.936 ns (0 allocations: 0 bytes)
n = 1000
  375.607 ns (0 allocations: 0 bytes)
  322.017 ns (0 allocations: 0 bytes)
n = 10000
  4.494 μs (0 allocations: 0 bytes)
  3.188 μs (0 allocations: 0 bytes)

```

i.e. `f2` is consistently faster, especially as the array gets larger and cache locality becomes an issue.

[Next page](https://discourse.julialang.org/t/moving-from-float64-to-float32-not-improving-performance/88516.md?page=2)