用OpenMP实现反向查找的并行化

Question

我试图找出与任意色度关系最密切的色温。也就是说，对于以下图中的任何(x，y)点，我想要的是属于Planckian轨迹的最接近的点，从那个点开始，我想要相关的黑体温度：

​

黑体曲线的参数函数是多项式：

#ifdef _OPENMP
#pragma omp declare simd
#endif
static inline void CCT_to_xy_blackbody(const float t, float *x, float *y)
{
  // Take correlated color temperature in K and find the closest blackbody illuminant in 1667 K - 250000 K
  float x_temp = 0.f;
  float y_temp = 0.f;

  if(t >= 1667.f && t <= 4000.f)
    x_temp = -0.2661239e9f / (t * t * t) - 0.2343589e6f / (t * t)  + 0.8776956e3f / t + 0.179910f;
  else if(t > 4000.f && t <= 25000.f)
    x_temp = -3.0258469e9f / (t * t * t) + 2.1070379e6f / (t * t)  + 0.2226347e3f / t + 0.240390f;

  if(t >= 1667.f && t <= 2222.f)
    y_temp = -1.1063814f * x_temp * x_temp * x_temp - 1.34811020f * x_temp * x_temp + 2.18555832f * x_temp - 0.20219683f;
  else if(t > 2222.f && t <= 4000.f)
    y_temp = -0.9549476f * x_temp * x_temp * x_temp - 1.37418593f * x_temp * x_temp + 2.09137015f * x_temp - 0.16748867f;
  else if(t > 4000.f && t <= 25000.f)
    y_temp =  3.0817580f * x_temp * x_temp * x_temp - 5.87338670f * x_temp * x_temp + 3.75112997f * x_temp - 0.37001483f;

  *x = x_temp;
  *y = y_temp;
}

因此，我们的想法是建立一个2D LUT：T -> x_bb, y_bb，测量(x，y)到每一组(x_bb，y_bb)之间的距离，找出这个距离的最小值，并且相应的指数会产生温度。

下面是我同时构建和搜索LUT的函数：

static inline float CCT_reverse_lookup(const float x, const float y)
{
  // Find out the closest correlated color temperature (closest point over the planckian locus)
  // for any arbitrary x, y chromaticity, by brute-force reverse-lookup.
  // Note that the LUT computation could be defered somewhere else, and computed once

  static const float T_min = 1700.f;
  static const float T_max = 25000.f;
  static const float T_range = T_max - T_min;
  static const size_t LUT_samples = 1<<16;

  // just init radius with something big.
  float radius = 2.f;
  float temperature = 0.f;

#ifdef _OPENMP
#pragma omp parallel for simd default(none) \
  firstprivate(x, y) shared(radius, temperature)\
  schedule(simd:static)
#endif
  for(size_t i = 0; i < LUT_samples; i++)
  {
    // we need more values for the low temperatures, so we scale the step with a power
    const float step = powf((float)i / (float)(LUT_samples - 1), 4.0f);

    // Current temperature in the lookup range
    const float T = T_min +  step * T_range;

    // Current x, y chromaticity
    float x_bb, y_bb;
    CCT_to_xy_blackbody(T, &x_bb, &y_bb);

    // Compute distance between current planckian chromaticity and input
    float radius_tmp = hypotf((x_bb - x), (y_bb - y));

    // If we found a smaller radius, save it
    const int match = (radius_tmp < radius);
    radius = (match) ? radius_tmp : radius;
    temperature = (match) ? T : temperature;
  }

  return temperature;
}

因此，在这里，我需要在线程之间共享radius和temperature，它比我想要的要慢。

我知道如果我对最小值感兴趣，我可以使用reduction(min:radius)，所以我想在这里使用类似的简化，以便使radius和temperature在每个线程中都是私有的，然后在最后，返回与所有线程的最小半径相关的温度。

这有可能吗？

Answer 1

当前的代码有一个严重的竞争状况。

// If we found a smaller radius, save it
const int match = (radius_tmp < radius);
radius = (match) ? radius_tmp : radius;
temperature = (match) ? T : temperature;

多个线程可以同时执行这些行，导致radius和temperature的值不同步。它应该是：

#ifdef _OPENMP
#pragma omp critical
#endif
if {radius_tmp < radius) {
  radius = radius_tmp;
  temperature = T;
}

而不是。

无论如何，OpenMP 4.0增加了用户定义的还原操作，所以如果您的编译器至少支持该版本，您可以尝试这样做。下面是一个使用结构包装多个值的示例：

#include <stdio.h>
#include <float.h>

struct pair {
  float radius;
  float temperature;
};

struct pair pair_min(struct pair r, struct pair n) {
  /* r is the current min value, n in the value to compare against it */
  if (n.radius < r.radius) {
    return n;
  } else {
    return r;
  }
}
#ifdef _OPENMP
// Define a new reduction operation
#pragma omp declare reduction \
  (pairmin:struct pair:omp_out=pair_min(omp_out,omp_in))    \
  initializer(omp_priv = { FLT_MAX, 0.0f })
#endif

int main(void) {
  struct pair min_radius = { FLT_MAX, 0.0f };
  struct pair values[4] = {
    {1.0f, 0.1f},
    {2.0f, 0.2f},
    {4.0f, 0.4f},
    {3.0f, 0.3f}
  };

#ifdef _OPENMP
#pragma omp parallel for reduction(pairmin:min_radius)
#endif
  for (int i = 0; i < 4; i++) {
    min_radius = pair_min(min_radius, values[i]);
  }

  printf("{%f, %f}\n", min_radius.radius, min_radius.temperature);
  return 0;
}

有关用户定义的裁减的更多信息，请参见OpenMP 5.0规范的2.19.5.7节(或编译器使用的版本在规范中的等效部分)。