Help with openMP implementation of Fortran Code

Seitz__Sebastian · ‎09-10-2018

Dear all,

I am new to programming with Fortran's openMP application (using Intel Compiler version 18). I recently have run into a problem to which I could not find an answer (after 1 month of digging around). The issue is (see below for the code) that I would like to solve a function where I assign some limiting conditions in the main body of the code, which also should be used in the function. For example, I specify that my function should allow for a current value A_cur=a0(i,j) where i and j are my do-loop running variables.

The issue is that in my main body of code, this assignment works fine. But as soon as I "call" my function f(x) where x not equal A_cur (but A_cur is a constant used in the function), my code reads A_cur=0 instead of A_cur=a0(i,j).

Hence, my question is how to solve this issue without defining the function f(x,a_in) where I set, a_in=A_cur? I would be really grateful for some insights on this issue as well as some recommendations.

This code is supposed to illustrate my issue with a very simple assignment (which however captures the issue). In this code, I specify that a_com=i, whereby i is the running variable in my OMP loop. In the main code it correctly takes the assigned value but it does not do so in the function utility. (Note that this is a simplified version of the assignment problem I mention above where the assignment issue between main body of the code and the function does not affect the outcome, while in my example above it would).

---------------------------------------------------------------------------------------------------------------------------------

module constants

integer, parameter:: i4b=selected_int_kind(15)
integer, parameter:: ndp=selected_real_kind(15)

real(ndp), parameter:: zero=0.0d0
real(ndp), parameter:: one=1.0d0
real(ndp), parameter:: two=2.0d0

integer(i4b):: t_in, a_com
real(ndp):: y_cur,a_cur

contains

double precision function utility(in)
real(ndp),intent(in):: in
real(ndp):: v1, cons


if((t_in)<two)then
cons=1+in
write(*,*) 't_in<2'
else
cons=2+in
write(*,*) 't_in>2'
endif

write(*,*) 'a_com', a_com

v1=sqrt(cons)
utility=v1

end function utility

end module constants

program test_openMP

use constants

!# Define some variables
integer(i4b):: i,t,j
integer(i4b), parameter:: i_max=10, t_max=3, j_max=3
real(ndp):: result(0:i_max,0:t_max,0:j_max), a_pr(0:i_max,0:t_max), &
income(0:j_max)
real(ndp):: a_in,a_max



a_max=100

do i=0,i_max
do t=0,t_max
a_pr(i,t)=(a_max)/((1+0.05*t)**i)
enddo
enddo

do j=0,j_max
income(j)=(1.07)**(1+j)
enddo


do t=0,t_max
t_in=t
write(*,*)'Run for t', t, 'of', t_max
!$OMP PARALLEL DO PRIVATE (a_in,y_cur,a_cur,a_com)
do i=0,i_max
write(*,*)'Run for i', i, 'of', i_max
a_in=a_pr(i,t)
a_com=i
do j=0,j_max
y_cur=income(j)
a_cur=a_in+y_cur
result(i,t,j)=utility(a_cur)
enddo
enddo
!$OMP END PARALLEL DO
!a_com=o
enddo


do i=0,i_max
t=2
j=1
write(*,*) 'i', '|', 'a_pr','|', 'Result','|', 'Utility'
write(*,*)'###############################################'
write(*,*) i,'|', a_pr(i,t),'|', result(i,t,j),'|', utility(a_pr(i,t)+income(j))
enddo

end program test_openMP

--------------------------------------------------------------------------------------------------------------------

Best,

Sebastian

John_Campbell · ‎09-18-2018

private variable a_com is different to the module variable a_com in utility, although it is not used. I would also recommend that array result is shared.

I changed your test example by : removing writes in Utility, declaring shared variables and including a record of thread in use.

The changed example produces similar results for both single thread and multi-thread, although I am not sure if it addresses your problem. It may give you some ideas ?

module constants

 integer, parameter:: i4b=selected_int_kind(8)
 integer, parameter:: ndp=selected_real_kind(15)

 real(ndp), parameter:: zero=0.0d0
 real(ndp), parameter:: one=1.0d0
 real(ndp), parameter:: two=2.0d0

 integer(i4b):: t_in, a_com
 real(ndp):: y_cur,a_cur

    contains
     
     double precision function utility(in)
         real(ndp),intent(in):: in
         real(ndp):: v1, cons
         
         
         if((t_in)<two)then
             cons=1+in
!             write(*,*) 't_in<2'
         else
             cons=2+in
!             write(*,*) 't_in>2'
         endif
         
!         write(*,*) 'a_com', a_com
         
         v1=sqrt(cons)
         utility=v1
         
     end function utility

end module constants

program test_openMP

    use constants
    USE omp_lib
     
     !# Define some variables
     integer(i4b):: i,t,j, id
     integer(i4b), parameter:: i_max=10, t_max=3, j_max=3
     real(ndp):: result(0:i_max,0:t_max,0:j_max), a_pr(0:i_max,0:t_max), &
                 income(0:j_max)
     real(ndp):: a_in,a_max
     integer(i4b):: last_thread_used(0:i_max)
     
     a_max=100
     
     do i=0,i_max
         do t=0,t_max
             a_pr(i,t)=(a_max)/((1+0.05*t)**i)
         enddo 
     enddo
     
     do j=0,j_max
         income(j)=(1.07)**(1+j)
     enddo
     
     
     do t=0,t_max
         t_in=t
         write(*,*)'Run for t', t, 'of', t_max
      !$OMP PARALLEL DO PRIVATE (a_in,y_cur,a_cur,a_com, id)   &
      !$OMP             SHARED  (a_pr, income, result, last_thread_used)
         do i=0,i_max
             id = -1
      !$     id = omp_get_thread_num()
             last_thread_used(i) = id
             write(*,*)'Run for i', i, 'of', i_max, id
             a_in=a_pr(i,t)
             a_com=i
             do j=0,j_max
                 y_cur=income(j)
                 a_cur=a_in+y_cur
                 result(i,t,j)=utility(a_cur)
             enddo
         enddo
      !$OMP END PARALLEL DO
      !a_com=o
     enddo
     
     
     do i=0,i_max
       t=2
       j=1
       write(*,*) 'i', '|', 'a_pr','|', 'Result','|', 'Utility'
       write(*,*)'###############################################'
       write(*,*)  i,'|', a_pr(i,t),'|', result(i,t,j),'|', utility(a_pr(i,t)+income(j)),last_thread_used(i)
     enddo

end program test_openMP

Greg_T_ · ‎09-18-2018

If additional information about OpenMP would be helpful in getting started, there is a very good set of notes and examples available at LLNL:

https://computing.llnl.gov/tutorials/openMP/

And the "Introduction to OpenMP" videos by Tim Mattson from Intel available on YouTube is excellent. This should be the URL to that collection of videos:

https://www.youtube.com/playlist?list=PLLX-Q6B8xqZ8n8bwjGdzBJ25X2utwnoEG

Regards, Greg

Seitz__Sebastian · ‎09-26-2018

Dear Greg, dear John,

thank you for your helpful comments. Unfortunately, my program is still not running. The above example now does, but the actual code is a little bit more complex than the easy case above. In general, it helped me to know that the variables within the module are regarded as different variables when using openMP as the issue seems to be be partially alleviated if I write the (private) variables as inputs into my subroutine.

That being said, even if this subroutine contains a function then, the program does not seem to write correctly into this function. As annoying as this is at the moment, I have decided to move on (I need to work on the open MPI implementation as well). Do there exist similar great sources as the ones you have mentioned, Greg?

Best,

Sebastian

jimdempseyatthecove · ‎09-26-2018

>> it helped me to know that the variables within the module are regarded as different variables when using openMP

module variables are SHARED unless stated differently (OpenMP clause or threadprivate)

Inside the function utility (in post #2), t_in and a_com are the module (shared) copy (also the main thread's values).
However, in the !$omp loop, t_in and a_com are the private copies, with exception that for the main thread they are the same as the module variables.

Jim Dempsey

Seitz__Sebastian · ‎09-26-2018

Dear Jim,

in my main code I actually used to define them as:

integer(i4b), private:: t_in, a_com

I thought this would then make this variable threadprivate in my module. Am I mistaken about that?

If there is another way of making them threadprivate (within the module where I define my function), I would be very grateful to know this way. The only other way I have found is to explicitly write these variables into my subroutine/ function and declare them as private in my !$OMP statement/call.

Best,

Sebastian

jimdempseyatthecove · ‎09-27-2018

>>integer(i4b), private:: t_in, a_com

This does not make the variable threadprivate. What it does is makes the variable hidden from view outside the scope of the module.

OpenMP has compiler directives that need not be associated with a parallel region.

module
...
!$omp threadprivate(t_in, a_com)
integer(i4b):: t_in, a_com
...
contains
...
! code follows

Jim Dempsey