使用可变范围时,循环未矢量化
Loop is not vectorized when variable extent is used
版本 A 代码未矢量化,而版本 B 代码已矢量化。
如何使版本 A 向量化并保留变量范围(不使用文字范围)?
嵌套循环用于广播乘法,如 python 和 matlab 的 numpy 库。 numpy库中广播的描述是here.
版本 A 代码(无 std::vector。无矢量化。)
这仅在.L169
中使用了imull (%rsi), %edx
,这不是SIMD指令。
#include <iostream>
#include <stdint.h>
typedef int32_t DATA_TYPE;
template <size_t N>
size_t cal_size(size_t (&Ashape)[N]){
size_t size = 1;
for(size_t i = 0; i < N; ++i) size *= Ashape[i];
return size;
}
template <size_t N>
size_t * cal_stride(size_t (&Ashape)[N] ) {
size_t size = cal_size(Ashape);
size_t * Astride = new size_t[N];
Astride[0] = size/Ashape[0];
for(size_t i = 1; i < N; ++i){
Astride[i] = Astride[i-1]/Ashape[i];
}
return Astride;
}
template <size_t N>
DATA_TYPE * init_data( size_t (&Ashape)[N] ){
size_t size = cal_size(Ashape);
DATA_TYPE * data = new DATA_TYPE[size];
for(size_t i = 0; i < size; ++i){
data[i] = i + 1;
}
return data;
}
template <size_t N>
void print_data(DATA_TYPE * Adata, size_t (&Ashape)[N] ){
size_t size = cal_size(Ashape);
for(size_t i = 0; i < size; ++i){
std::cout << Adata[i] << ", ";
}
std::cout << std::endl;
}
int main(void){
constexpr size_t nd = 3;
size_t Ashape[] = {20,3,4};
size_t Bshape[] = {1,3,1};
auto Astride = cal_stride(Ashape);
auto Bstride = cal_stride(Bshape);
auto Adata = init_data(Ashape);
auto Bdata = init_data(Bshape);
size_t c[nd] = {0,0,0};
///counter
size_t hint[nd] = {0,2,0};
//hint tells which are the broadcasting axes.
size_t A_i, B_i;
for(c[0] = 0; c[0] < Ashape[0]; ++c[0]){ // Ashape as hint[0] = 0
for(c[1] = 0; c[1] < Bshape[1]; ++c[1]){ //Bshape as hint[1] = 2
for(c[2] = 0; c[2] < Ashape[2];++c[2]){ //Asape as hint[2] = 0
A_i = c[0]*Astride[0] + c[1]*Astride[1] + c[2]*Astride[2];
B_i = c[1]*Bstride[1];
Adata[A_i] *= Bdata[B_i];
}
}
}
print_data(Adata, Ashape);
}
版本 B 代码(无 std::vector。文字范围,并且此向量化)
这在.L20
中使用了pmulld %xmm3, %xmm2
,这是一条SIMD指令。
#include <iostream>
#include <stdint.h>
typedef int32_t DATA_TYPE;
void print_data(DATA_TYPE * Adata, size_t size ){
for(size_t i = 0; i < size; ++i){
std::cout << Adata[i] << ", ";
}
std::cout << std::endl;
}
int main(void){
int32_t Adata[240];
int32_t Bdata[3];
size_t A_i, B_i, i,j,k;
for(i = 0; i < 20; ++i){
for(j = 0; j < 3; ++j){
for(k = 0; k < 4;++k){
A_i = i*12 + j*4 + k*1;
B_i = j*1;
Adata[A_i] *= Bdata[B_i];
}
}
}
print_data(Adata, 240);
}
boost multiarray vectorize 但为什么呢?
不确定它是否使用 simd 内存对齐。
#include "boost/multi_array.hpp"
#include <iostream>
int
main () {
// Create a 3D array that is 3 x 4 x 2
int d1,d2,d3;
typedef boost::multi_array<int, 3> array_type;
typedef array_type::index index;
array_type A(boost::extents[d1][d2][d3]);
array_type B(boost::extents[1][d2][1]);
// Assign values to the elements
for(index i = 0; i != d1; ++i)
for(index j = 0; j != d2; ++j)
for(index k = 0; k != d3; ++k)
A[i][j][k] *= B[0][j][0];
for(index i = 0; i != d1; ++i)
for(index j = 0; j != d2; ++j)
for(index k = 0; k != d3; ++k)
std::cout << A[i][j][k];
return 0;
}
2004 pdf 在 gcc.gnu.org 描述了 gcc 的一些循环优化。我希望 "Symbolic Chrecs"(对应于未分析的变量)意味着 gcc 可以融合具有可变范围的嵌套循环。
最后的办法是用元编程实现循环融合。
为了测试循环迭代内和循环迭代之间的内存访问之间的依赖关系,数组索引(在编译的多面体模型中)应采用 i*c0 + j*c1 + k*c2 + ...
形式,其中 i, j, k ...
是循环计数器和 c0, c1, c2 ...
是整数常量。
在您的示例中,显然问题出在 c[2] * Astride[2]
,因为 Astride[2]
不是常数。
确实,留下表达式
A_i = c[0]*Astride[0] + c[1]*Astride[1] + c[2];
允许对循环进行矢量化。
(编辑:注意 X = Adata + c[0]*Astride[0] + c[1]*Astride[1]
是循环不变量,所以在最内层的循环中我们只有数组 X
和索引 c[0]
)。
现在,问题是如何在 AStride[2]
的地方滴入一个常量。幸运的是,从你在cal_stride
中的计算来看,最后一个总是1。所以我们应该完成了。
版本 A 代码未矢量化,而版本 B 代码已矢量化。
如何使版本 A 向量化并保留变量范围(不使用文字范围)?
嵌套循环用于广播乘法,如 python 和 matlab 的 numpy 库。 numpy库中广播的描述是here.
版本 A 代码(无 std::vector。无矢量化。)
这仅在.L169
中使用了imull (%rsi), %edx
,这不是SIMD指令。
#include <iostream>
#include <stdint.h>
typedef int32_t DATA_TYPE;
template <size_t N>
size_t cal_size(size_t (&Ashape)[N]){
size_t size = 1;
for(size_t i = 0; i < N; ++i) size *= Ashape[i];
return size;
}
template <size_t N>
size_t * cal_stride(size_t (&Ashape)[N] ) {
size_t size = cal_size(Ashape);
size_t * Astride = new size_t[N];
Astride[0] = size/Ashape[0];
for(size_t i = 1; i < N; ++i){
Astride[i] = Astride[i-1]/Ashape[i];
}
return Astride;
}
template <size_t N>
DATA_TYPE * init_data( size_t (&Ashape)[N] ){
size_t size = cal_size(Ashape);
DATA_TYPE * data = new DATA_TYPE[size];
for(size_t i = 0; i < size; ++i){
data[i] = i + 1;
}
return data;
}
template <size_t N>
void print_data(DATA_TYPE * Adata, size_t (&Ashape)[N] ){
size_t size = cal_size(Ashape);
for(size_t i = 0; i < size; ++i){
std::cout << Adata[i] << ", ";
}
std::cout << std::endl;
}
int main(void){
constexpr size_t nd = 3;
size_t Ashape[] = {20,3,4};
size_t Bshape[] = {1,3,1};
auto Astride = cal_stride(Ashape);
auto Bstride = cal_stride(Bshape);
auto Adata = init_data(Ashape);
auto Bdata = init_data(Bshape);
size_t c[nd] = {0,0,0};
///counter
size_t hint[nd] = {0,2,0};
//hint tells which are the broadcasting axes.
size_t A_i, B_i;
for(c[0] = 0; c[0] < Ashape[0]; ++c[0]){ // Ashape as hint[0] = 0
for(c[1] = 0; c[1] < Bshape[1]; ++c[1]){ //Bshape as hint[1] = 2
for(c[2] = 0; c[2] < Ashape[2];++c[2]){ //Asape as hint[2] = 0
A_i = c[0]*Astride[0] + c[1]*Astride[1] + c[2]*Astride[2];
B_i = c[1]*Bstride[1];
Adata[A_i] *= Bdata[B_i];
}
}
}
print_data(Adata, Ashape);
}
版本 B 代码(无 std::vector。文字范围,并且此向量化)
这在.L20
中使用了pmulld %xmm3, %xmm2
,这是一条SIMD指令。
#include <iostream>
#include <stdint.h>
typedef int32_t DATA_TYPE;
void print_data(DATA_TYPE * Adata, size_t size ){
for(size_t i = 0; i < size; ++i){
std::cout << Adata[i] << ", ";
}
std::cout << std::endl;
}
int main(void){
int32_t Adata[240];
int32_t Bdata[3];
size_t A_i, B_i, i,j,k;
for(i = 0; i < 20; ++i){
for(j = 0; j < 3; ++j){
for(k = 0; k < 4;++k){
A_i = i*12 + j*4 + k*1;
B_i = j*1;
Adata[A_i] *= Bdata[B_i];
}
}
}
print_data(Adata, 240);
}
boost multiarray vectorize 但为什么呢? 不确定它是否使用 simd 内存对齐。
#include "boost/multi_array.hpp"
#include <iostream>
int
main () {
// Create a 3D array that is 3 x 4 x 2
int d1,d2,d3;
typedef boost::multi_array<int, 3> array_type;
typedef array_type::index index;
array_type A(boost::extents[d1][d2][d3]);
array_type B(boost::extents[1][d2][1]);
// Assign values to the elements
for(index i = 0; i != d1; ++i)
for(index j = 0; j != d2; ++j)
for(index k = 0; k != d3; ++k)
A[i][j][k] *= B[0][j][0];
for(index i = 0; i != d1; ++i)
for(index j = 0; j != d2; ++j)
for(index k = 0; k != d3; ++k)
std::cout << A[i][j][k];
return 0;
}
2004 pdf 在 gcc.gnu.org 描述了 gcc 的一些循环优化。我希望 "Symbolic Chrecs"(对应于未分析的变量)意味着 gcc 可以融合具有可变范围的嵌套循环。
最后的办法是用元编程实现循环融合。
为了测试循环迭代内和循环迭代之间的内存访问之间的依赖关系,数组索引(在编译的多面体模型中)应采用 i*c0 + j*c1 + k*c2 + ...
形式,其中 i, j, k ...
是循环计数器和 c0, c1, c2 ...
是整数常量。
在您的示例中,显然问题出在 c[2] * Astride[2]
,因为 Astride[2]
不是常数。
确实,留下表达式
A_i = c[0]*Astride[0] + c[1]*Astride[1] + c[2];
允许对循环进行矢量化。
(编辑:注意 X = Adata + c[0]*Astride[0] + c[1]*Astride[1]
是循环不变量,所以在最内层的循环中我们只有数组 X
和索引 c[0]
)。
现在,问题是如何在 AStride[2]
的地方滴入一个常量。幸运的是,从你在cal_stride
中的计算来看,最后一个总是1。所以我们应该完成了。