title: Speed in Matlab vs. Julia vs. Fortran
body: <p>I am playing around with different languages to solve a simple value function iteration problem where I loop over a state-space grid. I am trying to understand the performance differences and how I could tweak each code. For posterity I have posted full length working examples for each language below. However, I believe that most of the tweaking is to be done in the <code>while</code> loop. I am a bit confused what I am doing wrong in Fortran as the speed seems subpar.</p>
  <p><strong>Matlab ~2.7secs :</strong> I am avoiding a more efficient solution using the <code>repmat</code> function for now to keep the codes comparable. Code seems to be automatically multithreaded onto 4 threads</p>
  <pre><code>beta = 0.98;
  sigma = 0.5;
  R = 1/beta;
  
  a_grid = linspace(0,100,1001);
  tic
  [V_mat, next_mat] = valfun(beta, sigma, R ,a_grid);
  toc
  </code></pre>
  <p>where valfun()</p>
  <pre><code>function [V_mat, next_mat] = valfun(beta, sigma, R, a_grid)
      zeta = 1-1/sigma;
      len = length(a_grid);
      V_mat = zeros(2,len);
      next_mat = zeros(2,len);
  
      u = zeros(2,len,len);
      c = zeros(2,len,len);
  
      for i = 1:len
          c(1,:,i) = a_grid(i) - a_grid/R + 20.0;
          c(2,:,i) = a_grid(i) - a_grid/R;
      end
  
      u = c.^zeta * zeta^(-1);
      u(c&lt;=0) = -1e8;
  
      tol = 1e-4;
      outeriter = 0;
      diff = 1000.0;
  
       while (diff&gt;tol)  %&amp;&amp; (outeriter&lt;20000)
          outeriter = outeriter + 1;
          V_last = V_mat;
  
          for i = 1:len 
              [V_mat(1,i), next_mat(1,i)] = max( u(1,:,i) + beta*V_last(2,:));
              [V_mat(2,i), next_mat(2,i)] = max( u(2,:,i) + beta*V_last(1,:));
          end
          diff = max(abs(V_mat - V_last));
      end
  
      fprintf(&quot;\n Value Function converged in %i steps. \n&quot;, outeriter)  
  end
  </code></pre>
  <p><strong>Julia (after compilation) ~5.4secs (4 threads (9425469 allocations: 22.43 GiB)), ~7.8secs (1 thread (2912564 allocations: 22.29 GiB))</strong></p>
  <p><strong>[EDIT: after adding correct broadcasting and @views its only 1.8-2.1seconds now, see below!]</strong></p>
  <pre><code>using LinearAlgebra, UnPack, BenchmarkTools
  
  struct paramsnew
      β::Float64
      σ::Float64
      R::Float64
  end
  
  function valfun(params, a_grid)
      @unpack β,σ, R = params
      ζ = 1-1/σ
  
      len = length(a_grid)
      V_mat = zeros(2,len)
      next_mat = zeros(2,len)
      u = zeros(2,len,len)
      c = zeros(2,len,len)
  
      @inbounds for i in 1:len
          c[1,:,i] = @. a_grid[i] - a_grid/R .+ 20.0
          c[2,:,i] = @. a_grid[i] - a_grid/R
      end
  
      u = c.^ζ * ζ^(-1)
      u[c.&lt;=0] .= typemin(Float64)
  
      tol = 1e-4
      outeriter = 0
      test = 1000.0
  
      while test&gt;tol 
          outeriter += 1
          V_last = deepcopy(V_mat)
  
          @inbounds Threads.@threads for i in 1:len # loop over grid points
              V_mat[1,i], next_mat[1,i] = findmax( u[1,:,i] .+ β*V_last[2,:])
              V_mat[2,i], next_mat[2,i] = findmax( u[2,:,i] .+ β*V_last[1,:])
          end
  
          test = maximum( abs.(V_mat - V_last)[.!isnan.( V_mat - V_last )])
      end
  
      print(&quot;\n Value Function converged in &quot;, outeriter, &quot; steps.&quot;)
      return V_mat, next_mat
  end
  
  a_grid = collect(0:0.1:100)
  p1 = paramsnew(0.98, 1/2, 1/0.98);
  
  @time valfun(p1,a_grid)
  print(&quot;\n should be compiled now \n&quot;)
  
  @btime valfun(p1,a_grid)
  </code></pre>
  <p><strong>Fortran (O3, mkl, qopenmp) ~9.2secs:</strong> I also must be doing something wrong when declaring the <code>openmp</code> variables as the compilation will crash for some grid sizes when using <code>openmp</code> (SIGSEGV error).</p>
  <pre><code>module mod_calc
      use omp_lib
      implicit none
      integer, parameter :: dp = selected_real_kind(33,4931), len = 1001
      public :: dp, len
  
      contains
  
      subroutine linspace(from, to, array)
          real(dp), intent(in) :: from, to
          real(dp), intent(out) :: array(:)
          real(dp) :: range
          integer :: n, i
          n = size(array)
          range = to - from
  
          if (n == 0) return
  
          if (n == 1) then
              array(1) = from
              return
          end if
  
          do i=1, n
              array(i) = from + range * (i - 1) / (n - 1)
          end do
      end subroutine
  
      subroutine calc_val() 
          real(dp)::  bbeta, sigma, R, zeta, tol, test
          real(dp)::  a_grid(len), V_mat(2,len), V_last(2,len), &amp;
                      u(len,len,2), c(len,len,2)
          
          integer :: outeriter, i, sss, next_mat(2,len), fu
          character(len=*), parameter :: FILE_NAME = 'data.txt'   ! File name.
  
          call linspace(from=0._dp, to=100._dp, array=a_grid)
  
          bbeta = 0.98
          sigma = 0.5
          R = 1.0/0.98
          zeta = 1.0 - 1.0/sigma
  
          tol = 1e-4
          test = 1000.0
          outeriter = 0
  
          do i = 1,len
              c(:,i,1) = a_grid(i) - a_grid/R + 20.0
              c(:,i,2) = a_grid(i) - a_grid/R
          end do
  
          u = c**zeta * 1.0/zeta
  
          where (c&lt;=0)
              u = -1e6
          end where
  
          V_mat = 0.0
          next_mat = 0.0
  
          do while (test&gt;tol .and. outeriter&lt;20000)
              outeriter = outeriter+1
              V_last = V_mat
  
              !$OMP PARALLEL DEFAULT(NONE)          &amp;
              !$OMP SHARED(V_mat, next_mat,V_last, u, bbeta)  &amp;
              !$OMP PRIVATE(i)                      
              !$OMP DO SCHEDULE(static)
              do i=1,len 
                  V_mat(1,i)    = maxval(u(:,i,1) + bbeta*V_last(2,:))
                  next_mat(1,i) = maxloc(u(:,i,1) + bbeta*V_last(2,:),1)
                  V_mat(2,i)    = maxval(u(:,i,2) + bbeta*V_last(1,:))
                  next_mat(2,i) = maxloc(u(:,i,2) + bbeta*V_last(1,:),1)
              end do
              !$OMP END DO
              !$OMP END PARALLEL
  
              test = maxval(abs(log(V_last/V_mat)))
          end do
      end subroutine
  end module mod_calc
  
  program main
      use mod_calc
  
      implicit none
      integer:: clck_counts_beg,clck_rate,clck_counts_end
      
      call omp_set_num_threads(4)
  
      call system_clock ( clck_counts_beg, clck_rate )
      call calc_val()
      call system_clock ( clck_counts_end, clck_rate )
      write (*, '(&quot;Time = &quot;,f6.3,&quot; seconds.&quot;)')  (clck_counts_end - clck_counts_beg) / real(clck_rate)
  end program main
  </code></pre>
  <p>There should be ways to reduce the amount of allocations (Julia reports 32-45% gc time!) but for now I am too novice to see them, so any comments and tipps are welcome.</p>
  <h2>Edit:</h2>
  <p>Adding <code>@views</code> and correct broadcasting to the while loop improved the Julia speed considerably (as expected, I guess) and hence beats the Matlab loop now. With 4 threads the code now takes only 1.97secs. Specifically,</p>
  <pre><code>    @inbounds for i in 1:len
          c[1,:,i] = @views @. a_grid[i] - a_grid/R .+ 20.0
          c[2,:,i] = @views @. a_grid[i] - a_grid/R
      end
  
      u = @. c^ζ * ζ^(-1)
      @. u[c&lt;=0] = typemin(Float64)
  
      while test&gt;tol &amp;&amp; outeriter&lt;20000
          outeriter += 1
          V_last = deepcopy(V_mat)
  
          @inbounds Threads.@threads for i in 1:len # loop over grid points
              V_mat[1,i], next_mat[1,i] = @views findmax( @. u[1,:,i] + β*V_last[2,:])
              V_mat[2,i], next_mat[2,i] = @views findmax( @. u[2,:,i] + β*V_last[1,:])
          end
          test = @views maximum( @. abs(V_mat - V_last)[!isnan( V_mat - V_last )])
      end
  </code></pre>