The index of inner loops, Private or Shared?

Question

This is my code, I want to make it parallel with OpenMP. I have one main loop to make parallel and some inner loops.

Are the indices of inner loops, like p, i or Li private or shared?
What happend if I do not declare the variables as private or shared?

Do you suggest to use the allocatable variables for this parallel loop?

!$OMP PARALLEL DO
  do l = 1,n_rep      
   do p = 1,n_l - 1
    do q = 1,n_l - 1
     do r = 1,n_l - 1 
      Li = (p - 1)*(n_l - 1)**2 + (q - 1)*(n_l - 1) + r 

      alpha(Li) = pi*rand() 
      gamma(Li) = pi*rand() 
      beta(Li) = pi/2*rand()

      R_x(1,1) = 1.d0
      R_x(1,2) = 0.d0 
      R_x(1,3) = 0.d0

      R_x(2,1) = 0.d0
      R_x(2,2) = cos(alpha(Li)) 
      R_x(2,3) = sin(alpha(Li))

      R_x(3,1) = 0.d0
      R_x(3,2) = -sin(alpha(Li)) 
      R_x(3,3) = cos(alpha(Li))

      R_y(1,1) = cos(beta(Li))
      R_y(1,2) = 0.d0 
      R_y(1,3) = -sin(beta(Li))

      R_y(2,1) = 0.d0
      R_y(2,2) = 1.d0 
      R_y(2,3) = 0.d0

      R_y(3,1) = sin(beta(Li)) 
      R_y(3,2) = 0.d0
      R_y(3,3) = cos(beta(Li))

      R_z(1,1) = cos(gamma(Li))
      R_z(1,2) = sin(gamma(Li))
      R_z(1,3) = 0.d0 

      R_z(2,1) = -sin(gamma(Li))
      R_z(2,2) = cos(gamma(Li))
      R_z(2,3) = 0.d0

      R_z(3,1) = 0.d0
      R_z(3,2) = 0.d0
      R_z(3,3) = 1.d0

      R_xy = matmul(R_x,R_y)
      R_xyz = matmul(R_xy,R_z)


      do i = 1,n_f - 1
       do j = 1,n_f - 1
        do k = 1,n_f - 1
         Li = (i - 1)*(n_f - 1)**2 + (j - 1)*(n_f - 1) + k

         cf_x(i) = x_f(i) + (p - 1)*d_l - x_c(p)
         cf_y(j) = y_f(j) + (q - 1)*d_l - y_c(q)
         cf_z(k) = z_f(k) + (r - 1)*d_l - z_c(r)

         x_rotated = R_xyz(1,1)*cf_x(i) + R_xyz(1,2)*cf_y(j)        &
                   + R_xyz(1,3)*cf_z(k)
         y_rotated = R_xyz(2,1)*cf_x(i) + R_xyz(2,2)*cf_y(j)        &
                   + R_xyz(2,3)*cf_z(k)
         z_rotated = R_xyz(3,1)*cf_x(i) + R_xyz(3,2)*cf_y(j)        &
                   + R_xyz(3,3)*cf_z(k)

        enddo
       enddo
      enddo


     enddo 
    enddo
   enddo
  enddo

!$OMP END PARALLEL DO

OpenMP makes do loop indices private by default. Making them shared would create a race condition with unpredictable results. You have many variables which will need to be made private, possibly including allocatable arrays. — tim18

Holmz Holmz · Accepted Answer · 2017-08-05T11:48:55

Personally I would break this problem up a bit.

Size_of_Array = n_l * n_l * n_l
IF(ALLOCATED(Li))     DEALLOCATE(   Li   )
ALLOCATE(    Li     (Size_of_Array))
IF(ALLOCATED(Alpha))  DEALLOCATE(   Alpha)
ALLOCATE    (Alpha  (Size_of_Array))
IF(ALLOCATED(Beta))   DEALLOCATE(   Beta )
ALLOCATE(    Beta   (Size_of_Array))
IF(ALLOCATED(Gamma))  DEALLOCATE(   Gamma)
ALLOCATE(    Gamma  (Size_of_Array))

indexer = 0
do l = 1,n_rep      
  do p = 1,n_l - 1
    do q = 1,n_l - 1
      do r = 1,n_l - 1 
        indexer = indexer + 1
        Li(Indexer) = (p - 1)*(n_l - 1)**2 + (q - 1)*(n_l - 1) + r 
      ENDDO
    ENDDO
  ENDDO
ENDDO

alpha = pi*rand() 
gamma = pi*rand() 
beta  = pi/2*rand()
!?OMP DO PARALLEL 
DO I= 1, SIZE(Li)
  CALL Make_Array(Alpha(I), Beta(I), Gamma(I), MyArray(:,:,I) )
ENDDO
!etc

Basically moving the array to be inside of either an ELEMENTAL FUNCTION or a PURE SUBROUTINE. Then see what it does for speed with inlining and a single parallel do of some sort (OMP or other).

 PURE SUBROUTINE Make_Array(Alpha, Beta, Gamma, MyArray)
 IMPLICIT NONE
 DOUBLE,                INTENT(IN   ) :: Alpha
 DOUBLE,                INTENT(IN   ) :: Beta
 DOUBLE,                INTENT(IN   ) :: Gamma
 DOUBLE, DIMENSION(3,3) INTENT(INOUT) :: MyArray ! Maybe just intent(OUT)?

  R_x(:,:) = 0.d0 
  R_x(1,1) = 1.d0
  R_x(2,2) = cos(alpha) 
  R_x(2,3) = sin(alpha)

  R_x(3,2) = -sin(alpha) 
  R_x(3,3) = cos(alpha)

  R_y(1,1) = cos(beta)
  R_y(1,3) = -sin(beta)

  R_y(2,1) = 0.d0
  R_y(2,2) = 1.d0 
  R_y(2,3) = 0.d0

  R_y(3,1) = sin(beta(Li)) 
  R_y(3,2) = 0.d0
  R_y(3,3) = cos(beta(Li))

  R_z(1,1) = cos(gamma(Li))
  R_z(1,2) = sin(gamma(Li))
  R_z(1,3) = 0.d0 

  R_z(2,1) = -sin(gamma(Li))
  R_z(2,2) = cos(gamma(Li))

  END SUBROUTINE Make_Array

Etc... For other elemental functions or pure subroutines

  R_xy = matmul(R_x,R_y)
  R_xyz = matmul(R_xy,R_ ...

The index of inner loops, Private or Shared?

1 Answers