OpenCV 中的并行计算刍议

作者：江博文（OpenCV中国团队实习生，澳门大学硕士研究生）

在 CPU 主频遭遇瓶颈的当下，为提高软件性能，减少循环次数，需要对代码进行并行优化。一般而言，并行计算会在以下三个不同维度进行：

矢量并行化：利用 CPU 内的矢量寄存器执行 SIMD 运算，单条指令对矢量中的多个数据进行操作，其目的是提高 CPU 单个核心的运算能力。
线程并行化：将单个进程分散到多条合作线程中，共同处理任务，其目的是提高多核 CPU 的运算能力。
分布式内存并行化：利用 MPI 在多个代码之间进行数据通信，其目的是提高众核 CPU 集群的运算能力。

虽然在开启 O2 以上的优化选项的情况下，现代编译器会尝试分析代码中的循环，并启用自动矢量化与自动并行化在矢量与线程两个维度进行加速，但其加速往往比较保守，在追求极致性能的情况下需要程序开发人员手动显式调用矢量化与并行化指令。为降低软件的迁移成本，OpenCV 提供 OpenCV Universal Intrinsics 与 cv::parallel_for_ 等指令帮助软件开发人员在不同平台下优化软件性能。在此前也有文章介绍这两大指令（使用OpenCV中的universal intrinsics为算法提速；OpenCV中的并行计算parallel_for_（1）），本文将就 OpenCV Universal Intrinsics 发表些许浅见。

OpenCV Universal Intrinsics

SIMD，即单指令流多数据流，是用于划分计算机体系结构的费林分类法中的一种类型，可实现一条指令同时对矢量中的多个数据项执行相同操作，适用于大规模数据的简单循环数据并行处理。矢量的长度由支持 SIMD 指令的硬件决定，从理论上来说，相比于普通的标量计算的代码，借助 512 位位宽的 SIMD 指令实现高效矢量化的代码因其可以一次性操纵16个单精度浮点数，可以将面向单精度浮点数计算的性能提高16倍。

OpenCV Universal Intrinsics 旨在提供一个各类计算架构的统一编程模型，将程序开发与计算机底层的 SIMD 指令集相分离，使得应用开发者可以专注于算法。早在 OpenCV 3.0 时期，OpenCV 针对 128 位的 SSE2、NEON 和 VSX 指令进行封装，提出 Universal Intrinsics，在 OpenCV 4.0 时期，OpenCV 进一步提供 128 位的 RVV 指令、 256 位的 AVX2 与 512 位的 AVX512 指令的封装。

OpenCV Universal Intrinsics 的实现

OpenCV Universal Intrinsics 通过 #ifdef 和 Funtion-like Macros 两大预处理器指令为底层的矢量类型与操作提供封装，并在此基础上额外实现部分常用操作，使得不同指令集的操作接口趋于统一，从而达到在编译期自动选择目标平台最优的 SIMD 指令的目的。

在使用 OpenCV Universal Intrinsics 指令前必须先在编译器选项中开启相关 SIMD 指令集的支持（x86/64 平台 gcc 默认开启对 SSE2 指令集的支持）。以 AVX2 为例，如果在编译器选项中开启 -mavx2 (g++/clang++) 或 /Arch AVX2 (MSVC)，则会定义预处理器符号 __AVX__，OpenCV 在 /modules/core/include/opencv2/core/simd_intrinsics.hpp 对于 __AVX2__ 等预处理器符号进行判别，如果检测到相关的符号已被定义，则定义 CV_AVX2 宏，并包含相应的头文件。

// GCC/Clang: -mavx2
// MSVC: /arch:AVX2
#if defined __AVX2__
#  include <immintrin.h>
#  undef CV_AVX2
#  define CV_AVX2 1
#endif

在 g++/clang++ 中也可开启 -march=<arch> 选项，选择目标机器的指令集架构，开启所有支持的指令集，例如笔者的笔记本电脑 CPU 为 AMD Ryzen 4500U，对应的处理器架构为 Zen 2，则通过 -march=znver2 可一次性开启 SSE 4.2、AVX2、半精度浮点数等的支持。

gcc -c -Q -march=native --help=target 命令可查看编译机器所支持的全部指令集，其中 -march=native 选项将使用 CPUID 指令获取编译机器的 CPU 信息。（注意编译机器所支持的指令集架构与目标机器未必一致，如果目标机器调用本不支持的 SIMD 指令集，很可能在运行时产生 Undefined Behavior。）

关于不同 CPU 架构下编译器指令选项的更多内容，可参见https://gcc.gnu.org/onlinedocs/gcc/Submodel-Options.html) 与 MSVC /arch Options。

OpenCV Universal Intrinsics 的数据类型

OpenCV Universal Intrinsics 当前支持 8/16/32/64 位的有符号/无符号整型与 32/64 位的浮点数，不同数据类型支持的运算各不相同，具体可在https://docs.opencv.org/master/df/d91/groupcorehal__intrin.html 中查验是否支持相关操作。

半精度浮点数由于其存储大小上的优势，近年来多见于深度学习模型压缩等应用中，但现有的 SIMD 指令集中并没有直接使用半精度浮点数的运算指令。在开发此类应用时，需要使用 v_float32 v = vx_load_expand(const float16_t* ptr) 将半精度浮点数转换为单精度浮点数，调用单精度浮点数的相关的运算指令，最后使用 v_pack_store(float16_t* ptr, const v_float32& v) 将半精度浮点数转换为单精度浮点数，同时需要在编译器选项中开启 -f16c (g++/clang++)，相关的实现代码如下所示。

#  if defined __F16C__
#    undef CV_FP16
#    define CV_FP16 1
#  endifinline v_float32x8 v256_load_expand(const float16_t* ptr)
{
#if CV_FP16return v_float32x8(_mm256_cvtph_ps(_mm_loadu_si128((const __m128i*)ptr)));
#elsefloat CV_DECL_ALIGNED(32) buf[8];for (int i = 0; i < 8; i++)buf[i] = (float)ptr[i];return v256_load_aligned(buf);
#endif
}inline void v_pack_store(float16_t* ptr, const v_float32x8& a)
{
#if CV_FP16__m128i ah = _mm256_cvtps_ph(a.val, 0);_mm_storeu_si128((__m128i*)ptr, ah);
#elsefloat CV_DECL_ALIGNED(32) buf[8];v_store_aligned(buf, a);for (int i = 0; i < 8; i++)ptr[i] = float16_t(buf[i]);
#endif
}

OpenCV Universal Intrinsics 的常用指令

OpenCV Universal Intrinsics 中常用的数据读写指令分别为 cv::vx_load() 和 cv::vx_store()。

与 Intel Intrinsics 的函数命名方式不同的是默认的 cv::vx_load() 函数是支持非内存对齐的读取函数，而另外提供 cv::vx_load_aligned() 函数供严格内存对齐的数据读取使用。尽管矢量寄存器其内部的基本数据类型的元素一般是内存对齐的，但是在 SIMD 中的数据对齐还需要保证第一个元素的地址与矢量长度相对齐，例如 SSE 指令下存储器地址为 16 字节对齐，而 AVX2 指令下存储器地址为 32 字节对齐（即存储器首个元素地址最后一位为0，倒数第二位为偶数）。

现代处理器架构中，CPU 对于 Unali gned Data 的存取较为宽容，非内存对齐情况下读写也只有在CPU 单次访存的地址超出 Cache Line，发生 Cache Split 的情况下会导致性能损失，CPU 将需要多次读取缓存或内存。鉴于 SSE、AVX2 等指令集均有提供不要求严格内存对齐的读写指令，其性能开销可以接受，所以通常情况下还是建议优先使用 cv::vx_load()。如果想要在应用中严格实现内存对齐以获取更高性能，可以将 CV_STRONG_ALIGNMENT 置为 1 并使用 aligned 后缀的函数，并在使用 __attribute__ ((aligned())) 调用数据对齐，在此情况下将在运行时断言 CV_Assert(cv::isAlignedsizeof(cv::v_float32)(&buf)) 以检查存储器数据是否与矢量寄存器位长对齐。

除从存储器中读取位长与矢量寄存器位长相同的连续数据外，OpenCV Universal Intrinsics 还支持 cv::v_float64 fp(1.0, 2.0) 与 cv::vx_setall_()、cv::vx_setzero() 等寄存器初始化操作，以及从两个地址分别读取一半位长数据的 cv::vx_load_halves(loptr, hiptr)、仅读取一半位长数据的 vx_load_low()、将整型数据转换为其两倍位长的整型类型然后存入寄存器中的 vx_load_expanded()、将整型数据转换为其四倍位长的整型类型然后存入寄存器中的 vx_load_expanded_q() 等非常规的读取操作。值得一提的是，虽然 SIMD 支持从多个存储器地址加载 Gather 指令与将寄存器数据写入不同地址的 Scatter 指令，但在大多数应用中，想要发挥 SIMD 的最大性能，需要将数据以高速缓存友好型的形式整理，并以固定步长的连续读写指令进行数据的读取和存储，非连续高速缓存的 SIMD 内存操作更为复杂，需要更多的硬件资源，执行速度明显变慢。

以 vx 前缀的读写和初始化指令如 cv::vx_load() 相较于固定位长的指令，可根据当前已启用的 SIMD 指令集，重载为最大位长的相应指令（例如在启用 AVX2 指令集的情况下，cv::vx_load() 将重载为 cv::v256_load()），在充分利用矢量寄存器位宽的同时可确保软件具有良好的可移植性。如果直接使用 256 位或 512 位长的指令，则必须在使用前检验 CV_SIMD256 与 CV_SIMD512 等宏的定义，这将增加代码的冗余并提高软件开发的心智负担。

OpenCV Universal Intrinsics 的运算指令

OpenCV Universal Intrinsics 中的运算指令命名与 Intel Intrinsic 相仿，一般由三部分组成，每部分使用 "_" 分隔，第一部分通常为 v、v128、v256、v512，第二部分为操作函数名称，第三部分为操作的对象名及数据类型。

SIMD 中的算术运算可分为以下两种类型：

竖向运算（Vertical Operations）：由两个矢量生成一个相同位长的矢量的运算，包含常见的加减乘除以及位运算、叉乘运算；
横向运算（Horizontal Operations）：由一个矢量生成一个标量的 Reduce 与 Mask 归并运算。

一般情况下，优先使用竖向运算指令，横向运算由于引入寄存器内部数据的依赖，可能涉及数据重排操作，所以其速度明显慢于竖向运算。以下列将数组中所有元素相加的代码为例，在数组长度为 3200000 的情况下，fast_sum 用时 0.77ms，而 slow_sum 用时 1.22ms。

const std::size_t step = sizeof(cv::v_float32) / sizeof(float);
cv::v_float32 vp;float fast_sum(std::vector<float> r)
{t = (double)cv::getTickCount();cv::v_float32 v_sum = cv::vx_setzero_f32();for(auto i = 0; i < r.size(); i+=step){vp = cv::vx_load(&r[i]);v_sum += vp;}float fast_sum = cv::v_reduce_sum(v_sum);t = (double)cv::getTickCount() - t;std::cout << "fast_sum() Time : " << t * 1000 / cv::getTickFrequency() << " ms" << std::endl;return fast_sum;
}float slow_sum(std::vector<float> r)
{t = (double)cv::getTickCount();float slow_sum = 0;for(auto i = 0; i < r.size(); i+=step){vp = cv::vx_load(&r[i]);slow_sum += cv::v_reduce_sum(vp);}t = (double)cv::getTickCount() - t;std::cout << "slow_sum() Time : " << t * 1000 / cv::getTickFrequency() <<" ms" << std::endl;return slow_sum;
}

OpenCV Universal Intrinsics 实现的矩阵乘

以下代码通过 OpenCV Universal Intrinsics 实现二维矩阵乘，在启用 AVX2 指令、O3 优化选项后，相较于简单的三重标量循环，加速比接近8倍，与预估的理论性能相符。

cv::Mat matrix_multipy_simd(const cv::Mat &src1, const cv::Mat &src2)
{CV_Assert(src1.type() == src2.type() && src1.cols == src2.rows);const cv::Mat src2t = src2.t();const std::size_t step = sizeof(cv::v_float32) / sizeof(float);const std::size_t dst_h = src1.rows;const std::size_t dst_w = src2t.rows;std::size_t n = src1.cols;cv::Mat dst = cv::Mat::zeros(dst_h, dst_w, src1.type());for (std::size_t i = 0; i < dst_h; i++){for (std::size_t j = 0; j < dst_w; j++){cv::v_float32 v_sum = cv::vx_setzero_f32();for (std::size_t k = 0; k < n; k += step){cv::v_float32 vp1 = cv::vx_load(src1.ptr<float>(i, k));cv::v_float32 vp2 = cv::vx_load(src2t.ptr<float>(j, k));v_sum += vp1 * vp2;}dst.at<float>(i, j) = cv::v_reduce_sum(v_sum);}}return dst;
}

当循环运行次数不是矢量长度的整数倍时，轻则运算结果出错，重则导致程序崩溃引发 segmentation fault。在此情况下，我们有两种选择：

数据填充：重新安排数据结构，增加循环迭代次数，使得总循环运行数成为矢量长度的整数倍。
生成两个循环：其一是矢量化循环体，其二是标量余数。

cv::Mat matrix_multipy_padding(const cv::Mat &src1, const cv::Mat &src2)
{CV_Assert(src1.type() == src2.type() && src1.cols == src2.rows);double t = (double)cv::getTickCount();const std::size_t step = sizeof(cv::v_float32) / sizeof(float);cv::Mat dst;if (src1.cols % step == 0)dst = matrix_multipy_simd(src1, src2);else{int padding = step - (src1.cols % step);cv::Mat src1_padded;cv::Mat src2_padded;cv::copyMakeBorder(src1, src1_padded, 0, 0, 0, padding, cv::BORDER_CONSTANT, cv::Scalar::all(0));cv::copyMakeBorder(src2, src2_padded, 0, padding, 0, 0, cv::BORDER_CONSTANT, cv::Scalar::all(0));dst = matrix_multipy_simd(src1_padded, src2_padded);}t = (double)cv::getTickCount() - t;std::cout << "matrix_multipy_padding() Time :" << t * 1000 / cv::getTickFrequency() <<"ms" << std::endl;return dst;
}cv::Mat matrix_multipy_divide(const cv::Mat &src1, const cv::Mat &src2)
{double t = (double)cv::getTickCount();const std::size_t step = sizeof(cv::v_float32) / sizeof(float);const std::size_t dst_h = src1.rows;const std::size_t dst_w = src2.cols;std::size_t n = src1.cols;cv::Mat dst = cv::Mat::zeros(dst_h, dst_w, src1.type());const int r = n % step;if (r == 0){dst = matrix_multipy_simd(src1, src2);}else{cv::Mat src2t = src2.t();for (std::size_t i = 0; i < dst_h; i++){for (std::size_t j = 0; j < dst_w; j++){cv::v_float32 v_sum = cv::vx_setzero_f32();for (std::size_t k = 0; k < (n - r); k += step){（cv::v_float32 vp1 = cv::vx_load(src1.ptr<float>(i, k));cv::v_float32 vp2 = cv::vx_load(src2t.ptr<float>(j, k));v_sum += vp1 * vp2;}dst.at<float>(i, j) = cv::v_reduce_sum(v_sum);for (std::size_t k = n - r; k < n; k++)dst.at<float>(i, j) += src1.at<float>(i, k) * src2t.at<float>(j, k);}}}t = (double)cv::getTickCount() - t;std::cout << "matrix_multipy_divide() Time :" << t * 1000 / cv::getTickFrequency() << "ms" << std::endl;return dst;
}

实验结果表明，在嵌套循环结构中，如果内层的循环迭代数较小，使用被称为数据填充的方法一效率更高。（在 [4000, 35] 与 [35, 3000] 大小的矩阵乘中，采取数据填充方法的 matrix_multipy_divide() 用时仅 79ms，而采取双循环体方法的 matrix_multipy_divide() 用时 122ms。）而在内层循环的迭代数量足够大的情况下，二者的速度则相差无几。

OpenCV Universal Intrinsics 的介绍就此告一段落，由于 SIMD 属数据并行操作且一般不会出现数据依赖，不会引起同步和管理多个线程的开销，在应用过程中的难点主要是如何通过数据转换设计高速缓存友好型的数组结构以充分利用高速缓存的带宽（如上例中转置二维数组）。