0
 subroutine born_par(iconj)    !.false.=0, .true.=1
 use maxwell
 use fourierin
 implicit none
 real*8:: zkx,zky,zkz,zsx,zsy,zsz,zwx,zwy,zwz,ztim,second
 complex*16:: z,zak,zakx,zaky,zakz,zb,zg
 integer, intent(in):: iconj
 integer :: ik,im1,in1,il1,im2,in2,il2,ix,iy,iz

 real*8:: bx(0:mm+2),by(0:nn+2),bz(0:ll+2)
 bx=.0
 by=.0
 bz=.0


im1 = mm+1
in1 = nn+1
il1 = ll+1
im2 = mm+2
in2 = nn+2
il2 = ll+2


do ik = 1,3
print *,'1'
!$omp parallel   private(iy,iz,bx,by) shared(ll,nn,mm,b,im1)
!$omp do
do iz = 1,ll
 do iy = 1,nn
    bx=.0
    by=.0
    bx(1:mm)=real(b(1:mm,iy,iz,ik))       
    by(1:mm)= aimag(b(1:mm,iy,iz,ik))  
    call dst_fwd(im1,bx(0:im1))             
    call dst_fwd(im1,by(0:im1))
    b(1:mm,iy,iz,ik) =cmplx( bx(1:mm),by(1:mm))
  enddo
enddo
!$omp end do
!$omp end parallel
enddo
return
end subroutine born_par

The code snippet above crashes on runtime with no error message when I put OpenMP commands around the the two loops in order to compute Sine Transforms (subroutine dst_fwd ) in parallel. I am using ifort compiler and MKL Intel library for the discrete sine transform.

Does anyone could give me a hint what may cause the program to crash?

The following piece of code defines the module fourierin plus the subroutine dst_fwd

 include 'mkl_dfti.f90'
 include 'mkl_trig_transforms.f90'
 module fourierin
 use mkl_dfti
 use mkl_trig_transforms

 implicit none

  type(DFTI_DESCRIPTOR),pointer:: desc_h1


 contains

subroutine dst_fwd(n,f)
 implicit none
 integer:: n,istr,tt_type ,ipar(128),ir
integer,dimension(2):: istride
real*8:: scale, dpar(3*n/2+2),f(n+1)



f=f(1:n)*(n/2.0d0)


tt_type = MKL_SINE_TRANSFORM

call d_init_trig_transform(n,tt_type,ipar,dpar,ir)
if (ir.ne.0)then
  print *,'Error init fwd sine transform st'
stop
endif

call d_commit_trig_transform(f,desc_h1,ipar,dpar,ir)
if (ir.ne.0)then
  print *, 'Error COMMIT sine transform st'
stop
endif

call  d_forward_trig_transform (f,desc_h1,ipar,dpar,ir)
if (ir.ne.0)then
  print *, 'Error FORWARD  sine transform  st'
 stop
endif

call free_trig_transform(desc_h1,ipar,ir);
if (ir.ne.0)then
  print *,'Error free fwd sine transform '
  stop
 endif

end subroutine

Thanks!

Jon Clements
  • 138,671
  • 33
  • 247
  • 280
  • What about `ik`? You did not declare it so I assume it is `private`. But remember that private variables enter the `parallel` region uninitialized hence `ik` may be junk. – PetrH Aug 06 '14 at 09:52
  • Also with proper compiler flags you may get an error message. – PetrH Aug 06 '14 at 09:53
  • You should show where all those variables come from. – Vladimir F Героям слава Aug 06 '14 at 09:54
  • 1
    Are you sure it is allowed to call `dst_fwd` from a parallel region? I would bet it isn't. – Vladimir F Героям слава Aug 06 '14 at 09:55
  • 4
    @PetrH: by default variables are `shared`. Of course, if OP has failed to use `implicit none` then `ik` may not be declared, or initialised, at all and be, as you suggest, junk. If OP has omitted `implicit none` OP deserves all the pain and torment arising. – High Performance Mark Aug 06 '14 at 10:23
  • Thank you for your quick response! I guess it is the dst_fwd subroutine that causes problems? The subroutine is declared in the module fourierin. Also when I just perform manipulation in parallel on the variable b there is no crash.... PS: I read that the trigonometric transforms in MKL are thread safe... –  Aug 06 '14 at 14:48
  • ll,nn,mm are the dimensions of the array b(ll,mm,nn,3) and are declared in the module maxwell and intialized upon reading a configuration file. And of course the code works properly without the openMP instructions... :-) –  Aug 06 '14 at 14:53
  • What about `desc_h1`? It is accessed by all threads simultaneously, isn't it? I can imagine things getting really nasty if two threads are racing on the same pointer. – PetrH Aug 07 '14 at 13:32
  • This was the problem. I had defined it globally before. I moved the declaration into the body of the subroutine. Now it works fine. –  Aug 07 '14 at 13:50

1 Answers1

2

The problem is the desc_h1 pointer in module fourierin. When you bring the module in scope by use fourierin then desc_h1 becomes accessible by all threads causing a race condition. Adding desc_h1 in the list of private variables should fix it.

PetrH
  • 700
  • 1
  • 4
  • 12