OpenMP: gcc causa una suma extraña en caso de -march=native (-march=skylake-avx512) y -O3.

El siguiente código tendrá un comportamiento diferente, dependiendo de la optimización aplicada por gcc y de la arquitectura de destino:

#include <omp.h>
#include <stdlib.h>
#include <stdio.h>
#include <time.h>

static void malloc_testvals(int **testvals, int num_tests, int num_threads) {
  for (int i = 0; i < num_tests; i++) {
    testvals[i] = malloc(num_threads * sizeof(int));
  }
}

static void free_testvals(int **testvals, int num_tests) {
  for (int i = 0; i < num_tests; i++) {
    free(testvals[i]);
  }
}

static void make_random_testvals(int **testvals, int *sums, int num_tests, int num_threads) {
  srand(time(NULL));
  for (int i = 0; i < num_tests; i++) {
    sums[i] = 0;
    for (int j = 0; j < num_threads; j++) {
      testvals[i][j] = -100 + rand() % 201;
      sums[i] += testvals[i][j];
    }
  }
}

typedef struct ThreadCommunicator_s ThreadCommunicator;

typedef struct {
  long num_threads;
  ThreadCommunicator **threads;
} Communicator;

typedef struct ThreadCommunicators {
  Communicator *parent;
  long omp_longval;
} ThreadCommunicator;

static void ThreadCommunicator_init(ThreadCommunicator* self, Communicator* parent) {
  *self = (ThreadCommunicator) {
    .parent = parent,
    .omp_longval = 0
  };
}

static void Communicator_init(Communicator* self) {
  self->num_threads = omp_get_max_threads();
  self->threads = malloc(sizeof(ThreadCommunicator *) * self->num_threads);

  for (int rank = 0; rank < self->num_threads; rank++) {
    self->threads[rank] = malloc(sizeof(ThreadCommunicator));
    ThreadCommunicator_init(self->threads[rank], self);
  }
}

static void Communicator_deinit(Communicator* self) {
  for (int rank = 0; rank < self->num_threads; rank++) {
    free(self->threads[rank]);
  }
  free(self->threads);
}

//Sums over all thread-inherent numbers
static long ThreadCommunicator_allreduce_sum_l(ThreadCommunicator* self, long myval) {
  //share my result with others
  self->omp_longval = myval;
  #pragma omp barrier

  #pragma omp single
  {
    printf("self->parent->num_threads = %ld\n", self->parent->num_threads);
    printf("omp_get_num_threads()     = %d\n", omp_get_num_threads());
  }

  //------------------------------------------------------------------------------------------------------------------
  //Error will be miraculously gone if self->parent->num_threads is replaced by omp_get_num_threads().
  //------------------------------------------------------------------------------------------------------------------
  long sum = 0;
  for (int rank = 0; rank < self->parent->num_threads; rank++) {
    sum += self->parent->threads[rank]->omp_longval;
  }
  #pragma omp barrier

  return sum;
}

#define NUM_TESTS 1

int main() {
  Communicator communicator;
  Communicator_init(&communicator);

  int *testvals[NUM_TESTS];
  //solutions
  int sums[NUM_TESTS];
  malloc_testvals(testvals, NUM_TESTS, communicator.num_threads);

  make_random_testvals(testvals, sums, NUM_TESTS, communicator.num_threads);

  unsigned long error = 0;
  #pragma omp parallel
  {
    if (communicator.num_threads != omp_get_num_threads()) {
      printf("This is not supported in this test.\n");
      //due to make_random_testvals works with communicator.num_threads
      //but ThreadCommunicator_allreduce_sum_l does set only the first omp_get_num_threads() values of it
      abort();
    }

    ThreadCommunicator *thread_comm = communicator.threads[omp_get_thread_num()];

    for (int i = 0; i < NUM_TESTS; i++) {
      long thread_sum = ThreadCommunicator_allreduce_sum_l(thread_comm, testvals[i][omp_get_thread_num()]);
      #pragma omp atomic
        error += (unsigned long) labs(thread_sum - sums[i]);
    }
  }

  if (error != 0) {
    printf("Error occurred (error = %lu)!\n", error);
  }

  free_testvals(testvals, NUM_TESTS);

  Communicator_deinit(&communicator);
}

Compilando mediante

gcc -Wall -std=c99 -fopenmp -O3 -march=skylake-avx512

o

gcc -Wall -std=c99 -fopenmp -O3 -march=native

en

Intel(R) Xeon(R) Gold 6230 CPU
y con
gcc (GCC) 8.3.1 20191121 (Red Hat 8.3.1-5)

producirá ejemplarmente esta salida:

self->parent->num_threads = 16
omp_get_num_threads()     = 16
Error occurred (error = 8070309797393041808)!

Curiosamente, este error desaparece si se aplican uno de los siguientes cambios:

• Reemplazar -O3 por -O2
• Eliminar -march=... de la lista de opciones
• Reemplazar self->parent->num_threads con omp_get_num_threads() como se indica en el código.
• Compilar en una máquina diferente con -march=native (aunque por supuesto no tengo una visión completa de qué sistemas se ven afectados y cuáles no).

Pregunto si esto es un error del compilador o si mi código no cumple con el estándar C o la especificación de OpenMP, por ejemplo, debido a una carrera de datos. ¡Cualquier ayuda es muy apreciada!</time.h></stdio.h></stdlib.h></omp.h>