0

I'm getting some weird behaviour (most likely caused by me) when running a section of code in Fortran.

This is the section of code:

  SUBROUTINE BMCalc(a_lat_in, unit_vector_in, atoms_in_i, atoms_in_dp, &
    rcut_in, rcut_nl_in, &
    expansion_in, perturb_id, perturb_xyz, &
    result_arr)

    REAL(8), INTENT(IN) :: a_lat_in
    REAL(8), INTENT(IN) :: rcut_in
    REAL(8), INTENT(IN) :: rcut_nl_in
    REAL(8), INTENT(IN), DIMENSION(1:3,1:3) :: unit_vector_in
    INTEGER, INTENT(IN), DIMENSION(:) :: atoms_in_i
    REAL(8), INTENT(IN), DIMENSION(:,:) :: atoms_in_dp
    INTEGER, INTENT(IN), DIMENSION(1:3) :: expansion_in
    INTEGER, INTENT(IN), DIMENSION(:) :: perturb_id
    REAL(8), INTENT(IN), DIMENSION(:,:) :: perturb_xyz
    REAL(8), INTENT(OUT), DIMENSION(:) :: result_arr

    INTEGER :: n
    REAL(8), DIMENSION(1:3,1:3) :: strain
    REAL(8), DIMENSION(1:3,1:3) :: unit_vector
    INTEGER :: threads, thread_id
    REAL(8), DIMENSION(1:size(result_arr,1)) :: t_result_arr

    strain(1,1) = 1.00D0
    strain(1,2) = 0.0D0
    strain(1,3) = 0.0D0
    strain(2,1) = 0.0D0
    strain(2,2) = 1.00D0
    strain(2,3) = 0.0D0
    strain(3,1) = 0.0D0
    strain(3,2) = 0.0D0
    strain(3,3) = 1.00D0


    !$OMP PARALLEL DO PRIVATE(strain, unit_vector, t_result_arr) 
    SHARED(a_lat_in, unit_vector_in)
    Do n=1,17
      strain(1,1) = 1.0D0 + (n-8) * 0.001D0
      strain(2,2) = 1.0D0 + (n-8) * 0.001D0
      strain(3,3) = 1.0D0 + (n-8) * 0.001D0
      unit_vector = matmul(strain,unit_vector_in)
      Call FullCalc(a_lat_in, unit_vector, atoms_in_i, atoms_in_dp, rcut_in, rcut_nl_in, &
      expansion_in, perturb_id, perturb_xyz, 1, t_result_arr)

      result_arr(9000+n) = strain(1,1)
      result_arr(9100+n) = t_result_arr(1)
      print *, t_result_arr(1)
    End Do
    !$OMP END PARALLEL DO

    Do n=1,17
      print *, n, result_arr(9000+n), result_arr(9100+n)
    End Do
  END SUBROUTINE BMCalc

This then prints out:

1  0.99299999999999999        78.651434125940156     
2  0.99399999999999999        78.651434125940156     
"  
16   1.0080000000000000        78.651434125940156     
17   1.0089999999999999        78.651434125940156

If I change the first loop to print out the value, it then prints out the correct values.

Amended code:

  result_arr(9000+n) = strain(1,1)
  result_arr(9100+n) = t_result_arr(1)
  print *, t_result_arr(1)
End Do
!$OMP END PARALLEL DO

Values:

1  0.99299999999999999       -860.32575653894082     
2  0.99399999999999999       -860.44626262594090     
"    
16   1.0080000000000000       -860.08077932360152     
17   1.0089999999999999       -859.91597159744981

I've had this happen in the past and it's most likely something I'm doing wrong. I'm not sure how to describe the problem, but I'm looking for any help or advice.

Vladimir F Героям слава
  • 57,977
  • 4
  • 76
  • 119
benpalmer
  • 21
  • 3
  • 1
    Please post a [mcve] that can actually be compiled and tested. Compile your code with all debugging options like `gfortran -g -fbacktrace -Wall -fcheck=all`. – Vladimir F Героям слава Jun 25 '18 at 11:45
  • 1
    You've made `strain` private. OpenMP doesn't guarantee that private variables start, when the parallel region opens, with the values that the variables of the same name had prior to the opening of the parallel region. As your code is written it seems that the off-diagonal elements of `strain` might have any value at all. I'm not sure that's the cause of the error you report, but it probably doesn't help. – High Performance Mark Jun 25 '18 at 12:34
  • Hi Vladimir and Mark. Thank you for your suggestions. I've set the value of the strain array inside each loop within !$OMP and it now works as expected. – benpalmer Jun 25 '18 at 14:56
  • 1
    `firstprivate` might be part of a better solution, see the first answer to the duplicate question. – High Performance Mark Jun 25 '18 at 15:34

0 Answers0