What is the best way to reduce an array of arrays using OpenMP?

Question

I am using OpenMP with Fortran. I have boiled down my use case to a very simple example. I have an array of objects with a custom derived type, and each object contains an array with a different size. I want to make sure that whatever happens in the loop, I apply a reduction to all the values array components of the vector objects:

program main

  implicit none

  integer :: i

  type vector
     real,allocatable :: values(:)
  end type vector

  type(vector) :: vectors(3)

  allocate(vectors(1)%values(3))
  vectors(1)%values = 0

  allocate(vectors(2)%values(6))
  vectors(2)%values = 0

  allocate(vectors(3)%values(9))
  vectors(3)%values = 0

  !$OMP PARALLEL REDUCTION(+:vectors%values)

  !$OMP DO

  do i=1,1000
     vectors(1)%values = vectors(1)%values + 1
     vectors(2)%values = vectors(2)%values + 2
     vectors(3)%values = vectors(3)%values + 3
  end do

  !$OMP END DO

  !$OMP END PARALLEL

  print*,sum(vectors(1)%values)
  print*,sum(vectors(2)%values)
  print*,sum(vectors(3)%values)

end program main

In this case, REDUCTION(+:vectors%values) doesn't work because I get the following errors:

test2.f90(22): error #6159: A component cannot be an array if the encompassing structure is an array.   [VALUES]
  !$OMP PARALLEL REDUCTION(+:vectors%values)
-------------------------------------^
test2.f90(22): error #7656: Subobjects are not allowed in this OpenMP* clause; a named variable must be specified.   [VECTORS]
  !$OMP PARALLEL REDUCTION(+:vectors%values)
-----------------------------^
compilation aborted for test2.f90 (code 1)

I tried overloading the meaning of + for the vector type and then specifying REDUCTION(+:vectors), but then I still get:

test.f90(43): error #7621: The data type of the variable is not defined for the operator or intrinsic specified on the OpenMP* REDUCTION clause.   [VECTORS]
  !$OMP PARALLEL REDUCTION(+:vectors)
-----------------------------^

What is the recommended way to deal with derives types such as these and getting the reduction to work?

Just for reference, the correct output when compiling without OpenMP is

3000.000    
12000.00    
27000.00

Vladimir F Героям слава · Accepted Answer · 2015-07-20 15:55:54Z

3

This is not just OpenMP problem, you cannot reference vectors%values as a one entity if values is an allocatable array component because rules of Fortran 2003 forbid this. That is because such an array would not have any regular strides in memory, the allocatable components are stored at random adresses.

If the number of elements of the encompassing array is small you can do

  !$OMP PARALLEL REDUCTION(+:vectors(1)%values,vectors(2)%values,vectors(3)%values)

  !$OMP DO

  do i=1,1000
     vectors(1)%values = vectors(1)%values + 1
     vectors(2)%values = vectors(2)%values + 2
     vectors(3)%values = vectors(3)%values + 3
  end do

  !$OMP END DO

  !$OMP END PARALLEL

otherwise you must make another loop, let's say j and make the reduce just vectors(j)%values.

If the compiler does not accept structure components in the reduction clause (have to study the latest standard to see if it hasn't been relaxed), you can make a workaround

  !$OMP PARALLEL
  do j = 1, size(vectors)
    call aux(vectors(j)%values)
  end do
  !$OMP END PARALLEL

contains
  subroutine aux(v)
    real :: v(:)

    !$OMP DO REDUCTION(+:v)
    do i=1,1000
      v = v + j
    end do
    !$OMP END DO
  end subroutine

Associate or pointers would be simpler, but they are not allowed either.

edited Jul 20, 2015 at 15:55

answered Jul 19, 2015 at 16:49

Vladimir F Героям слава

60.7k4 gold badges82 silver badges131 bronze badges

Sign up to request clarification or add additional context in comments.

3 Comments

astrofrog Over a year ago

The vectors(1)%values doesn't seem to be recognized as valid syntax. If I try and compile this with gfortran 4.9.3 I get Syntax error in OpenMP variable list. Does this work for you, and if so, what compiler do you use?

Vladimir F Героям слава Over a year ago

That is because the syntax checker requires a variable name and not expression. You can use associate as a workaround in gfortran-4.8 but later versions don't like it.

Vladimir F Героям слава Over a year ago

@astrofrog see the workaround using an internal procedure.

Community · Accepted Answer · 2017-05-23 12:29:48Z

As an alternative to Vladimir's answer, you can always implement your own reduction using a temporary array and a critical section:

program main

  implicit none

  integer :: i

  type vector
     real,allocatable :: values(:)
  end type vector

  type(vector) :: vectors(3)
  type(vector),allocatable :: tmp(:)

  allocate(vectors(1)%values(3))
  vectors(1)%values = 0

  allocate(vectors(2)%values(6))
  vectors(2)%values = 0

  allocate(vectors(3)%values(9))
  vectors(3)%values = 0

  !$OMP PARALLEL PRIVATE(TMP)

  ! Use a temporary array to hold the local sum
  allocate( tmp(size(vectors)) )
  do i=1,size(tmp)
    allocate( tmp(i)%values( size(vectors(i)%values )) )
    tmp(i)%values = vectors(i)%values
  enddo ! i

  !$OMP DO

  do i=1,1000
     tmp(1)%values = tmp(1)%values + 1
     tmp(2)%values = tmp(2)%values + 2
     tmp(3)%values = tmp(3)%values + 3
  end do

  !$OMP END DO

  ! Get the global sum one thread at a time
  !$OMP CRITICAL
  vectors(1)%values = vectors(1)%values + tmp(1)%values
  vectors(2)%values = vectors(2)%values + tmp(2)%values
  vectors(3)%values = vectors(3)%values + tmp(3)%values
  !$OMP END CRITICAL

  deallocate(tmp)
  !$OMP END PARALLEL

  print*,sum(vectors(1)%values)
  print*,sum(vectors(2)%values)
  print*,sum(vectors(3)%values)

end program main

This snippet could be arranged more efficiently by a loop over all elements of vectors. Then, tmp could be a scalar.

Collectives™ on Stack Overflow

What is the best way to reduce an array of arrays using OpenMP?

2 Answers 2

3 Comments

Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

3 Comments

Comments

Your Answer

Sign up or log in

Post as a guest

Related