multiple openmp reduction on stl::vector gets wrong results

[Sean_D_]

Hi,

In the following test code, I use omp declare reduction on stl::vector.  Using icpc 17.0, omp parallel reduction works correctly on a single stl::vector.  However, when I attempt to reduce to two separate vectors, v1 & v2, it returns v1=sum(v1)+sum(v2) and v2=0.0 (where sum is across threads).  By my reading of the standard, I would expect it to return v1=sum(v1) and v2=sum(v2).  This is what GCC 6.3 does.

Can my omp directives be tweaked to work with icpc 17.0?  Or is there another problem?

Thanks in advance,

Sean

(was advised to resubmit this post to C++ forum - original is here:  https://software.intel.com/en-us/forums/intel-moderncode-for-parallel-architectures/topic/733357)

// following on from http://stackoverflow.com/questions/43168661/openmp-and-reduction-on-stdvector
//
#include <algorithm>
#include <functional>
#include <iostream>
#include <vector>
#include <omp.h>

#pragma omp declare reduction(+ : std::vector<double> : std::transform(omp_in.begin(),omp_in.end(),omp_out.begin(),omp_out.begin(),std::plus<double>())) initializer (omp_priv=omp_orig)


int check1(int size){
  std::vector<double> result(size,0.0);

  int npart=100;
  int n;
  {
    int i;
// Works in GCC 6.3, Intel 17.0:
#pragma omp parallel for private(i,n) reduction(+:result)
    for (n=0; n<npart; n++){
      for (i=0; i<size; i++){
	result  += i;
      }
    }
  }
  int fail=0;
  for (int i=0; i<size; i++) {
    if (result != i*100){
      fail=1;
    } 
    std::cout <<"i="<< i<<" "<<result<<std::endl;
  }
  return fail;
}

int check2(int size){
  std::vector<double> result(size,0.0);
  std::vector<double> result2(size,0.0);

  int npart=100;
  int n;
  {
    int i;
// Works in GCC 6.3, fails in Intel 17.0:
//#pragma omp parallel for private(i,n) reduction(+:result),reduction(+:result2)
// Works in GCC 6.3, fails in Intel 17.0:
#pragma omp parallel for private(i,n) reduction(+:result,result2)
    for (n=0; n<npart; n++){
      for (i=0; i<size; i++){
	result  += i;
	result2 += i;
      }
    }
  }
  int fail=0;
  for (int i=0; i<size; i++) {
    if (result != i*100 || result2 != i*100){
      fail=1;
    }
    std::cout <<"i="<< i<<" "<<result<<" "<<result2<<std::endl;
  }
  return fail;
}



int main(int argc, char *argv[]) {

  int size;

  if (argc < 2)
    size = 10;
  else
    size = atoi(argv[1]);

  int fail1=check1(size);
  int fail2=check2(size);

  return fail1+fail2;
}