OpenMP declare reduction: initializer refers to variable which is not 'omp_priv' nor 'omp_orig'

Leaving out non-essential code:

double* max_of_two(double *r,double *n);

int npoints = 1000000;
double* xy = (double*)malloc(2*npoints*sizeof(double));

double zero_coord[2] = {0.0,0.0};
double max_coord[2] = {0.0,0.0};

#pragma omp declare reduction\
  (takemax:double*:omp_out=max_of_two(omp_out,omp_in)) \
  initializer(omp_priv=zero_coord)

 #pragma omp parallel for reduction(takemax:max_coord)
  for (int i=0; i<npoints; i++) {
    max_coord = max_of_two(max_coord,xy+2*i);
  }

This gives the error:

error: '#pragma omp declare reduction' initializer refers to variable 'zero_coord' which is not 'omp_priv' nor 'omp_orig'

Which is strange because in simpler code I can initialize omp_priv=-1 and such.

There are several issues: OpenMP allows to reduce arrays element-wise. For stack arrays, the compiler wants to use an element-wise reduction operation and looks for a reduction for type double. If you want to define the reduction on a pair of doubles, you should use such type.

For such type you can simply initialize with omp_priv = {0,0} or omp_priv = omp_orig (works for max, but not for the general case).

The main reason for discarding variables different from omp_priv and omp_orig is a question of variable scoping. General variables are not necessarily visible in the scope of the initialization. If you want to initialize from a (global) variable, you would need to wrap the variable in a function and call this function from the initializer clause initializer(init(&omp_priv)) or initializer(omp_priv = init()).

The following compiles and executes for me using gcc and clang:

#include <math.h>
#include <stdio.h>
#include <stdlib.h>

typedef struct {
  double x;
  double y;
} pair;

pair max_of_two(pair r, pair n) {
  r.x = r.x > n.x ? r.x : n.x;
  r.y = r.y > n.y ? r.y : n.y;
  return r;
}

int main() {
  int npoints = 100;
  pair *xy = (pair *)malloc(npoints * sizeof(pair));

  pair max_coord = {0.0, 0.0};

#pragma omp declare reduction(takemax                                          \
:pair : omp_out = max_of_two(omp_out, omp_in)) initializer(omp_priv = {0, 0})

#pragma omp parallel for reduction(takemax : max_coord)
  for (int i = 0; i < npoints; i++) {
    xy[i].x = sin(i + .5);
    xy[i].y = cos(i + .5);
    max_coord = max_of_two(max_coord, xy[i]);
  }

  printf("max_coord = %lf, %lf\n", max_coord.x, max_coord.y);
}