计算科学 - 由直线程序给出的快速评估函数 - 吾爱随笔录

由直线程序给出的快速评估函数

计算科学多项式编译复数部件

2021-12-21 18:55:19

我有一个简单但很长的函数，它接受一个向量 x[10]，并输出一个向量 y[100]。它是一个自动生成的多元多项式评估函数，即只涉及（复杂的）乘法和加法（其中 5000 个），其中有许多变量来缓存可重用的中间结果。结构如下：

static inline void
evaluate(const complex x[10], complex y[100])
{
  const complex &X0 = x[0];   // shortcuts to input
  const complex &X1 = x[1];
  ...
  const complex &X9 = x[9];

  static constexpr complex C0 = 1;   // constants
  static constexpr complex C1 = -1;

  const complex G0 = C1 * X8;        // multiply. G standds for 'gate'
  const complex G1 = C0 + G0;        // add.
  const complex G2 = G1 * X9;
  ...
  const complex G5000 = G3427 + G3430 + G3433 + G3436 + G3438;

  y[0] = G104; // output
  y[1] = G135;
  y[2] = G195;
  ....
}

如何在英特尔处理器上将此功能优化到最大？例如，减少乘法的数量或近似它们会有所帮助，因为复杂的乘法是昂贵的。我试过的：

我-O3 -ffast-math -march=native对 GCC 5-7 使用和其他标志
我对此进行了分析，它主要有 L1 数据写入未命中和一些 L1 数据读取未命中。
想法：门（G）形成一个 DAG，所以如果我将每组自变量递归地打包成向量，那么一个人可能能够利用 SIMD 指令。编译器在 SLP 自动矢量化期间执行此操作，请参阅最近的工作 ( 1 )( 2 )，但我不确定他们处理了多少。

这是一个小例子的 DAG 可视化：

  G2 = G1 * X9;
  G3 = X8 * X6;
  G5 = G1 * X0;
  G6 = X8 * X7;
  G8 = X4 * X4;
  G9 = X5 * X5;
  G11 = X6 * X6;
  G12 = G11 + X6;

这是其中 100 个门的可视化示例：每个等级/线（以及更多）可以一次矢量化。现代系统是否对此进行了优化？

注意：变量const意味着这个直线程序（无循环）采用静态单分配（SSA）形式，这是编译器优化中使用的一个术语。

从 DAG 对节点进行矢量化排序（选择合适的拓扑排序）后，我得到 21 层/秩的可矢量化操作（表中没有的 1 秩对应于种子：常量或向量 X）：该表中的每个条目都可以用 SIMD 向量化，从左到右运行。我将在我的应用程序中执行这个函数 10 亿次。

3个回答

好的，你有一个非常好的问题，我尝试运行一些基准测试。首先，我没有你的参数，所以我用了你的小例子。其次，由于您没有指定语言，因此我使用了 C + GSL（因为我对 C++ 不太熟悉）

#if __STDC_VERSION__ >= 199901L
#define _XOPEN_SOURCE 700
#else
#define _XOPEN_SOURCE 600
#endif /* __STDC_VERSION__ */
#include <time.h>
#include <gsl/gsl_vector.h>
#include <stdio.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_blas.h>

gsl_vector_complex *A;
gsl_vector_complex *B;
gsl_vector_complex *C;

void evaluate(gsl_vector_complex *in, gsl_vector_complex *out)
{
    gsl_complex n_unity;
    gsl_complex unity;
    gsl_complex zero;

    GSL_SET_COMPLEX (&n_unity, -1, 0);
    GSL_SET_COMPLEX (&unity, 1, 0);
    GSL_SET_COMPLEX (&zero, 0, 0);

    const gsl_complex neg_x8 = gsl_complex_mul (n_unity, gsl_vector_complex_get (in, 8));

    gsl_vector_complex_set (A, 0, neg_x8);
    gsl_vector_complex_set (A, 1, neg_x8);
    gsl_vector_complex_set (A, 2, neg_x8);
    gsl_vector_complex_set (A, 3, gsl_vector_complex_get (in, 8));
    gsl_vector_complex_set (A, 4, gsl_vector_complex_get (out, 4));
    gsl_vector_complex_set (A, 5, neg_x8);
    gsl_vector_complex_set (A, 6, gsl_vector_complex_get (in, 8));
    gsl_vector_complex_set (A, 7, gsl_vector_complex_get (out, 7));
    gsl_vector_complex_set (A, 8, gsl_vector_complex_get (in, 4));
    gsl_vector_complex_set (A, 9, gsl_vector_complex_get (in, 5));
    gsl_vector_complex_set (A, 10, gsl_vector_complex_get (out, 10));
    gsl_vector_complex_set (A, 11, gsl_vector_complex_get (in, 6));
    gsl_vector_complex_set (A, 12, gsl_vector_complex_get (in, 6));

    gsl_vector_complex_set (B, 0, unity);
    gsl_vector_complex_set (B, 1, unity);
    gsl_vector_complex_set (B, 2, gsl_vector_complex_get (in, 9));
    gsl_vector_complex_set (B, 3, gsl_vector_complex_get (in, 6));
    gsl_vector_complex_set (B, 4, unity);
    gsl_vector_complex_set (B, 5, gsl_vector_complex_get (in, 0));
    gsl_vector_complex_set (B, 6, gsl_vector_complex_get (in, 7));
    gsl_vector_complex_set (B, 7, unity);
    gsl_vector_complex_set (B, 8, gsl_vector_complex_get (in, 4));
    gsl_vector_complex_set (B, 9, gsl_vector_complex_get (in, 5));
    gsl_vector_complex_set (B, 10, unity);
    gsl_vector_complex_set (B, 11, gsl_vector_complex_get (in, 6));
    gsl_vector_complex_set (B, 12, gsl_vector_complex_get (in, 6));

    gsl_vector_complex_set_basis (C, 1);
    gsl_vector_complex_set (C, 2, gsl_vector_complex_get (in, 9));
    gsl_vector_complex_set (C, 5, gsl_vector_complex_get (in, 0));
    gsl_vector_complex_set (C, 12, gsl_vector_complex_get (in, 6));

    gsl_vector_complex_mul (A, B);
    gsl_blas_zaxpy (unity, C, A);
    gsl_vector_complex_memcpy (out, A);

}

void evaluate_naive(gsl_vector_complex *in, gsl_vector_complex *out)
{
    gsl_complex n_unity;
    gsl_complex unity;
    gsl_complex zero;

    GSL_SET_COMPLEX (&n_unity, -1, 0);
    GSL_SET_COMPLEX (&unity, 1, 0);
    GSL_SET_COMPLEX (&zero, 0, 0);

    const gsl_complex neg_x8 = gsl_complex_mul (n_unity, gsl_vector_complex_get (in, 8));

    gsl_vector_complex_set (out, 0, neg_x8);
    gsl_vector_complex_set (out, 1, gsl_complex_add (unity, neg_x8));
    gsl_vector_complex_set (out, 2, gsl_complex_sub (gsl_vector_complex_get (in, 9), gsl_complex_mul (neg_x8, gsl_vector_complex_get (in, 9))));
    gsl_vector_complex_set (out, 3, gsl_complex_mul (gsl_vector_complex_get (in, 8), gsl_vector_complex_get (in, 6)));
    gsl_vector_complex_set (out, 5, gsl_complex_sub (gsl_vector_complex_get (in, 0), gsl_complex_mul (neg_x8, gsl_vector_complex_get (in, 0))));
    gsl_vector_complex_set (out, 6, gsl_complex_mul (gsl_vector_complex_get (in, 8), gsl_vector_complex_get (in, 7)));
    gsl_vector_complex_set (out, 8, gsl_complex_mul (gsl_vector_complex_get (in, 4), gsl_vector_complex_get (in, 4)));
    gsl_vector_complex_set (out, 9, gsl_complex_mul (gsl_vector_complex_get (in, 5), gsl_vector_complex_get (in, 5)));
    gsl_vector_complex_set (out, 11, gsl_complex_mul (gsl_vector_complex_get (in, 6), gsl_vector_complex_get (in, 6)));
    gsl_vector_complex_set (out, 12, gsl_complex_add (gsl_vector_complex_get (out, 11), gsl_vector_complex_get (in, 6)));
}

int main()
{
    struct timespec start, end;
    size_t time_spent;

    gsl_vector_complex *x;
    gsl_vector_complex *y;

    double x_,y_;
    gsl_complex val;

    const gsl_rng_type * T;
    gsl_rng * r;

    gsl_rng_env_setup();

    T = gsl_rng_default;
    r = gsl_rng_alloc (T);

    x = gsl_vector_complex_calloc (13);
    y = gsl_vector_complex_calloc (13);
    A = gsl_vector_complex_calloc (13);
    B = gsl_vector_complex_calloc (13);
    C = gsl_vector_complex_calloc (13);


    for (size_t iteration = 0; iteration < 100000000; iteration++)
    {
        for (size_t index = 0; index < 13; index++)
        {
            x_ = gsl_ran_gaussian_ziggurat (r, 1.0);
            y_ = gsl_ran_gaussian_ziggurat (r, 1.0);
            GSL_SET_COMPLEX(&val, x_, y_);

            gsl_vector_complex_set (x, index, val);
        }

        for (size_t index = 0; index < 13; index++)
        {
            x_ = gsl_ran_gaussian_ziggurat (r, 1.0);
            y_ = gsl_ran_gaussian_ziggurat (r, 1.0);
            GSL_SET_COMPLEX(&val, x_, y_);

            gsl_vector_complex_set (y, index, val);
        }


        clock_gettime (CLOCK_MONOTONIC, &start);
#ifdef NAIVE
        evaluate_naive (x,y);
#else
        evaluate (x, y);
#endif
        clock_gettime (CLOCK_MONOTONIC, &end);
        time_spent = (end.tv_sec - start.tv_sec) * 1000000000 + (end.tv_nsec - start.tv_nsec);
        printf("%e\n", 1.0 * time_spent);

    }

    gsl_rng_free (r);
    gsl_vector_complex_free(x);
    gsl_vector_complex_free(y);
    gsl_vector_complex_free(A);
    gsl_vector_complex_free(B);
    gsl_vector_complex_free(C);
    return 0;
}

我尝试了矢量化和幼稚的实现，这是我的结果：

gcc, no optimizations, naive version, mean = 270.6; std = 190.6 ns
gcc, no optimizations, vector version, mean = 318.1; std = 381.1 ns
gcc, optimized, naive version, mean = 106.9; std = 160.3 ns
gcc, optimized, vector version, mean = 154.5; std = 109.6 ns
gcc, optimized, machine-specific, naive version, mean = 105.2; std = 96.2 ns
gcc, optimized, machine-specific, vector version, mean = 154.8; std = 144.8 ns

icc, optimized, naive version, mean = 75.6; std = 99.4 ns
icc, optimized, vector version, mean = 164.6; std = 209.9 ns
icc, optimized, machine-specific, naive version, mean = 75.8; std = 189.5 ns
icc, optimized, machine-specific, vector version, mean = 166.9; std = 1850.2 ns // this std is somehow bad, however, I used my laptop to run the tests
icc, no optimizations, naive version, mean = 130.9; std = 153.6 ns
icc, no optimizations, vector version, mean = 266.6; std = 238.8 ns

QED 如果您不确定该怎么做，让编译器 (icc) 发挥它的黑魔法。这种情况可能是特定于平台的，此外，当您增加数组的数量/大小时，时间可能会发生变化。PS 还有一种可能的解决方案可以尝试：Intel 提供了自己的元素向量操作函数，可以改善结果，但是，您必须重写代码。PPS 实际上，我完全忘记了默认情况下 icc 不遵循 IEEE 标准（相当于 -ffast-math），并且可以极快地为您提供极不正确的结果。意识到。

一旦您处于一长串表达式的级别，除了希望编译器在汇编级别找到优化机会之外，您几乎无能为力。这可能会给您带来 10% 或 20% 的收益，但不会产生有意义的影响。

你应该做的是回到你试图实现的原始公式。可能有更有效的方法来看待整个问题，例如通过使用循环、迭代、重新表述数学语句等。这种方法可能能够将你的函数加速十倍——重要的因素。但是在您展示的代码中，这些都无法再访问，并且没有编译器能够撤消为产生您展示的代码而进行的转换。

是否可以并行执行 10 亿次运行？即使没有办法并行化方法体的内部，根据您的问题，您可能会通过并行调用方法来达到某种形式的并发性。如果您不重新使用输出 y[100] 作为函数的输入 x[10]，那么没有人会阻止您在不同的线程/内核/机器等中执行它。

伪代码：

open_one_zillion_threads(); x = get_threadlocal_x(); y = get_threadlocal_y(); evaluate(const x[10], y[100]); store_y_in_proper_place(y[100]); rejoin_threads();

还：

您对输出的精度限制是什么？结果是否需要双精度？你真的需要输出的所有方面吗？在许多问题中，精度和速度之间存在合理的权衡，甚至从双精度切换到单精度可能会提高性能。

另一个想法：

应该可以使用矩阵和向量的形式来表达操作。如果都是简单的操作，您可以将输入数组解释为向量，并对其应用一系列矩阵操作来表示层。这些矩阵将不是方阵：

y [100] = A (B (C (D x [10])))

$y[100] = A( B( C( D ~x[10]) ) )$

这可能是反直觉的，因为您的矩阵中会有大量的零，但它会打开您的问题以使用功能强大且经过良好优化的 LA 包（Lapack / Blas 等），这可能会为您带来开箱即用的并行化在那些矩阵向量产品上，它们可能更擅长使用 SIMD / AVX 等。如果您可以在矩阵向量操作中充分利用 AVX 寄存器，那么其中几个包含零不会对您造成太大伤害。Matrix 方法会将您置于拥有大量在线资料和图书馆的地方，如果您想更改该方法的某些部分，您会更加灵活。

如果您使用的是 intel cpu，那么使用他们自制的 intel 编译器可能会给您另外50%的收益。

其它你可能感兴趣的问题

上一篇用于科学计算的消费类硬件？下一篇高维欠定线性方程的数值解