# OpenMP stochastic simulation

**URL:** <https://fortran-lang.discourse.group/t/openmp-stochastic-simulation/8606>\
**Category:** Uncategorized\
**Created:** [September 16, 2024, 4:20pm UTC](https://fortran-lang.discourse.group/t/openmp-stochastic-simulation/8606 "2024-09-16T16:20:00Z")\
**Posts on this page:** 12\
**Page:** 2

<div class="post-metadata">

**Author:** ![hkvzjal](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/hkvzjal/32/3055_2.png) [@hkvzjal](https://fortran-lang.discourse.group/u/hkvzjal)\
**Post date:** [September 17, 2024, 10:08pm UTC](https://fortran-lang.discourse.group/t/openmp-stochastic-simulation/8606/21 "2024-09-17T22:08:21Z")

</div>

Looking at your initialization of arrays, it looks like you might be hitting a “first touch” policy problem which can possibly hinder performance. You might try to initialize your arrays such that each thread initializes a subsection… No guarantees thought but it might be worth taking a look at that.

---

<div class="post-metadata">

**Author:** ![urbanjost](https://avatars.discourse-cdn.com/v4/letter/u/0ea827/32.png) [@urbanjost](https://fortran-lang.discourse.group/u/urbanjost)\
**Post date:** [September 17, 2024, 11:13pm UTC](https://fortran-lang.discourse.group/t/openmp-stochastic-simulation/8606/22 "2024-09-17T23:13:57Z")

</div>

As a proof-of-concept you could use the simple random-number generator used here with OpenMP and if the results are encouraging look at some more complex PRNGs?

Using the supplied test cases as-is and then with the intrinsics would help confirm the root cause.  
[https://people.sc.fsu.edu/~jburkardt/f\_src/random\_openmp/random\_openmp.f90](https://people.sc.fsu.edu/~jburkardt/f_src/random_openmp/random_openmp.f90)

Are there MSWIndows tools to show the mutex calls? That might be useful. Even something like a MSWindows version of strace? Not a MSWindows user, so not sure what it has.

---

<div class="post-metadata">

**Author:** ![martin](https://avatars.discourse-cdn.com/v4/letter/m/73ab20/32.png) [@martin](https://fortran-lang.discourse.group/u/martin)\
**Post date:** [September 18, 2024, 5:12am UTC](https://fortran-lang.discourse.group/t/openmp-stochastic-simulation/8606/23 "2024-09-18T05:12:30Z")

</div>

Sorry, looking at the code within the forum is not that easy. I missed that the allocatable variable `row_sum(:)` is not used. But it is actually declared within the parallel region. So the deallocate is automatically invoked. However, I think for a non-allocated variable it should not have much impact (currently I cannot check in windows, the Linux perf tool does not show much interesting stuff, except that random\_number eats up all the resources if used.)

---

<div class="post-metadata">

**Author:** ![martin](https://avatars.discourse-cdn.com/v4/letter/m/73ab20/32.png) [@martin](https://fortran-lang.discourse.group/u/martin)\
**Post date:** [September 18, 2024, 5:29am UTC](https://fortran-lang.discourse.group/t/openmp-stochastic-simulation/8606/24 "2024-09-18T05:29:07Z")

</div>

BTW, not having 64bit unsigned integer numbers with well defined overflow behaviour makes it difficult to implement well-established random number generators in fortran. There are a couple generators which support split and merge operations, which helps to avoid any syncronisation within a parallel region. gfortran implements the xoshiro256\*\* generator and its “jump” function for splitting (see libgfortran/intrinsics/random.c of in its current source tree). But I do not see how that is actually used by the compiler in conjunction with openmp code.

---

<div class="post-metadata">

**Author:** ![PierU](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/pieru/32/1848_2.png) [@PierU](https://fortran-lang.discourse.group/u/PierU)\
**Post date:** [September 18, 2024, 6:38am UTC](https://fortran-lang.discourse.group/t/openmp-stochastic-simulation/8606/25 "2024-09-18T06:38:00Z")

</div>

The only “first touch” that I can see out of the parallel region is `sim_ind_a(1,:) = init_a`, which is a very limited one. I doubt moving it to the parallel region make a significant difference.

However, this makes me think that most of the code is _maybe_ memory bound, as it has to write a very large amound of data (12GB estimated with the parameters as posted). If so, most of code being memory bound and the call to `RANDOM_NUMBER()` being serialized by a mutex, there’s no chance to get a significant speed-up with multithreading.

---

<div class="post-metadata">

**Author:** ![aledinola](https://avatars.discourse-cdn.com/v4/letter/a/3ec8ea/32.png) [@aledinola](https://fortran-lang.discourse.group/u/aledinola)\
**Post date:** [September 18, 2024, 10:13am UTC](https://fortran-lang.discourse.group/t/openmp-stochastic-simulation/8606/26 "2024-09-18T10:13:59Z")

</div>

Thanks! I actually thought of replacing the intrinsic random\_number with a simpler generator from the book “Modern Fortran explained 2018” (solution to exercise 6 chapter 9) and test this against the intrinsic.

---

<div class="post-metadata">

**Author:** ![septc](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/septc/32/77_2.png) [@septc](https://fortran-lang.discourse.group/u/septc)\
**Post date:** [September 18, 2024, 1:02pm UTC](https://fortran-lang.discourse.group/t/openmp-stochastic-simulation/8606/27 "2024-09-18T13:02:55Z")

</div>

On my PC (Ryzen5700X(8-core, max-freq set to 3.4GHz) + Ubuntu22), the timing with gfortran-12 -O3 -march=native -fopenmp is like this:

```auto
1-thread : stdout -> 11.6 sec (time command -> 8.50user 4.55system)
2-thread : stdout -> 6.8 sec (time command -> 8.57user 4.58system)
4-thread : stdout -> 4.2 sec (time command -> 8.63user 4.54system)
8-thread : stdout -> 3.0 sec (time command -> 8.72user 4.64system)

```

If I replace `random_number()` with `u = 0.5d0`, the result is:

```auto
1-thread : stdout -> 6.9 sec (time command -> 3.74user 4.57system)
2-thread : stdout -> 4.4 sec (time command -> 3.83user 4.58system)
4-thread : stdout -> 3.1 sec (time command -> 3.81user 4.64system)
8-thread : stdout -> 2.4 sec (time command -> 3.85user 4.84system)

```

For serial runs, my PC is ~ 2 times slower than your laptop (= Core i7 12800H, max 4.8 GHz?, the memory (DDR4 vs DDR5?) may also be different):  
[https://www.cpubenchmark.net/cpu.php?cpu=Intel+Core+i7-12800H&id=4778](https://www.cpubenchmark.net/cpu.php?cpu=Intel+Core+i7-12800H&id=4778)

But with threading the program becomes faster to some extent, although the scaling is not very good (probably because memory bound? (as mentioned above))

These pages may also be related:

> [@Learning coarrays, collective subroutines and other parallel features of Modern Fortran](https://fortran-lang.discourse.group/t/learning-coarrays-collective-subroutines-and-other-parallel-features-of-modern-fortran/1146/18):
>
> The ifort RANDOM\_NUMBER is known to use an exclusive lock in threaded applications, reducing performance. It might not be the best choice if you’re comparing performance.

> **[MSVC mutex is slower than you might expect](https://stoyannk.wordpress.com/2016/04/30/msvc-mutex-is-slower-than-you-might-expect/)**
>
> TLDR; The 2015 MSVC C++ runtime’s std::mutex (and potentially other) implementation is significantly slower than the equivalent code written by hand. The reason is that the runtime is built w…

Though this is not about ifort + Windows, similar things might be happening…? (But it also seems (to me) that the speed-down is too large for the above laptop.)

---

<div class="post-metadata">

**Author:** ![martin](https://avatars.discourse-cdn.com/v4/letter/m/73ab20/32.png) [@martin](https://fortran-lang.discourse.group/u/martin)\
**Post date:** [September 18, 2024, 5:56pm UTC](https://fortran-lang.discourse.group/t/openmp-stochastic-simulation/8606/28 "2024-09-18T17:56:37Z")

</div>

Note that the consumption part is not parallelised and takes a significant amount of runtime. But at least it should not get slower the more threads are used.

---

<div class="post-metadata">

**Author:** ![aledinola](https://avatars.discourse-cdn.com/v4/letter/a/3ec8ea/32.png) [@aledinola](https://fortran-lang.discourse.group/u/aledinola)\
**Post date:** [September 18, 2024, 6:37pm UTC](https://fortran-lang.discourse.group/t/openmp-stochastic-simulation/8606/29 "2024-09-18T18:37:25Z")

</div>

Thanks for pointing this out. Yes, I should move the consumption part inside the parallel region

---

<div class="post-metadata">

**Author:** ![PierU](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/pieru/32/1848_2.png) [@PierU](https://fortran-lang.discourse.group/u/PierU)\
**Post date:** [September 18, 2024, 8:47pm UTC](https://fortran-lang.discourse.group/t/openmp-stochastic-simulation/8606/31 "2024-09-18T20:47:19Z")

</div>

> [@aledinola](#):
>
> In your first experiment, you are also using the ifort random\_number but nonetheless you get a speed up.

Unless I missed something in his post, it seems to me that @septc used only gfortran for tests. It seems that in contrast to ifort, `RANDOM_NUMBER()` has no mutex in gfortran.

---

<div class="post-metadata">

**Author:** ![aledinola](https://avatars.discourse-cdn.com/v4/letter/a/3ec8ea/32.png) [@aledinola](https://fortran-lang.discourse.group/u/aledinola)\
**Post date:** [September 18, 2024, 9:03pm UTC](https://fortran-lang.discourse.group/t/openmp-stochastic-simulation/8606/32 "2024-09-18T21:03:48Z")

</div>

You’re right, I had missed it

---

<div class="post-metadata">

**Author:** ![aledinola](https://avatars.discourse-cdn.com/v4/letter/a/3ec8ea/32.png) [@aledinola](https://fortran-lang.discourse.group/u/aledinola)\
**Post date:** [September 19, 2024, 10:07pm UTC](https://fortran-lang.discourse.group/t/openmp-stochastic-simulation/8606/33 "2024-09-19T22:07:12Z")

</div>

I installed gfortran and I ran a speed comparison between `ifort` and `gfortran`, serial vs parallel code. The results confirm what others have pointed out. Here is a table with running time in seconds

| | Num cores | Time gfortran | Time ifort |
| --- | --- | --- | --- |
| OS | | Windows 11 | Windows 11 |
| version | | GNU Fortran (GCC) 12.3.0 | Version 2021.11.0 |
| flags | | -O3 -march=native -fopenmp  
 | -O3 -Qipo -QxHost -Qopemp |
| serial | 0 | 7.2 | 4.88 |
| parallel | 16 | 1.016 | 50.82 |

I’m quite happy about the speed gain with gfortran (~7 times). Instead, **ifort is faster in serial but the performance on multiple cores is abysmal**. I think for simulations I’ll go with gfortran (given that I’d like to use the intrinsic random number subroutine). As I side note, this small example is yet another instance that shows how gfortran (an open-source compiler) manages to stay competitive

P.S. If someone can suggest some compiler flag I could use with gfortran to further improve speed (esp. in serial)… thanks 🙂

I report below a slightly updated version of the MWE

```fortran
module mymod
    implicit none
    
    private
    public :: linspace, markov_simul, myerror
    
    contains
    
    subroutine myerror(string)

    implicit none
    character(len=*), intent(in) :: string
    write(*,*) "ERROR: ", string
    write(*,*) "Program will terminate.."
    stop

    end subroutine myerror

    function linspace(a, b, num) result(array)
        ! Input parameters
        real(8), intent(in) :: a ! The starting value
        real(8), intent(in) :: b ! The stopping value
        integer, intent(in) :: num ! The number of points

        ! Output array
        real(8) :: array(num)

        ! Local variables
        integer :: i
        real(8) :: step

        ! Error handling for num <= 1
        if (num == 1) then
            array(1) = a
            return
        endif

        ! Calculate the step size
        step = (b - a) / real(num - 1, 8)

        ! Populate the array
        do i = 1, num
            array(i) = a + (i - 1) * step
        enddo
    end function linspace
    
    function simul_iid(z_prob,dbg) result(ind_sim)
    implicit none
    
    !##### DESCRIPTION #############
    !Simulate a random draw from a discretized distribution 
    ! Note: this is fast if the size of z_prob is small
    ! When z_prob is large (say, more than 100 elements), sample_pmf is faster
    
    !##### INPUTS/OUTPUTS #############
    real(8), intent(in) :: z_prob(:) !discretized prob vector, must be >=0 and sum to 1
    integer, intent(in) :: dbg !0-1 flag
    integer :: ind_sim !index for the draw
    
    !##### LOCALS #############
    integer :: n, i
    real(8) :: u ! random number (a scalar)
    real(8) :: check, prob_sum
    
    !##### ROUTINE CODE #############
    n = size(z_prob)
    
    if (dbg==1) then
        if ( any(z_prob<0d0) .or. any(z_prob>1d0) ) then
            error stop "simul_iid: 0<=prob<=1 violated"
        endif
        check = sum(z_prob)
        if (abs(check-1d0)>1d-6) then
            error stop "simul_iid: prob does not sum to one"
        endif
    endif
    
    call RANDOM_NUMBER(u)
    !u = 0.5d0

    do i = 1,n-1
        !Cumulative sum p1+p2+..+pi
        prob_sum = sum(z_prob(1:i))
        if (u <= prob_sum ) then
            ind_sim = i
            return
        endif !END IF
    enddo !end for

    ! Else, choose the last value
    ind_sim = n
    
    end function simul_iid

    subroutine markov_simul(ysim,Ptran,Tsim,y1,dbg_in)
    implicit none
    
    !##### DESCRIPTION #############
    ! Routine for simulating a Markov chain with time-independent
    ! transition matrix (standard case). To increase speed, set
    ! the debug flag to zero.
    ! This subroutine calls the function simulate_iid.
    
     !##### INPUTS/OUTPUTS #############
    real(8), intent(in) :: Ptran(:,:) ! Transition matrix, dim: (ns,ns)
    integer, intent(in) :: Tsim ! Lenght of the time series to simulate
    integer, intent(in) :: y1 ! Initial condition 
    integer, intent(out) :: ysim(:) ! Simulated time series s.t. y_sim(1)=y1, dim: (Tsim)
    integer, intent(in),optional :: dbg_in ! 0-1 debug flag
    
    !##### LOCALS #############
    integer :: n, i, ns, t, dbg
    real(8) :: check, prob_sum
    
    !##### ROUTINE CODE #############
    if (present(dbg_in)) then
        dbg = dbg_in
    else
        dbg = 0
    endif
    
    ns = size(Ptran,dim=1)
    
    ! Check if inputs are correct (only if debug flag = 1)
    if (dbg==1) then
        if (size(Ptran,1)/=size(Ptran,2)) then
            call myerror("markov_simul: Ptran is not square")
        endif
        if ( any(Ptran<0d0) .or. any(Ptran>1d0) ) then
            call myerror("markov_simul: 0<=prob<=1 violated")
        endif
        if (Tsim < 2) then
            call myerror("markov_simul: Tsim must be at least 2")
        endif
        if ( y1<1 .or. y1>ns) then
            call myerror("markov_simul: initial condition is out of bounds")
        endif
        if ( size(ysim)/=Tsim ) then
            call myerror("markov_simul: ysim has wrong size")
        endif
    endif
    
    ! Simulate series of length Tsim
    ysim(1) = y1
    
    do t = 1,Tsim-1
        ysim(t+1) = simul_iid(Ptran(ysim(t),:),dbg)    
    enddo
    
    
end subroutine markov_simul

end module mymod
!================================================================================================
program main
    ! THIS CODE USES OPENMP
    use mymod, only: linspace, markov_simul 
    use omp_lib 

    implicit none
    
    ! ------------------ Declarations ------------------
    integer :: par_fortran
    integer, parameter :: nthreads = 16
    integer, parameter :: N_sim = 14000 
    integer, parameter :: T_sim = 30000
    integer, parameter :: n_a=500, n_z=2
    real(8) :: t1,t2
    real(8) :: Pi(2,2), z_grid(2), b, grid_max, r, w
    integer, allocatable :: pol_ap_ind(:,:), sim_ind_a(:,:), sim_ind_z(:,:)
    real(8), allocatable :: a_grid(:), sim_a(:,:), sim_c(:,:)
    ! Variables for the simulation
    integer :: hh_c, t, init_a, init_z
    ! ------------------ End declarations ------------------
    
    write(*,*) "Starting program..."
    
    write(*,*) "max # of threads = ",omp_get_max_threads()
    CALL OMP_SET_NUM_THREADS(nthreads)
    !$omp parallel if (par_fortran==1)
    write(*,*) "hello from thread ",omp_get_thread_num() !," of ",nthreads
    !$omp end parallel 
    
    write(*,*) "Type 1 to run in parallel, 0 to run in serial"
    read(*,*) par_fortran

    ! ---------------------- Set economic parameter values ----------------------!
    Pi(1,:) = [0.6d0, 0.4d0]
    Pi(2,:) = [.05d0, 0.95d0]
    z_grid = [0.5d0, 1.0d0]
    r = 0.04d0
    w = 1.0d0
    b = 0.0d0
    grid_max = 4.0d0
    
    ! Set grid for assets
    allocate(a_grid(n_a))
    a_grid = linspace(b,grid_max,n_a)
    
    ! Set policy function for assets (index) 
    allocate(pol_ap_ind(n_a,n_z))
    pol_ap_ind(:,1) = n_a/5 ! Low shock
    pol_ap_ind(:,2) = n_a/2 ! High shock
    
    write(*,*) "============================================================="
    write(*,*) "Starting SIMULATION..."
    write(*,*) "============================================================="
    
    t1 = omp_get_wtime()
    
    allocate(sim_a(T_sim+1,N_sim),sim_c(T_sim,N_sim),sim_ind_a(T_sim+1,N_sim), sim_ind_z(T_sim,N_sim))
    
    ! Initial conditions for state variables "a" and "z"
    init_a = 1
    init_z = n_z/2
    ! Recall: first index is time, second index is individual
    sim_ind_a(1,:) = init_a
    
    ! Loop over time periods to simulate assets
    !t1 = omp_get_wtime()
    
    ! Outer loop over individuals
    !$omp parallel if (par_fortran==1) default(shared) private(hh_c,t)
    !$omp do schedule(dynamic)
    do hh_c=1,N_sim
        !write(*,*) "Doing hh_c = ", hh_c
        ! Simulate a sequence of shocks (length Tsim) for individual hh_c
        call markov_simul(sim_ind_z(:,hh_c),Pi,T_sim,init_z)
        ! Inner loop over time periods
        do t=2,T_sim+1       
            ! Update assets
            sim_ind_a(t,hh_c) = pol_ap_ind(sim_ind_a(t-1,hh_c),sim_ind_z(t-1,hh_c))
            sim_a(t,hh_c) = a_grid(sim_ind_a(t,hh_c))
        enddo  
        ! Update consumption
        do t=1,T_sim
            sim_c(t,hh_c) = (1d0+r)*sim_a(t,hh_c)+w*z_grid(sim_ind_z(t,hh_c))-sim_a(t+1,hh_c)
        enddo  
    enddo
    !$omp enddo 
    !$omp end parallel

    t2 = omp_get_wtime()
    
    write(*,*) "============================================================="
    write(*,*) 'SIMULATION runs for ',t2-t1,' seconds.'
    write(*,*) "============================================================="
    
    write(*,*) "============================================================="
    write(*,*) 'MODEL MOMENTS'
    write(*,'(X,A,F10.4)') "Average assets, simulation: ", sum(sim_a)/size(sim_a)
    write(*,'(X,A,F10.4)') "Average consumption, simulation: ", sum(sim_c)/size(sim_c)
    write(*,*) 'EXAMPLE SIMULATED SERIES FOR 10 INDIVIDUALS at t=1000'
    write(*,'(X,A)') "hh_c, sim_a, sim_c"
    do hh_c=1,10
        write(*,'(I5,10F10.4)') hh_c, sim_a(1000,hh_c), sim_c(1000,hh_c)
    enddo
    
    deallocate(a_grid,pol_ap_ind,sim_a,sim_c,sim_ind_z)
    
    write(*,*) "Program terminated OK!" 
    
end program main

```

[Previous page](https://fortran-lang.discourse.group/t/openmp-stochastic-simulation/8606.md?page=1)