# OpenMP Parallel Loop **URL:** \ **Category:** Help\ **Created:** [July 28, 2023, 9:13am UTC](https://fortran-lang.discourse.group/t/openmp-parallel-loop/6244 "2023-07-28T09:13:02Z")\ **Posts on this page:** 20\ **Page:** 1
**Author:** ![phangelder](https://avatars.discourse-cdn.com/v4/letter/p/35a633/32.png) [@phangelder](https://fortran-lang.discourse.group/u/phangelder)\ **Post date:** [July 28, 2023, 9:13am UTC](https://fortran-lang.discourse.group/t/openmp-parallel-loop/6244/1 "2023-07-28T09:13:02Z")
Would appreciate how to use OpenMP to parallel the following code ```auto maxvab=100000 ; thispos=0 ! jwork(random) contains 1 or 0 in random do ivab=1,maxvab if(jwork(ivab) == 0)cycle jWork(ivab)=thisPos thisPos=thisPos+1 end do ``` ---
**Author:** ![everythingfunctional](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/everythingfunctional/32/176_2.png) [@everythingfunctional](https://fortran-lang.discourse.group/u/everythingfunctional)\ **Post date:** [July 28, 2023, 12:21pm UTC](https://fortran-lang.discourse.group/t/openmp-parallel-loop/6244/3 "2023-07-28T12:21:56Z")
> [@phangelder](#): > > ```auto > jWork(ivab)=thisPos > thisPos=thisPos+1 > > ``` This pair of lines mean that you cannot run the loop in parallel without getting different answers. If you do the iteration where `ivab == maxvab` first, `jwork(maxvab)` will get a different value than if you did it last. This is called a race condition. ---
**Author:** ![phangelder](https://avatars.discourse-cdn.com/v4/letter/p/35a633/32.png) [@phangelder](https://fortran-lang.discourse.group/u/phangelder)\ **Post date:** [July 28, 2023, 1:29pm UTC](https://fortran-lang.discourse.group/t/openmp-parallel-loop/6244/5 "2023-07-28T13:29:02Z")
you are right. What would you do to solve this problem which appears all the time in sparse solver? ---
**Author:** ![everythingfunctional](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/everythingfunctional/32/176_2.png) [@everythingfunctional](https://fortran-lang.discourse.group/u/everythingfunctional)\ **Post date:** [July 28, 2023, 1:42pm UTC](https://fortran-lang.discourse.group/t/openmp-parallel-loop/6244/6 "2023-07-28T13:42:52Z")
This looks like it needs a parallel `scan` operation. I.e. ```auto pure function scan(arr, op, init) result(scanned) type(T), intent(in) :: arr(:) type(U), intent(in) :: init interface pure function op(lhs, rhs) type(T), intent(in) :: lhs type(U), intent(in) :: rhs type(U) :: op end function end interface type(U) :: scanned(size(arr)) end function ``` Then you can do ```auto positions = scan(jWork, increment_if_positive, 0) where (jWork > 0) jWork = positions ``` Writing that `scan` to work in parallel is hard™. The conditional assignment in parallel is easy. The compiler can probably do it for you. ---
**Author:** ![phangelder](https://avatars.discourse-cdn.com/v4/letter/p/35a633/32.png) [@phangelder](https://fortran-lang.discourse.group/u/phangelder)\ **Post date:** [July 28, 2023, 1:46pm UTC](https://fortran-lang.discourse.group/t/openmp-parallel-loop/6244/7 "2023-07-28T13:46:26Z")
i think a sort routine for jwork after the loop will work but would it perform better than the original serial code? ---
**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:** [July 28, 2023, 1:47pm UTC](https://fortran-lang.discourse.group/t/openmp-parallel-loop/6244/8 "2023-07-28T13:47:44Z")
There are some workarounds, but with a much more complex code. And given the pretty low workload here (100000 iterations with that few instruction is a low workload), it will very likely be even slower than the simple serial version (at least on usual hardware). ---
**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:** [July 28, 2023, 2:38pm UTC](https://fortran-lang.discourse.group/t/openmp-parallel-loop/6244/9 "2023-07-28T14:38:48Z")
That said, OpenMP 5 has introduced the `scan` directive. I’m not sure how it can be used here, but it could be worth investigating. ---
**Author:** ![JohnCampbell](https://avatars.discourse-cdn.com/v4/letter/j/5daacb/32.png) [@JohnCampbell](https://fortran-lang.discourse.group/u/JohnCampbell)\ **Post date:** [July 29, 2023, 4:27am UTC](https://fortran-lang.discourse.group/t/openmp-parallel-loop/6244/10 "2023-07-29T04:27:21Z")
> [@phangelder](#): > > Would appreciate how to use OpenMP to parallel the following code Typically this DO loop would be part of a larger OMP computation. This DO should be allocated to a single thread, while other “jwork” vectors are allocated to other threads. Typically there are many vectors, which can be allocated dynamically to the next available thread. You need to monitor the status of each jwork vector, as to if they are: - not yet processed, - being processed, or - have completed processing. Often they need to be allocated sequentially, but during their processing, the other vectors need to be monitored. Processing can not be completed until earlier vectors are completed, which can require waiting. This can be as simple as monitoring “last\_vector\_processed”, which can be updated as an atomic construct. This can still be very efficient, as the vector size is typically much greater than the number of active threads. ---
**Author:** ![phangelder](https://avatars.discourse-cdn.com/v4/letter/p/35a633/32.png) [@phangelder](https://fortran-lang.discourse.group/u/phangelder)\ **Post date:** [July 29, 2023, 6:51am UTC](https://fortran-lang.discourse.group/t/openmp-parallel-loop/6244/11 "2023-07-29T06:51:33Z")
Dear John Would be a very neat approach to solve this problem. Can you share an example of the implementation of processed/being processed/completed ? regards ---
**Author:** ![JohnCampbell](https://avatars.discourse-cdn.com/v4/letter/j/5daacb/32.png) [@JohnCampbell](https://fortran-lang.discourse.group/u/JohnCampbell)\ **Post date:** [July 30, 2023, 1:06pm UTC](https://fortran-lang.discourse.group/t/openmp-parallel-loop/6244/12 "2023-07-30T13:06:56Z")
The practical example I have is with an !$OMP skyline direct solver. As equations must be reduced sequentially, I use "last\_vector\_processed”, and monitor this when reducing other equation columns in other threads. Again, as the bandwidth (vector size) is much larger than the number of threads, substantial computation can be done before the threads interact, waiting for "last\_vector\_processed” approaches this equation number. I do use a function “wait\_a\_while” to suspend processing, but monitor how often and for how long this is used. It does not provide a significant delay to the calculation. The following code appears to achieve what is wanted. it must be compiled with -O0 or -O1 to avoid being optimised away ! ```auto ! These routines cause a thread to wait if the equation "Jeq" is not complete ! The delay must be short so as not to delay the restart of this column ! and also not get into a tight loop in this thread ! The call to QueryPerformance_tick > QueryPerformanceCounter is used to ! hopefully prevent this occurring. ! ! the performance of this routine varies with different processors ! subroutine wait_on_jeq ( Jeq, Last_eqn, op_wait ) ! integer*4 Jeq ! equation that is being waited for integer*4 Last_eqn ! last equation to be completed integer*8 op_wait ! count of equations that are waiting ! !z logical NA_wait(*) ! shared wait status for each equation : .false. when complete ! integer*4 iw ! loop count integer*8 tick ! dummy argument to overcome /O1 ! iw = 0 DO if ( Last_eqn >= Jeq ) return call take_a_small_delay (iw, op_wait, tick) !$OMP FLUSH ( Last_eqn ) END DO ! end subroutine wait_on_jeq subroutine take_a_small_delay (iw, op_wait, tick) integer*4 :: iw ! count of delay calls for this equation integer*8 :: op_wait ! count of equations that are waiting integer*8 :: tick ! dummy argument to overcome /O1 ! integer*8, external :: QueryPerformance_tick ! ! The following code order does not work well on i7-4790K ! I shall retest to see if I can reproduce this result ! retest with different compile options -O1 helped if ( iw == 0 ) op_wait = op_wait + 1 ! first call if ( mod(iw,4) == 0 ) tick = QueryPerformance_tick () iw = iw + 1 if ( iw < 0 ) iw = 0 ! end subroutine take_a_small_delay ``` ---
**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:** [August 1, 2023, 1:09pm UTC](https://fortran-lang.discourse.group/t/openmp-parallel-loop/6244/13 "2023-08-01T13:09:14Z")
I’m not sure about what you are doing, but it seems to me that the OpenMP task directive with the depend clause can handle this kind of problem. ---
**Author:** ![phangelder](https://avatars.discourse-cdn.com/v4/letter/p/35a633/32.png) [@phangelder](https://fortran-lang.discourse.group/u/phangelder)\ **Post date:** [August 2, 2023, 12:31am UTC](https://fortran-lang.discourse.group/t/openmp-parallel-loop/6244/14 "2023-08-02T00:31:27Z")
I am solving large scale eigenvalue problems with Arpack Sparse solver on a 20 cores PC. Need to form three matrices in CSR format. Typical size of each matrix is 5millions x 5 millions. Typical outer loop is over a million of points and inner loop is over a million of cells. The inner loop is as listed below which i have tried to parallelize ``` do jelem=1,nelz ! ( size(nelz) >1 million)) kount=0 do nod=1,lnomax ! (size(8)) ! hot spot if(lnods(jelem,nod) == inode)kount=kount+1 end do if(kount == 0)cycle do nod=1,lnomax jnode=lnods(jelem,nod) nf2=(jnode-1)*ndfmax do idof=1,ndfmax ! (size(ndfmax) = 6) jWork(nf2+idof)=jWork(nf2+idof)+1 end do end do end do ``` Thank you in advance for any suggestion. ---
**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:** [August 2, 2023, 12:15pm UTC](https://fortran-lang.discourse.group/t/openmp-parallel-loop/6244/15 "2023-08-02T12:15:39Z")
```auto !$OMP PARALLEL DO PRIVATE(kount,nod,jnode,nf2) REDUCTION(+:jWork) do jelem=1,nelz ! ( size(nelz) >1 million)) kount = count(lnods(jelem,1:lnomax) == inode) if(kount == 0)cycle do nod=1,lnomax jnode=lnods(jelem,nod) nf2=(jnode-1)*ndfmax jWork(nf2+1:nf2+ndfmax) = jWork(nf2+1:nf2+ndfmax) + 1 end do end do ``` In the case `jnode` is guaranteed to be different for different `jelem` values, then the reduction clause can be dropped ---
**Author:** ![phangelder](https://avatars.discourse-cdn.com/v4/letter/p/35a633/32.png) [@phangelder](https://fortran-lang.discourse.group/u/phangelder)\ **Post date:** [August 3, 2023, 12:37am UTC](https://fortran-lang.discourse.group/t/openmp-parallel-loop/6244/16 "2023-08-03T00:37:06Z")
The use of reduction clause will create a copy of JWORK for each thread. The size of JWORK is the same as the number of equation and will cause stack overflow with larger problem. Would it be possible to parallel this section with critical, atomic and/or lock? If so, how and where would be the best location to give the best performance? ---
**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\ **Post date:** [August 3, 2023, 8:36am UTC](https://fortran-lang.discourse.group/t/openmp-parallel-loop/6244/17 "2023-08-03T08:36:46Z")
You can use the Intel Advisor to model threading designs based on the measured serial execution properties. See [Model Threading Design | Intel® Advisor User Guide](https://www.intel.com/content/www/us/en/docs/advisor/user-guide/2023-2/model-threading-designs.html#GUID-78F396B1-BCBC-4697-AA8B-A5F32EDB5F74). The idea is you instrument your code with special annotations: ```fortran use advisor_annotate ... call annotate_site_begin("assembly") do jelem=1,nelz ! ( size(nelz) >1 million)) call annotate_task_begin("count_and_reduce") kount=0 do nod=1,lnomax ! (size(8)) ! hot spot if(lnods(jelem,nod) == inode)kount=kount+1 end do if(kount == 0)cycle do nod=1,lnomax jnode=lnods(jelem,nod) nf2=(jnode-1)*ndfmax do idof=1,ndfmax ! (size(ndfmax) = 6) call annotate_lock_acquire(0) jWork(nf2+idof)=jWork(nf2+idof)+1 call annotate_lock_release(0) end do end do call annotate_task_end end do call annotate_site_end ``` The tool can then help you explore different threading models (for Fortran you’ll typically only look at OpenMP). By the way, I’d recommend reading the following resources from the European Performance Optimization and Productivity (PoP) Center of Excellence: - [Parallelization of indirect reductions on large data structures](https://co-design.pop-coe.eu/patterns/indirect-reductions-large-data-structure.html) - [Indirect reductions with atomic or critical constructs on large data structures](https://co-design.pop-coe.eu/patterns/indirect-reductions-atomic-large-data-structure.html) - [Multidependencies in indirect reductions on large data structures](https://co-design.pop-coe.eu/best-practices/multideps-indirect-reductions-large-data-structure.html) The following picture that I borrowed from the Barcelona Supercomputing Centre [facebook page](https://www.facebook.com/BSCCNS/posts/openmp-architecture-has-recently-been-enhanced-with-the-support-of-ompss-multide/1539916496085543/) illustrates the different parallelization concepts: ![image](https://global.discourse-cdn.com/free1/uploads/fortran_lang/original/2X/e/e52e7d36441710392d34b2a3f3dec17d4f6aac60.jpeg) _Addendum_: here are a few more resources from the experts, - [Using OpenMP Tasking | YouTube](https://www.youtube.com/watch?v=Wx4eQQihP6I), lecture by Christian Terboven and Dirk Schmidl - [How to Get OpenMP Tasks To Do Your Work (OpenMP Webinar) | YouTube](https://www.youtube.com/watch?v=zdNqfG3yBew), lecture by Ruud van der Pas - [EUROfusion OpenMP on GPUs Webinars (1 to 7)](https://indico.euro-fusion.org/event/688/) with Michael Klemm and Christian Terboven; [Webinar 2](https://indico.euro-fusion.org/event/688/attachments/854/2179/OpenMP_Webinar_2_2021-02-16.pdf) (PDF, 2.9 MB) covers tasking including dependencies. Ruud van der Pas and Christian Terboven are also authors of the book [_Using OpenMP - The Next Step: Affinity, Accelerators, Tasking, and SIMD_](https://direct.mit.edu/books/book/4482/Using-OpenMP-The-Next-StepAffinity-Accelerators), The MIT Press, 2017. (Caveat - the book focuses on OpenMP 4.5, but the multidependencies were only introduced in OpenMP 5.0.) ---
**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\ **Post date:** [August 3, 2023, 10:33am UTC](https://fortran-lang.discourse.group/t/openmp-parallel-loop/6244/18 "2023-08-03T10:33:19Z")
> [@PierU](#): > > ```fortran > kount = count(lnods(jelem,1:lnomax) == inode) > if(kount == 0)cycle > do nod=1,lnomax > jnode=lnods(jelem,nod) > nf2=(jnode-1)*ndfmax > jWork(nf2+1:nf2+ndfmax) = jWork(nf2+1:nf2+ndfmax) + 1 > end do > > ``` This could also be modified as: ```fortran if (any(lnods(jelem,1:lnomax) == inode)) then do nod = 1, lnomax ! ... end do end if ``` The discontigous access on the second dimension, since the first dimension of `lnods` is of size `nelzs` (over a million), is not ideal IMO. Hopefully the compiler recognizes it’s a constant offset at least, and then uses the cached result (the row of `lnods(jelem,:)`) in the second loop. ---
**Author:** ![phangelder](https://avatars.discourse-cdn.com/v4/letter/p/35a633/32.png) [@phangelder](https://fortran-lang.discourse.group/u/phangelder)\ **Post date:** [August 3, 2023, 11:05am UTC](https://fortran-lang.discourse.group/t/openmp-parallel-loop/6244/19 "2023-08-03T11:05:02Z")
even with the use of any, i still can’t get my parallel region to work correctly. !$omp parallel & !$omp private (jelem, inode, nod, jnode, nf2, idof) & !$omp shared (nelz, lnods, lnomax, ndfmax, jwork) !$omp do ``` do jelem=1,nelz if (any(lnods(jelem,1:lnomax) == inode)) then do nod=1,lnomax jnode=lnods(jelem,nod) nf2=(jnode-1)*ndfmax do idof=1,ndfmax jWork(nf2+idof)=jWork(nf2+idof)+1 end do end do endif end do ``` !$omp end do !$omp end parallel ---
**Author:** ![phangelder](https://avatars.discourse-cdn.com/v4/letter/p/35a633/32.png) [@phangelder](https://fortran-lang.discourse.group/u/phangelder)\ **Post date:** [August 3, 2023, 11:24am UTC](https://fortran-lang.discourse.group/t/openmp-parallel-loop/6244/20 "2023-08-03T11:24:42Z")
thankyou for the reading recommendation ---
**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\ **Post date:** [August 3, 2023, 12:25pm UTC](https://fortran-lang.discourse.group/t/openmp-parallel-loop/6244/21 "2023-08-03T12:25:31Z")
Try making the `jwork` update atomic. ---
**Author:** ![phangelder](https://avatars.discourse-cdn.com/v4/letter/p/35a633/32.png) [@phangelder](https://fortran-lang.discourse.group/u/phangelder)\ **Post date:** [August 3, 2023, 12:30pm UTC](https://fortran-lang.discourse.group/t/openmp-parallel-loop/6244/22 "2023-08-03T12:30:58Z")
already try and not luck. the issue is several threads executing the if statement in parallel ``` if (any(lnods(jelem,1:lnomax) == inode)) then .... endif ``` [Next page](https://fortran-lang.discourse.group/t/openmp-parallel-loop/6244.md?page=2)