# Something faster than for loops

**URL:** <https://discourse.julialang.org/t/something-faster-than-for-loops/24001>\
**Category:** General Usage\
**Created:** [May 8, 2019, 4:45pm UTC](https://discourse.julialang.org/t/something-faster-than-for-loops/24001 "2019-05-08T16:45:05Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![Aquaman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aquaman/32/6586_2.png) [@Aquaman](https://discourse.julialang.org/u/Aquaman)\
**Post date:** [May 8, 2019, 4:45pm UTC](https://discourse.julialang.org/t/something-faster-than-for-loops/24001/1 "2019-05-08T16:45:06Z")

</div>

Dear All,

I am an “old” programmer used to Fortran and Matlab. I can’t get rid of For and While loops, but I know Julia can do things much faster!

I am using the following code

```julia
n=5
p=5
m=5

A=rand(n,p)
B=rand(p,m)
C=zeros(Float64,n,m)
for i in 1:n

  for j in 1:m
     C[i,j]=0.0
     for k in 1:p

        C[i,j] = C[i,j]+A[i,k]*B[k,j]
     end
   end
end

```

If n is big, it slows down. e.g. n \> 10000.

How can the code be made faster in Julia?

---

<div class="post-metadata">

**Author:** ![rdeits](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rdeits/32/286_2.png) [@rdeits](https://discourse.julialang.org/u/rdeits)\
**Post date:** [May 8, 2019, 4:49pm UTC](https://discourse.julialang.org/t/something-faster-than-for-loops/24001/2 "2019-05-08T16:49:05Z")

</div>

Loops are just as fast in Julia as they are in Fortran, and a well-written loop is often the single fastest way to implement an algorithm.

What is your code actually supposed to do? Your inner loop over `k in 1:p` just does the same operation over and over. Removing that loop would obviously speed up your code, but I suspect that isn’t what you want.

Also, have you read the Julia [performance tips](https://docs.julialang.org/en/v1/manual/performance-tips/index.html)? That’s the first place to start to learn how to make your Julia code faster.

---

<div class="post-metadata">

**Author:** ![Aquaman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aquaman/32/6586_2.png) [@Aquaman](https://discourse.julialang.org/u/Aquaman)\
**Post date:** [May 8, 2019, 4:57pm UTC](https://discourse.julialang.org/t/something-faster-than-for-loops/24001/3 "2019-05-08T16:57:59Z")

</div>

Hi,

it is just an example.

if n=m=p= **5000**

```julia
n=5000
p=5000
m=5000

A=rand(n,p);
B=rand(p,m);
C=zeros(Float64,n,m)
for i in 1:n

  for j in 1:m
     C[i,j]=0.0
     for k in 1:p

        C[i,j] = C[i,j]+A[i,j]*B[i,j]
     end
   end
end

```

The simulation takes **so long** …

---

<div class="post-metadata">

**Author:** ![davidbp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/davidbp/32/463_2.png) [@davidbp](https://discourse.julialang.org/u/davidbp)\
**Post date:** [May 8, 2019, 4:58pm UTC](https://discourse.julialang.org/t/something-faster-than-for-loops/24001/4 "2019-05-08T16:58:17Z")

</div>

I would add that… are you timing your code as it is shown or do you have it inside a function?

```julia

function f()
   n=5
   p=5
   m=5

   A=rand(n,p)
   B=rand(p,m)
   C=zeros(Float64,n,m)
   for i in 1:n
      for j in 1:m
         C[i,j]=0.0
         for k in 1:p

            C[i,j] = C[i,j]+A[i,j]*B[i,j]
         end
       end
    end
end

@time f()

```

This is what I get

```julia
julia> @time f()
  0.033329 seconds (59.47 k allocations: 3.133 MiB)

julia> @time f()
  0.000006 seconds (7 allocations: 1024 bytes)

```

---

<div class="post-metadata">

**Author:** ![Aquaman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aquaman/32/6586_2.png) [@Aquaman](https://discourse.julialang.org/u/Aquaman)\
**Post date:** [May 8, 2019, 5:01pm UTC](https://discourse.julialang.org/t/something-faster-than-for-loops/24001/5 "2019-05-08T17:01:16Z")

</div>

> [@Aquaman](#):
>
> if n=m=p= **5000**

?

---

<div class="post-metadata">

**Author:** ![e3c6](https://avatars.discourse-cdn.com/v4/letter/e/e79b87/32.png) [@e3c6](https://discourse.julialang.org/u/e3c6)\
**Post date:** [May 8, 2019, 5:05pm UTC](https://discourse.julialang.org/t/something-faster-than-for-loops/24001/6 "2019-05-08T17:05:52Z")

</div>

You should try to switch the order of the loops. In Julia matrices are stored column-major, which means that to access `A[i,j]` linearly in memory (which is fastest), `i` be the variable that changes fastest.

Therefore, instead of

```julia
for i in 1:n
  for j in 1:m
    C[i,j]=0.0
    ....
  end
end

```

try this:

```julia
for j in 1:m
  for i in 1:n
    C[i,j]=0.0
    ....
  end
end

```

(the only change is that `i` is traversed in the inner loop now).

---

<div class="post-metadata">

**Author:** ![davidbp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/davidbp/32/463_2.png) [@davidbp](https://discourse.julialang.org/u/davidbp)\
**Post date:** [May 8, 2019, 5:08pm UTC](https://discourse.julialang.org/t/something-faster-than-for-loops/24001/7 "2019-05-08T17:08:38Z")

</div>

```julia
function f(n,p,m)

          A=rand(n,p)
          B=rand(p,m)
          C=zeros(Float64,n,m)
          for i in 1:n
             for j in 1:m
                C[i,j]=0.0
                for k in 1:p

                   C[i,j] = C[i,j]+A[i,j]*B[i,j]
                end
              end
           end
       end

julia> @time f(1000,1000,1000)
1.756890 seconds (10 allocations: 22.889 MiB, 2.62% gc time)

julia> @time f(2000,2000,2000)
27.841935 seconds (10 allocations: 91.553 MiB, 0.36% gc time)

julia> @time f(5000,5000,5000)
429.574190 seconds (10 allocations: 572.205 MiB, 0.02% gc time)

```

---

<div class="post-metadata">

**Author:** ![Aquaman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aquaman/32/6586_2.png) [@Aquaman](https://discourse.julialang.org/u/Aquaman)\
**Post date:** [May 8, 2019, 5:18pm UTC](https://discourse.julialang.org/t/something-faster-than-for-loops/24001/8 "2019-05-08T17:18:51Z")

</div>

> [@davidbp](#):
>
> end @time f()

**in case of n= 5000 …**

---

<div class="post-metadata">

**Author:** ![davidbp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/davidbp/32/463_2.png) [@davidbp](https://discourse.julialang.org/u/davidbp)\
**Post date:** [May 8, 2019, 5:21pm UTC](https://discourse.julialang.org/t/something-faster-than-for-loops/24001/9 "2019-05-08T17:21:13Z")

</div>

I just gave you the result!

```julia
429.574190 seconds (10 allocations: 572.205 MiB, 0.02% gc time)

```

🙂

---

<div class="post-metadata">

**Author:** ![Aquaman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aquaman/32/6586_2.png) [@Aquaman](https://discourse.julialang.org/u/Aquaman)\
**Post date:** [May 8, 2019, 5:23pm UTC](https://discourse.julialang.org/t/something-faster-than-for-loops/24001/10 "2019-05-08T17:23:08Z")

</div>

Thanks .🙂

---

<div class="post-metadata">

**Author:** ![lstagner](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lstagner/32/448_2.png) [@lstagner](https://discourse.julialang.org/u/lstagner)\
**Post date:** [May 8, 2019, 5:30pm UTC](https://discourse.julialang.org/t/something-faster-than-for-loops/24001/11 "2019-05-08T17:30:43Z")

</div>

You also want to try multithreading. Just adding `Threads.@threads` to the outermost loop significantly improved performance.

```julia
function f2(n,p,m)
    A=rand(n,p)
    B=rand(p,m)
    C=zeros(Float64,n,m)
    Threads.@threads for j in 1:m
        for i in 1:n
            C[i,j]=0.0
            for k in 1:p
                C[i,j] = C[i,j]+A[i,j]*B[i,j]
            end
         end
     end
 end

julia> @time f2(5000,5000,5000); # with JULIA_NUM_THREADS=12 
 52.803032 seconds (11 allocations: 572.205 MiB, 0.21% gc time)

```

F.Y.I. You need to set the environmental variable `JULIA_NUM_THREADS` before you start julia

---

<div class="post-metadata">

**Author:** ![davidbp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/davidbp/32/463_2.png) [@davidbp](https://discourse.julialang.org/u/davidbp)\
**Post date:** [May 8, 2019, 5:33pm UTC](https://discourse.julialang.org/t/something-faster-than-for-loops/24001/12 "2019-05-08T17:33:58Z")

</div>

Even without multithreading if you avoid going to the arrays 125 billion times, you can improve the performance a lot

```julia
function f3(n,p,m)
   A=rand(n,p)
   B=rand(p,m)
   C=zeros(Float64,n,m)
   for i in 1:n
      for j in 1:m
         aux=0.0
         A_i_j = A[i,j]
         B_i_j = B[i,j]
         for k in 1:p
            aux += A_i_j *B_i_j
         end
         C[i,j]=aux
       end
    end
end

```

```julia
@time f3(5000,5000,5000)
127.116506 seconds (10 allocations: 572.205 MiB, 0.10% gc time)

```

---

<div class="post-metadata">

**Author:** ![kristoffer.carlsson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kristoffer.carlsson/32/22_2.png) [@kristoffer.carlsson](https://discourse.julialang.org/u/kristoffer.carlsson)\
**Post date:** [May 8, 2019, 5:36pm UTC](https://discourse.julialang.org/t/something-faster-than-for-loops/24001/13 "2019-05-08T17:36:33Z")

</div>

> [@lstagner](#):
>
> You also want to try multithreading

And `@inbounds` and `@simd`.

However, the toy example is starting to wear out its usefulness. There is nothing faster than loops since loops are just a way of jumping back to an earlier instruction. Instead one should talk about the operation that should be performed and how that can be done as efficiently as possible.

---

<div class="post-metadata">

**Author:** ![davidbp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/davidbp/32/463_2.png) [@davidbp](https://discourse.julialang.org/u/davidbp)\
**Post date:** [May 8, 2019, 5:39pm UTC](https://discourse.julialang.org/t/something-faster-than-for-loops/24001/14 "2019-05-08T17:39:34Z")

</div>

Just for fun

```julia
function f3(n,p,m)
   A=rand(n,p)
   B=rand(p,m)
   C=zeros(Float64,n,m)
   @inbounds for i in 1:n
      for j in 1:m
         aux=0.0
         A_i_j = A[i,j]
         B_i_j = B[i,j]
         @simd for k in 1:p
            aux += A_i_j *B_i_j
         end
         C[i,j]=aux
       end
    end
end

```

Single threaded code in my 2013 laptop

```julia
@time f3(5000,5000,5000)
 14.279610 seconds (10 allocations: 572.205 MiB, 0.69% gc time)

```

---

<div class="post-metadata">

**Author:** ![Aquaman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aquaman/32/6586_2.png) [@Aquaman](https://discourse.julialang.org/u/Aquaman)\
**Post date:** [May 8, 2019, 5:41pm UTC](https://discourse.julialang.org/t/something-faster-than-for-loops/24001/15 "2019-05-08T17:41:06Z")

</div>

Great, this is the answer which I need !

---

<div class="post-metadata">

**Author:** ![davidbp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/davidbp/32/463_2.png) [@davidbp](https://discourse.julialang.org/u/davidbp)\
**Post date:** [May 8, 2019, 5:54pm UTC](https://discourse.julialang.org/t/something-faster-than-for-loops/24001/16 "2019-05-08T17:54:56Z")

</div>

Remember that you can get a nice speed boost if you don’t need Float64 (going down to 7 sec in Float32). Even bigger if your CPU has AVX512.

```julia
function f5(n,p,m)
  A=rand(Float32, n,p)
  B=rand(Float32, p,m)
  C=zeros(Float32, n,m)
  @inbounds for i in 1:n
     for j in 1:m
        aux=Float32(0.0)
        A_i_j = A[i,j]
        B_i_j = B[i,j]
        @simd for k in 1:p
           aux += A_i_j *B_i_j
        end
        C[i,j]=aux
      end
   end
end

@time f5(5000,5000,5000)
  7.271317 seconds (10 allocations: 286.103 MiB, 1.38% gc time)

```

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [May 8, 2019, 5:56pm UTC](https://discourse.julialang.org/t/something-faster-than-for-loops/24001/17 "2019-05-08T17:56:44Z")

</div>

Why not `aux = p* A_i_j *B_i_j` instead?

But I think some things are wrong with this code. I don’t know what it’s supposed to do.  
It looks like the code is simply doing:

```julia
C .= p .* A .* B

```

But the dimensions don’t match up correctly. We just so happen to have that `n == m == p` in the examples. If they don’t, we could get segfaults, etc.

---

<div class="post-metadata">

**Author:** ![davidbp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/davidbp/32/463_2.png) [@davidbp](https://discourse.julialang.org/u/davidbp)\
**Post date:** [May 8, 2019, 5:57pm UTC](https://discourse.julialang.org/t/something-faster-than-for-loops/24001/18 "2019-05-08T17:57:41Z")

</div>

knowing what multiplication means is indeed useful. I guess the question was about how to use “julia syntax” not how to “rethink” the code to go from O(n^3) to O(n^2).

---

<div class="post-metadata">

**Author:** ![StefanKarpinski](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stefankarpinski/32/24_2.png) [@StefanKarpinski](https://discourse.julialang.org/u/StefanKarpinski)\
**Post date:** [May 8, 2019, 7:45pm UTC](https://discourse.julialang.org/t/something-faster-than-for-loops/24001/19 "2019-05-08T19:45:28Z")

</div>

Is this example trying to do matrix multiplication? That’s what the sizes of the input and output matrices would suggest…

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [May 8, 2019, 8:24pm UTC](https://discourse.julialang.org/t/something-faster-than-for-loops/24001/20 "2019-05-08T20:24:02Z")

</div>

I think so. That would explain the three loops.

In that case, @Aquaman:  
Just use the BLAS library; it will be many times faster than the loops.

```julia
# If C has not already been allocated
C = A * B
# If you pre-allocated C
using LinearAlgebra
mul!(C, A, B)

```

The BLAS libraries also use loops – but [there are actually five of them, not three](http://www.cs.utexas.edu/users/flame/BLISRetreat2018/images/BLIS_gemm.png)!\*  
They loop over an optimized microkernel. The optimized kernel is architecture dependent, so OpenBLAS (and other optimized BLAS libraries) will choose the correct one based on your CPU.  
They’re also multithreaded.

\*Although two of the loops are probably going to be completely unrolled.

[Next page](https://discourse.julialang.org/t/something-faster-than-for-loops/24001.md?page=2)