- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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
Link Copied
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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.
![](/skins/images/D2683F18326913BBA0436CB7114DD569/responsive_peak/images/icon_anonymous_message.png)
- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page