Intel® Fortran Compiler
Build applications that can scale for the future with optimized code designed for Intel® Xeon® and compatible processors.
29246 Discussions

Thread data sharing using IMSL and OpenMP

Christian_B_4
Beginner
‎05-20-2013 08:24 PM
625 Views

Hello,

First of all I should say that I am quite new to Fortran and OpenMP, so I apologize in advance for the lack of proper vocabulary and/or knowledge....

I have to run a series of minimizations for the same function but with different argument values. I achieve this in a do-loop and I am trying to parallelize that process with OpenMP. The minimization is done with the DBCPOL routine of the ISML library. I manage to do that, but when debugging the code with the /debug:parallel option to analize the Thread Data Sharing (concerned about the claim that IMSL is thread-safe) there are several events which I cannot either diagnose or dismiss.A minimal working sample code showing the problem is found below (this is not the actual problem, but it is an attempt to understand where the problem is, setting aside the complexity of the real functions and providing a MWE for you to try):

Mail.f90

[fortran] program AAA
use XX
INCLUDE 'link_fnl_static.h' !(To use the static IMSL library form of the libraries)
!DEC$ OBJCOMMENT LIB:'libiomp5md.lib'
call MPublic()
pause
end program AAA [/fortran]

Modules.f90

[fortran] module data_type
    implicit none
    integer, parameter :: IB=4, RP=8
end module data_type

module XX
    use data_type
    implicit none
    real(kind=RP),private, dimension(2) :: ysub
    private :: obj_func, minimum
    contains
    function minimum(yin)
        implicit none
        real(kind=RP), dimension(2) ,intent(in):: yin
        integer(kind=IB), parameter :: Np=2, IBTYPE=0
        integer(kind=IB)            :: max_eval
        real(kind=RP), dimension(2) :: minimum
        real(kind=RP), dimension(2) :: x_opt, x_min, x_max
        real(kind=RP), dimension(2), parameter :: x0 = (/ 0._rp , 0._rp /)
        real(kind=RP)               :: tol, obj_opt
        ysub = yin
        tol      = 1.d-6
        max_eval = 8000
        x_max    = (/ 800._rp , 800._rp /)
        x_min    = (/ -600._rp , -600._rp /)
        call dbcpol(obj_func, Np, x0, IBTYPE, x_min, x_max, tol, max_eval, x_opt, obj_opt)
        minimum=x_opt
    end function minimum
    subroutine obj_func(Np, x, val)
        implicit none
        integer(kind=IB), intent(in)             :: Np
        real(kind=RP), dimension(Np), intent(in) :: x
        real(kind=RP) , intent(out)              :: val
        real(kind=RP), dimension(2) :: yin
        yin=ysub
        val=(x(1)-yin(1))**2+(x(2)-yin(2))**2
    end subroutine obj_func
    subroutine MPublic()
    use omp_lib
    implicit none
    integer :: start = 1
    integer ,parameter ::  Nvar = 200
    integer :: i ,k
    real(kind=RP) ,dimension(Nvar,2) :: Minf
    real(kind=RP), dimension(2) :: y,yt
    integerNTHREADS, TID

    !$omp parallel do private(y,TID,ysub)
    do i=start,Nvar
    !$ TID = OMP_GET_THREAD_NUM()  !     Obtain thread id   
    y=(/ i , i /)
    Minf(i,:)=minimum(y)
    write(*,fmt="(I3,xx,2(f9.4,x),I3)")  i, (y/Minf(i,:)),TID
    enddo
    !$omp end parallel do

    end subroutine
end module XX [/fortran]

I get a sharing event with the VAR variable, used in the objective function definition required by DBCPOL. Also you should note that I need the value of the y variable there, even when it is not in its arguments.

Without the Data sharing diagnose, I should get as an answer a set of 1.0000 1.0000. It can be seen that some of the values, depending on the run, are not correct, I guess due to the "sharing" problem.

I am not sure whether it is a problem with the implementation, or with the library, but I would really appreciate any insight you can provide about it.

I am using the Intel Fortran compiler 11.1, IMSL FNL 7.0 for Intel(R) Fortran Compiler 11.1 - IA32 and Visual Studio 2008.

In an attempt to fix this problem, I have tried another code for 2 threads, where I use two modules, each with duplicate functions (with different names) in order to avoid the sharing. Then, instead of a PARALLEL DO, I use just PARALLEL and with an IF clause, I manually distribute the work for each thread. The same problem was observed there regarding the data sharing events, but in this case the results seem to be ok. For what I have read, not all of those events imply a potential error caused by different threads "working" on the same variables, but as a newbie, I am quite confused about it.

Main.f90

[fortran] program AAA
INCLUDE 'link_fnl_static.h' !(To use the static library form of the libraries)
!DEC$ OBJCOMMENT LIB:'libiomp5md.lib'
use MINIMIZE
call MPublic()
pause
end program AAA [/fortran]

Modules.f90

[fortran] module data_type
    implicit none
    integer, parameter :: IB=4, RP=8
end module data_type

module XX1
    use data_type
    implicit none
    real(kind=RP), dimension(2),private :: ysub1
    private :: obj_func
    contains
    function minimum(yin1)
        implicit none
        real(kind=RP), dimension(2) ,intent(in):: yin1
        integer(kind=IB), parameter :: Np=2, IBTYPE=0
        integer(kind=IB)            :: max_eval
        real(kind=RP), dimension(2) :: minimum,x_opt, x_min, x_max
        real(kind=RP), dimension(2), parameter :: x0 = (/ 0._rp , 0._rp /)
        real(kind=RP)               :: tol, obj_opt
        ysub1 = yin1
        tol      = 1.d-6
        max_eval = 10000
        x_max    = (/ 1000._rp , 1000._rp /)
        x_min    = (/ -500._rp , -500._rp /)
        call dbcpol(obj_func, Np, x0, IBTYPE, x_min, x_max, tol, max_eval, x_opt, obj_opt)
        minimum=x_opt
    end function minimum

     subroutine obj_func(Np, x, val1)
         implicit none
        integer(kind=IB), intent(in)             :: Np
        real(kind=RP), dimension(Np), intent(in) :: x
        real(kind=RP) , intent(out)              :: val1
        real(kind=RP), dimension(2) :: yin1
        yin1=ysub1
        val1=(x(1)-yin1(1))**2+(x(2)-yin1(2))**2
    end subroutine obj_func
end module XX1

module XX2
    use data_type
    implicit none
    real(kind=RP), dimension(2),private :: ysub2
    private :: obj_func2
    contains
    function minimum2(yin2)
        implicit none
        real(kind=RP), dimension(2) ,intent(in):: yin2
        integer(kind=IB), parameter :: Np=2, IBTYPE=0
        integer(kind=IB)            :: max_eval
        real(kind=RP), dimension(2) :: minimum2,x_opt2, x_min, x_max
        real(kind=RP), dimension(2), parameter :: x0 = (/ 0._rp , 0._rp /)
        real(kind=RP)               :: tol, obj_opt2
        ysub2 = yin2
        tol      = 1.d-6
        max_eval = 10000
        x_max    = (/ 1000._rp , 1000._rp /)
        x_min    = (/ -500._rp , -500._rp /)
        call dbcpol(obj_func2, Np, x0, IBTYPE, x_min, x_max, tol, max_eval, x_opt2, obj_opt2)
        minimum2=x_opt2
    end function minimum2

    subroutine obj_func2(Np, x2, val2)
         implicit none
        integer(kind=IB), intent(in)             :: Np
        real(kind=RP), dimension(Np), intent(in) :: x2
        real(kind=RP) , intent(out)              :: val2
        real(kind=RP), dimension(2) :: yin2
        yin2=ysub2
        val2=(x2(1)-yin2(1))**2+(x2(2)-yin2(2))**2
    end subroutine obj_func2
end module XX2

module MINIMIZE
    use data_type
    use XX1
    use XX2
    implicit none

    contains

    subroutine MPublic()
        use omp_lib
        use XX1
        use XX2
        implicit none
        integer :: start = 1
        integer ,parameter ::  Nvar = 400
        integer :: i
        real(kind=RP) ,dimension(Nvar,2) :: Minf
        real(kind=RP), dimension(2) :: y
        integer NTHREADS, TID
        NTHREADS=2
        !$omp parallel  default(none) private(y,TID,i) shared(start,Minf,NTHREADS)
        !$ TID = OMP_GET_THREAD_NUM()  !     Obtain thread id
          
        IF (TID .EQ. 0) THEN     !     Only master thread does this
            do i=start,(Nvar*(1+TID)/NTHREADS)
                y=(/ i , i /)
                Minf(i,:)=minimum(y)
                write(*,fmt="(I3,xx,2(f9.4,x),I3)")  i, y/Minf(i,:),TID
            enddo
        else IF (TID .EQ. 1) THEN     !     
            do i=(Nvar*(TID)/NTHREADS)+1,(Nvar*(1+TID)/NTHREADS)
                y=(/ i , i /)
                Minf(i,:)=minimum2(y)
                write(*,fmt="(I3,xx,2(f9.4,x),I3)")  i, y/Minf(i,:),TID
            enddo
        end if
        !$omp end parallel
    end subroutine
end module MINIMIZE [/fortran]

Thank you very much for having taking, at least, the time to read this gigantic question. Any other comments or suggestions are welcome as well.

Best regards,

Christian

0 Kudos

Link Copied

1 Reply
Steven_L_Intel1
Employee
‎05-21-2013 06:37 AM
625 Views

Since you are using IMSL 7.0, which you did not obtain from Intel,. I recommend that you seek advice in Rogue Wave's IMSL forums at http://forums.roguewave.com/forumdisplay.php?110-IMSL-Numerical-Libraries In addition to being a place to get support for IMSL 7, you also get the benefit of direct interaction with IMSL developers and support engineers.

I will observe that many IMSL routines aren't really thread-safe in the sense that you can't call the same IMSL routine from multiple threads and always get correct behavior, though a given IMSL routine will function properly in a multithread environment. Many IMSL routines use static storage and if you have simultaneous calls they will interfere with each other. Particularly problematic are the IMSL functions that call back to user-passed routines.

0 Kudos
Copy link
Reply

Community support is provided Monday to Friday. Other contact methods are available here.

Intel does not verify all solutions, including but not limited to any file transfers that may appear in this community. Accordingly, Intel disclaims all express and implied warranties, including without limitation, the implied warranties of merchantability, fitness for a particular purpose, and non-infringement, as well as any warranty arising from course of performance, course of dealing, or usage in trade.

For more complete information about compiler optimizations, see our Optimization Notice.