 program test_barrier
!
!  this program tests the use of BARRIER in an OMP PARALLEL DO region
!
   use omp_lib
   implicit none

   integer :: i, n, id, nstep
   integer*4, parameter :: ni=16, nj=14, nk=27
   real*8 a(ni,nj), b(nj,nk), c(ni,nk)
!
   n  = omp_get_num_threads()
   id = omp_get_thread_num()
   write (*,11) 'outside parallel',n,id
!
   a = 1
   b = 2
!
   call omp_set_num_threads (4)        
!
   do nstep = 8,9
!
     if ( nstep==8) write (*,*) 'using barrier where all threads have 2 passes'
     if ( nstep==9) write (*,*) 'using barrier where only 1 thread has 3 passes'
!
  !$omp parallel default(none) private(i,n,id,c) shared (nstep,a,b)
     n  = omp_get_num_threads()
     id = omp_get_thread_num()
     write (*,11) 'in parallel',n,id
   !$omp do
     do i = 1, nstep
       n  = omp_get_num_threads()
       id = omp_get_thread_num()
       write (*,12) 'in loop',n,id,i
       call matmul_omp ( c,a,b, ni,nj,nk)
       write (*,12) 'work done',n,id,i
     end do
  !$omp end do
 !$omp end parallel

     write (*,*) 'nstep =',nstep,' now complete'
   end do

 11 format (a,' : threads = ',i0,' : thread = ',i0)
 12 format (a,' : threads = ',i0,' : thread = ',i0,' : loop = ',i0)
 end program test_barrier

 subroutine matmul_omp (c,a,b, ni,nj,nk)
   integer*4 ni,nj,nk
   real*8 :: c(ni,nk), a(ni,nj), b(nj,nk)
!
   integer*4 ii,jj,kk
   real*8    s

!  start all threads at same time
 !$omp barrier
   do ii = 1,ni
     do kk = 1,nk
       s = 0
       do jj = 1,nj
         s = s + a(ii,jj)*b(jj,kk)
       end do
       c(ii,kk) = s
     end do
   end do

 end subroutine matmul_omp   
 