Amusing bug. This computes the pearson correlation of 2 vectors correctly on a small test case if the omp pragma is present, and incorrectly if it is removed. Do you see why?

cor.cpp

template <typename REAL>
REAL cor(const cpuvec<REAL> &x, const cpuvec<REAL> &y)
{
  const len_t n = x.size();
  if (n != y.size())
    throw std::runtime_error("non-conformal arguments");
  
  const REAL *x_d = x.data_ptr();
  const REAL *y_d = y.data_ptr();
  
  const REAL meanx = x.sum() / n;
  const REAL meany = y.sum() / n;
  REAL cp = 0, normx = 0, normy = 0;
  
  #pragma omp parallel for reduction(+: cp, normx, normy)
  for (len_t i=0; i<n; i++)
  {
    const REAL xi_mm = x_d[i] - meanx;
    const REAL yi_mm = y_d[i] - meany;
    
    cp = xi_mm * yi_mm;
    normx += xi_mm * xi_mm;
    normy += yi_mm * yi_mm;
  }
  
  return cp / sqrt(normx * normy);
}