OpenMP reduction on operator

Help
1

Somewhere I got the impression that Fortran/OpenMP could do a reduction on operators, just like in C++. The following is from a code that is in my “working examples” folder:

Module Typedef
  Type inttype
     integer :: value = 0
  end type inttype
  Interface operator(+)
     module procedure addints
  end Interface operator(+)
contains
  function addints(i1,i2) result(isum)
    implicit none
    Type(inttype),intent(in) :: i1,i2
    Type(inttype) :: isum
    isum%value = i1%value+i2%value
  end function addints
end Module Typedef

used as

  Type(inttype),dimension(nsize) :: intarray
  Type(inttype) :: intsum = inttype(0)

  !$OMP parallel do reduction(+:intsum)
  do i=1,nsize
     intsum = intsum + intarray(i)
  end do
  !$OMP end parallel do

but right now every compiler gives me some variant of:

   51 |   !$OMP parallel do reduction(+:intsum)
      |                                1
Error: !$OMP DECLARE REDUCTION + not found for type TYPE(inttype) at (1)

Am I mistaken that such a reduction ought to work?

2

Hi Victor @victotronics ,
Welcome to the forum! I believe you also need the declare reduction specifier.

e.g.:

!$omp declare reduction(+:inttype:omp_out=omp_out+omp_in)

See also: omp declare reduction (ufrj.br)

3

Thanks. That worked. But it’s an extra step over the C++ equivalent where declaring the operator is enough.

4

The rule you are referring to for C++ sounds like this one:

If the reduction identifier is an implicitly declared reduction identifier or otherwise not an id-expression then it is implicitly converted to one by prepending the keyword operator (for example, + becomes operator+). This conversion is valid for the +, *, /, && and || operators.

One hiccup in Fortran is there is no guarantee the operator(+) is available to begin with. For example a user could forget to import it:

program foo
  use Typedef, only: inttype
  integer, parameter :: nsize = 50
  type(inttype) :: intsum, intarray(nsize)
  intsum = inttype(0)
  !$omp parallel do reduction(+:intsum)  !! + not available
  do i=1,nsize
     intsum = intsum + intarray(i)       !! + not available
  end do
  !$omp end parallel do
end program

The OpenMP 5.0 standard addresses this as follows:

If the reduction-identifier is the same as the name of a user-defined operator or an extended operator, or the same as a generic name that is one of the allowed intrinsic procedures, and if the operator or procedure name appears in an accessibility statement in the same module, the accessibility of the corresponding declare reduction directive is determined by the accessibility attribute of the statement.

The overhead we are talking about is at most five lines for the set of operators (+, *, /, .and., .or.) where C++ will lookup the corresponding operatorX if it is defined. In Fortran you will need to add a line as @lkedward has shown, preferably close to the interface itself:

  interface operator(+)
    module procedure add
  end interface
  !$omp declare reduction(+:cmplx:omp_out=omp_out+omp_in)

For anyone else curious about this, here is another example:

module cmplx_type
  implicit none
  private
  public :: cmplx, operator(+)
  type :: cmplx
    real :: re = 0, im = 0
  end type
  interface operator(+)
    module procedure add
  end interface
  !$omp declare reduction(+:cmplx:omp_out=omp_out+omp_in)
contains 
  pure function add(a,b) result(c)
    type(cmplx), intent(in) :: a, b
    type(cmplx) :: c
    c = cmplx(a%re + b%re, a%im + b%im)
  end function
end module

program cmplx_demo
  use cmplx_type
  implicit none

  type(cmplx) :: a(5), asum

  integer, parameter :: sp = kind(0.0)
  complex(sp) :: b(5), bsum

  integer :: i

  call random_number(a%re)
  call random_number(a%im)

  asum = cmplx(0.,0.)
  !$omp parallel do reduction(+: asum)
  do i = 1, size(a)
    asum = asum + a(i)
  end do

  b%re = a%re
  b%im = a%im

  bsum = 0
  !$omp parallel do reduction(+: bsum)
  do i = 1, size(b)
    bsum = bsum + b(i)
  end do

  print *, asum
  print *, bsum
end program

I’ve tested it with both gfortran (10 and higer) and ifort (2021.4 and higher). For ifx (version 2022.2.0) it didn’t work and I’ve reported an error.

2 Likes