如何编写可以自动矢量化或向量化代码

可以自动矢量化的代码

---

使用gcc version 4.3.4的gcc编译器，编译下面的示例代码：

int main(){

  const unsigned int ArraySize = 10000000;

  float* a = new float[ArraySize];

  float* b = new float[ArraySize];

  float* c = new float[ArraySize];

  for (unsigned int j = 0; j< 200 ; j++) // some repetitions

    for ( unsigned int i = 0; i < ArraySize; ++ i)

        c[i] = a[i] * b[i];

}

编译并运行上面的代码：

g++ vectorisable_example.cpp -o vectorisable_example -O2 time ./vectorisable_example

简单，现在我们看下如何告诉 compiler 自动对代码进行矢量化:

g++ vectorisable_example.cpp -o vectorisable_example_vect -O2 -ftree-vectorize time ./vectorisable_example_vect

第二种编译方式，产生的可执行程序比第一种产生的可执行程序运行的更快？我们深入研究下，看看为什么他就这么快？

简介

---

从上世纪九十年代后期开始，英特尔就已经将单指令多数据（SIMD）机器指令集成到其商品CPU产品线中。这些SIMD指令允许一次将相同的数学运算（乘，加...）应用于2,4或甚至更多的数据值上。这组数据值也被称为“vector”（不要与代数中的vector混淆）。理论上矢量计算的性能增益等于CPU可以容纳的矢量单位的数量。

尽管SIMD指令集已经集成到普通CPU中相当长一段时间了，但是这些扩展功能（SIMD指令集）只能通过在C或者C++代码中通过准汇编语言进行使用。像GNU编译器系列（GCC）这样的现代编译器现在可以将常规的C或C ++源代码转换成向量操作。这个过程被称为自动矢量化（auto-vectorization）。
这篇文章，我们将会介绍必要的软件工具和编程技术，并进一步提供利用GCC自动矢量化功能的示例源代码

硬件

---

要估算自动矢量化代码可以带来多大的性能提升，首先需要了解您正在使用的CPU的能力以及CPU可以处理多少个Vector。因此，运行以下命令：

$ cat /proc/cpuinfo...
flags: fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca cmov pat pse36 clflush dts acpi mmx fxsr sse sse2ss ht tm pbe syscall nx rdtscp lm constant_tsc arch_perfmon pebs bts rep_good nopl xtopology nonstop_tsc aperfmperf pni pclmulqdq dtes64 monitor ds_cpl vmx est tm2ssse3cx16 xtpr pdcmsse4_1 sse4_2x2apic popcnt xsaveavxlahf_lm ida arat epb xsaveopt pln pts dts tpr_shadow vnmi flexpriority ept vpid...

在CPU属性中，您将找到有关SIMD功能的信息（以粗体突出显示）。若信息中包含下表列出的字符串，则您的CPU支持(SIMD)功能：

Name	Register Size	Amount of floats	Amount of doubles
mmx (*)	64 bit	0	0
sse2	128 bit	4	2
ssse3	128 bit	4	2
sse4_1	128 bit	4	2
sse4_2	128 bit	4	2
avx	256 bit	8	4

• 英特尔处理器上的第一个矢量单元只能容纳两个迭代器。

所有Intel和AMD 64位CPU都至少支持sse2指令集。 AVX指令集已于2011年推出，采用Sandy Bridge架构，可在最新一代Mac Book Pro笔记本电脑上使用。
AVX指令集可以同时对四个双精度浮点值应用数学运算。因此，与未使用矢量单位的版本相比，应用程序可以实现的最大性能增益是4倍。
Vector单位大小的持续增大，反应了是计算机微体系结构演化的新趋势。十年前，计算机的时钟频率几乎每年翻一番。计算机的时钟频率已经达到了极限，制造商转而开始提供越来越大的矢量单位，并且在每个芯片中提供越来越多的核心数量。

哪些代码可以被auto-vectorized ?

---

为了将计算分布到CPU的矢量单元，编译器必须对源代码的依赖性和相互影响有一个很好的理解。但最重要的是，编译器必须能够检测出那些代码段可以使用SIMD指令进行优化。最简单的情况是对数组执行计算的循环：

for ( unsigned int i = 0; i < ArraySize; ++ i)
{
      c[i] = a[i] * b[i];
}

使用等于或者高于gcc461(在slc5_amd64_gcc461 scram体系结构的服务器中可用)版本的gcc编译器，编译的时候，一定要带上is-ftree-vectorize。

注意:从自动矢量化中获益最多的循环是包含数学运算的循环。只是迭代对象集合的循环不会获利，在大多数情况下甚至不能自动矢量化。

可以在这里找到可以自动矢量化的循环结构列表：http://gcc.gnu.org/projects/tree-ssa/vectorization.html
一旦你使用选项-ftree-vectorizer-verbose = 7编译你的代码，GCC将会给你一个关于你的程序中所有循环的详细报告，以及它们是否已经被自动矢量化了。以下报告是成功对循环进行向量化的结果：

autovect.cpp:66: note: vect_model_load_cost: aligned.

autovect.cpp:66: note: vect_get_data_access_cost: inside_cost = 1, outside_cost = 0.

autovect.cpp:66: note: vect_model_load_cost: aligned.

autovect.cpp:66: note: vect_get_data_access_cost: inside_cost = 2, outside_cost = 0.

autovect.cpp:66: note: vect_model_store_cost: aligned.

autovect.cpp:66: note: vect_get_data_access_cost: inside_cost = 3, outside_cost = 0.

autovect.cpp:66: note: vect_model_load_cost: aligned.

autovect.cpp:66: note: vect_model_load_cost: inside_cost = 1, outside_cost = 0 .

autovect.cpp:66: note: vect_model_load_cost: aligned.

autovect.cpp:66: note: vect_model_load_cost: inside_cost = 1, outside_cost = 0 .

autovect.cpp:66: note: vect_model_simple_cost: inside_cost = 1, outside_cost = 0 .

autovect.cpp:66: note: vect_model_simple_cost: inside_cost = 1, outside_cost = 1 .

autovect.cpp:66: note: vect_model_store_cost: aligned.

autovect.cpp:66: note: vect_model_store_cost: inside_cost = 1, outside_cost = 0 .

autovect.cpp:66: note: Cost model analysis:

Vector inside of loop cost: 5

Vector outside of loop cost: 11

Scalar iteration cost: 5

Scalar outside cost: 0

prologue iterations: 0

epilogue iterations: 2

Calculated minimum iters for profitability: 3

autovect.cpp:66: note:  Profitability threshold = 3

autovect.cpp:66: note: Profitability threshold is 3 loop iterations.

autovect.cpp:66: note: LOOP VECTORIZED.

若一个循环不能被向量化, GCC将会给出原因:

autovect.cpp:133: note: not vectorized: control flow in loop.

上面的输出的含义是：在循环中调用了tostd :: cout，引入了一个控制流，从而无法对该循环进行向量化。

让GCC自动对你的代码进行向量化

---

想要让GCC对自己的代码进行自动向量化，需要使用最新的编译器，建议使用至少GCC 4.6.1版本的编译器，通常来说，版本越新越好。

编译器标示

---

想要打开auto-vectorization，使用标示：

-ftree-vectorize

如果您使用优化级别-O3or进行编译，则隐式启用此选项。
想要得到哪些loop是已经被auto-vectorized和哪些loops没有被成功矢量化以及原因，可以使用选项：

-ftree-vectorizer-verbose=7

关于如何阅读这个输出，请看下面的Real-World代码示例部分。
某些循环结构（例如浮点数的减法）只能在允许编译器改变数学运算顺序的情况下进行矢量化。要做到这一点，需要使用选项：

-ffast-math

如果您使用优化级别-Ofast，则会隐式启用该选项。请注意，fast-math选项会修正修改浮点运算中的错误操作。详情请看http://gcc.gnu.org/onlinedocs/gcc-4.1.1/gcc/Optimize-Options.html
此外，您可以明确指出编译器在生成二进制文件时应使用哪一个SIMD指令集。在编译x86-64架构时，GCC默认使用SSE2指令集。如果要启用AVX指令集，请使用编译器标志：

-mavx

需要注意的是，要运行二进制文件的机器必须支持AVX指令集。您还可以让通过下面的选项，让编译器自己决定使用机器中的哪个指令集：

-march=native

若你想要使用 C++11 的特性，例如：lambda 表达式,一定要确认开启了新的C++ Standard：

-std=gnu++0x

启用auto-vectorization最佳实践

---

使用数组结构

---

想要启用自动向量化，编译器不仅需要能够分析循环迭代，还需要能够分析访问内存的模式。嵌套数组的复杂类型很难或不可能由编译器自动进行矢量化。建议使用简单的c数组orstd :: arrays（在C ++ 11中引入）来保存数据，并允许编译器以连续的方式访问数据。这也将使您的程序能够利用CPU的各种缓存。

如果你需要可变数组，建议使用std :: vectorand 获取指向第一个元素的指针来访问循环中的元素：

std::vectormyData(23);

double * pMyData = &myData.front();

...

pMyData[i] += 5;

数组结构的完整示例可以在下面的“综合代码示例”一节中找到。

限制分支

---

尽力限制for循环内的分支。 GCC能够将一些if语句翻译成向量化的代码，但很有限。这样做时，将评估所有可能的分支，并丢弃未采用分支的值。

尽量简化代码

---

如果在for循环中有复杂的计算，您希望GCC对其进行矢量化，可以考虑使用C ++ 11中的lambda表达式，将它们分解为更简单的代码。在这里可以找到C ++11新功能的介绍:http://en.wikipedia.org/wiki/Anonymous_function#C.2B.2B

这里是一个lambda函数的例子，用来计算一个double变量的平方。

auto kernel_square =

[] // capture nothing

( double const& a) // take 1 parameter by reference

->double // lambda function returns the square

{

return ( a * a );

};

在循环中调用lambda函数：

for ( unsigned int i = 0; i < size; ++ i)

{

ptr_square[i] = kernel_square( ptr_x[i] ) + kernel_square( ptr_y[i] ) + kernel_square( ptr_z[i] ) ;

}

请注意，此代码将由GCC自动矢量化。对lambda函数的调用不会像常规函数调用那样，造成特别大的开销，因为代码会被完全内联。

另一个更全面的例子是lambda函数中的for循环。这个循环也会被GCC自动矢量化：

// Defining a compute kernel to encapsulate a specific computation

auto kernel_multiply =

[ &cFixedMultiply ] // capture the constant by reference of the scope of the lambda expression

( DataArray const& a, DataArray const& b, DataArray & res ) // take 3 parameters by reference

->void // lambda function returns void

{

// simple loop vectorized

for( unsigned int i = 0; i < a.size(); ++ i)

{

res[i] = a[i] * b[i];

}

};

不要调用外部函数

---

调用外部函数，如exp（），log（）等等，会导致循环不能被向量化。因此，您必须决定循环中的数学运算是否足够，从而将循环分解。这意味着，在第一个循环中，计算调用exp（）时，需要使用的参数，并将此结果存储在临时数组中。GCC会自动矢量化这个循环。第二个循环将简单地执行对exp（）的调用并存储结果。

如果要调用的函数是由您控制的，请尽量尝试将此函数指定为C ++ 11 lambda表达式。可以在这里找到该特性的介绍：http://en.wikipedia.org/wiki/Anonymous_function#C.2B.2B

在循环中使用普通的整数计数器

---

使用普通的整数计数器来构建for循环，而不是std :: iterator，因为这些使得GCC难以分析内存访问和循环的属性，如迭代计数。

for (vector::iterator it = y.begin(); it != y.end(); it++)

{

(*it) += 23;

}

变为：

for (unsigned int i = 0; i < y.size(); ++i)

{

y[i] += 23;

}

综合代码示例

---

下面的代码给出了一些关于如何用GCC编写自动矢量化C ++代码的例子。复制并粘贴源代码并使用下面命令进行编译即可：

g++ -Ofast -ftree-vectorizer-verbose=7 -march=native -std=c++11 -o autovect autovect.cpp

请保证使用的gcc版本至少是：GCC 4.7

/*

   Compile with ( use at least gcc 4.7 ):
   g++ -Ofast -ftree-vectorizer-verbose=7 -march=native -std=c++11 -o autovect autovect.cpp

*/

#include <math.h>

#include <string>
#include <iostream>
#include <array>
#include <vector>

// Sturcture-Of-Array to hold coordinates
struct Vector3
{
   std::vector<double> x;
   std::vector<double> y;
   std::vector<double> z;
   
   // final result of the distance calcualtion
   std::vector<double> distance;

   void add( double _x, double _y, double _z )
   {
      x.push_back( _x );
      y.push_back( _y );
      z.push_back( _z );
      distance.push_back( 0.0f );
   }

};

   
int main()
{

// Fixed Size Arrays

   typedef std::array<double, 10> DataArray;

   DataArray vect_a = { 0,1,2,3,4,5,6,7,8,9 };
   DataArray vect_b = {0.5,1,2,3,4,5,6,7,8,9 };
   DataArray vect_res_plain = { 0,0,0,0,0,0,0,0,0,0};   
   DataArray vect_res_lambda = { 0,0,0,0,0,0,0,0,0,0};

   constexpr double cFixedMultiply = 23.5f;

   // simple loop vectorized
   // -- auto-vectorized --
   for( unsigned int i = 0; i < vect_a.size(); ++ i)
   {
      vect_res_plain[i] = vect_a[i]  + vect_b[i];
   }

   // Defining a compute kernel to encapsulate a specific computation
   auto kernel_multiply = 
      [ &cFixedMultiply ] // capture the constant by reference of the scope of the lambda expression
      ( DataArray const& a, DataArray const& b, DataArray & res ) // take 3 parameters by reference
      ->void // lambda function returns void
   {
      // simple loop vectorized
      // -- auto-vectorized --
      for( unsigned int i = 0; i < a.size(); ++ i)
      {
         res[i] = a[i] * b[i] * cFixedMultiply;
      }
   };

   // call the lambda function
   // this call is autovectorized
   kernel_multiply ( vect_a, vect_b, vect_res_lambda );


   // This kernel will be called multiple times and performs the quadrature
   auto kernel_square = 
      [] // capture nothing
      ( double const& a) // take 1 parameter by reference
      ->double // lambda function returns the square
   {
      return ( a * a );
   };

   // create struct and fill with dummy values
   Vector3 v3;
   for ( unsigned int i = 0; i < 50 ; ++ i)
   {
      v3.add( i * 1.1, i * 2.2,  i* 3.3 );   
   }

   // store the size in a local variable. This is needed to GCG
   // can estimate the loop iterations
   const unsigned int size = v3.x.size();

   // -- auto-vectorized --
   for ( unsigned int i = 0; i < size; ++ i)
   {
      v3.distance[i] = sqrt( kernel_square( v3.x[i] ) + kernel_square( v3.y[i] ) + kernel_square( v3.z[i] )) ;
   }

   // output the result, so GCC won't optimize the calculations away
   std::cout << std::endl << "Computation on std::array" << std::endl;
   for( unsigned int i = 0; i < vect_a.size(); ++ i)
   {
      std::cout << vect_res_plain[i] << std::endl;
      std::cout << vect_res_lambda[i] << std::endl;
   }

   std::cout << std::endl << "Computation on Structure-of-Array with variable sized std::vectors" << std::endl;
   for( unsigned int i = 0; i < v3.x.size(); ++ i)
   {
      std::cout << "sqrt( " << v3.x[i] << "^2 + " << v3.y[i] << "^2 + " << v3.z[i] << "^2 ) = " 
                << v3.distance[i] << std::endl;
   }

   return 0;
}

自动矢量化的实际应用

---

在这里，您可以找到一部分利用GCC编译器的自动矢量化功能的应用列表：
The VDT mathematical library

参考: https://twiki.cern.ch/twiki/bin/view/CMSPublic/WorkBookWritingAutovectorizableCode