优化复数的内存访问
Optimizing memory access for complex numbers
我有一个对复数进行运算的内核，我正在加载这样的值：
thrust::complex<float> x = X[tIdx];
其中 X 在全局内存中。当我使用 nvvp 分析这个内核时，我发现它是内存带宽限制的，分析器建议我改进内存访问模式：
Global Load L2 Transactions/Access=8, Ideal Transactions/Access=4
反汇编确认该行确实被分成两个 32 位负载，产生跨步访问模式：
LDG.E R9, [R16];
LDG.E R11, [R16+0x4];
如何将其编译成单个 64 位加载？
可能的解决方案
我意识到这与 this earlier question 密切相关，但建议的解决方案（更改全局内存布局或使用共享内存）似乎不如 64 位加载理想。
NVidia developer blog 建议 reinterpret_cast 为矢量数据类型，例如 float2，但我不太清楚这如何适应指针别名规则。
我还必须承认，这在某种程度上是一个理论问题。对于这个特定的内核，我受到设备内存带宽的限制，因此将 L2 事务的数量减半不会显着提高整体性能。但我预计将来会处理更复杂的数字，如果有一个简单的解决方案，那么我想现在就开始使用它。
这里的基本问题是编译器在生成向量加载和存储指令之前似乎需要对类型进行显式对齐规范。考虑以下简单示例：
class __align__(8) cplx0
{
public:
__device__ __host__ cplx0(float _re, float _img) : re(_re), img(_img) {};
float re, img;
};
class cplx1
{
public:
__device__ __host__ cplx1(float _re, float _img) : re(_re), img(_img) {};
float re, img;
};
template<typename T>
__global__ void memsetkernel(T* out, const T val, int N)
{
int tid = threadIdx.x + blockIdx.x * blockDim.x;
int stride = blockDim.x * gridDim.x;
#pragma unroll 8
for(; tid < N; tid += stride) out[tid] = val;
}
template<typename T>
__global__ void memcpykernel(const T* __restrict__ in, T* __restrict__ out, int N)
{
int tid = threadIdx.x + blockIdx.x * blockDim.x;
int stride = blockDim.x * gridDim.x;
#pragma unroll 8
for(; tid < N; tid += stride) out[tid] = in[tid];
}
template<typename T>
void memcpy(const T* in, T* out, int Nitems)
{
int nthreads = 1024;
int nblocks = 13 * 2; // GTX 970 with 13 SM
memcpykernel<T><<<nblocks, nthreads>>>(in, out, Nitems);
cudaDeviceSynchronize();
}
template<typename T>
void memset(T* in, const T value, int Nitems)
{
int nthreads = 1024;
int nblocks = 13 * 2; // GTX 970 with 13 SM
memsetkernel<T><<<nblocks, nthreads>>>(in, value, Nitems);
cudaDeviceSynchronize();
}
int main(void)
{
const int Nitems = 1 << 24;
typedef cplx0 fcomplex0;
typedef cplx1 fcomplex1;
{
fcomplex0* in;
fcomplex0* out;
cudaMalloc((void **)&in, Nitems * sizeof(fcomplex0));
cudaMalloc((void **)&out, Nitems * sizeof(fcomplex1));
for(int i=0; i<10; i++) {
memset<fcomplex0>(in, fcomplex0(1.0f,1.0f), Nitems);
memcpy<fcomplex0>(in, out, Nitems);
}
cudaFree(in);
cudaFree(out);
}
{
fcomplex1* in;
fcomplex1* out;
cudaMalloc((void **)&in, Nitems * sizeof(fcomplex1));
cudaMalloc((void **)&out, Nitems * sizeof(fcomplex1));
for(int i=0; i<10; i++) {
memset<fcomplex1>(in, fcomplex1(1.0f,1.0f), Nitems);
memcpy<fcomplex1>(in, out, Nitems);
cudaDeviceSynchronize();
}
cudaFree(in);
cudaFree(out);
}
cudaDeviceReset();
return 0;
}
这里我们有两种自制的复杂类型，一种具有显式对齐规范，一种没有。否则它们是相同的。在此测试工具中让它们通过一个简单的 mempcy 和 memset 内核，使我们能够检查每种类型的工具链的代码生成行为并对性能进行基准测试。
首先是代码。对于具有显式 8 字节对齐的 cplx0 class，编译器在两个内核中发出矢量化加载和存储：
memcpykernel
ld.global.nc.v2.f32 {%f5, %f6}, [%rd17];
st.global.v2.f32 [%rd18], {%f5, %f6};
memsetkernel
st.global.v2.f32 [%rd11], {%f1, %f2};
而对于 cplx1 情况，它不会：
memcpykernel
ld.global.nc.f32 %f1, [%rd16];
ld.global.nc.f32 %f2, [%rd16+4];
st.global.f32 [%rd15+4], %f2;
st.global.f32 [%rd15], %f1;
memsetkernel
st.global.f32 [%rd11+4], %f2;
st.global.f32 [%rd11], %f1;
在性能方面，memset 案例（CUDA 8 发布工具包，带有 Linux 367.48 驱动程序的 GTX 970）在性能上存在重要差异：
$ nvprof ./complex_types
==29074== NVPROF is profiling process 29074, command: ./complex_types
==29074== Profiling application: ./complex_types
==29074== Profiling result:
Time(%) Time Calls Avg Min Max Name
33.04% 19.264ms 10 1.9264ms 1.9238ms 1.9303ms void memcpykernel<cplx1>(cplx1 const *, cplx1*, int)
32.72% 19.080ms 10 1.9080ms 1.9055ms 1.9106ms void memcpykernel<cplx0>(cplx0 const *, cplx0*, int)
19.15% 11.165ms 10 1.1165ms 1.1120ms 1.1217ms void memsetkernel<cplx1>(cplx1*, cplx1, int)
15.09% 8.7985ms 10 879.85us 877.67us 884.13us void memsetkernel<cplx0>(cplx0*, cplx0, int)
Thrust 模板化复杂类型没有明确的对齐定义（尽管它可能可以通过专门化来实现，尽管这会在某种程度上破坏目的）。所以你在这里唯一的选择是要么制作你自己的具有显式对齐的 Thrust 类型版本，要么使用另一个复杂类型（如 CUBLAS 和 CUFFT 使用的 cuComplex 类型）。
我有一个对复数进行运算的内核，我正在加载这样的值：
thrust::complex<float> x = X[tIdx];
其中 X 在全局内存中。当我使用 nvvp 分析这个内核时，我发现它是内存带宽限制的，分析器建议我改进内存访问模式：
Global Load L2 Transactions/Access=8, Ideal Transactions/Access=4
反汇编确认该行确实被分成两个 32 位负载，产生跨步访问模式：
LDG.E R9, [R16];
LDG.E R11, [R16+0x4];
如何将其编译成单个 64 位加载？
可能的解决方案
我意识到这与 this earlier question 密切相关，但建议的解决方案（更改全局内存布局或使用共享内存）似乎不如 64 位加载理想。
NVidia developer blog 建议 reinterpret_cast 为矢量数据类型，例如 float2，但我不太清楚这如何适应指针别名规则。
我还必须承认，这在某种程度上是一个理论问题。对于这个特定的内核，我受到设备内存带宽的限制，因此将 L2 事务的数量减半不会显着提高整体性能。但我预计将来会处理更复杂的数字，如果有一个简单的解决方案，那么我想现在就开始使用它。
这里的基本问题是编译器在生成向量加载和存储指令之前似乎需要对类型进行显式对齐规范。考虑以下简单示例：
class __align__(8) cplx0
{
public:
__device__ __host__ cplx0(float _re, float _img) : re(_re), img(_img) {};
float re, img;
};
class cplx1
{
public:
__device__ __host__ cplx1(float _re, float _img) : re(_re), img(_img) {};
float re, img;
};
template<typename T>
__global__ void memsetkernel(T* out, const T val, int N)
{
int tid = threadIdx.x + blockIdx.x * blockDim.x;
int stride = blockDim.x * gridDim.x;
#pragma unroll 8
for(; tid < N; tid += stride) out[tid] = val;
}
template<typename T>
__global__ void memcpykernel(const T* __restrict__ in, T* __restrict__ out, int N)
{
int tid = threadIdx.x + blockIdx.x * blockDim.x;
int stride = blockDim.x * gridDim.x;
#pragma unroll 8
for(; tid < N; tid += stride) out[tid] = in[tid];
}
template<typename T>
void memcpy(const T* in, T* out, int Nitems)
{
int nthreads = 1024;
int nblocks = 13 * 2; // GTX 970 with 13 SM
memcpykernel<T><<<nblocks, nthreads>>>(in, out, Nitems);
cudaDeviceSynchronize();
}
template<typename T>
void memset(T* in, const T value, int Nitems)
{
int nthreads = 1024;
int nblocks = 13 * 2; // GTX 970 with 13 SM
memsetkernel<T><<<nblocks, nthreads>>>(in, value, Nitems);
cudaDeviceSynchronize();
}
int main(void)
{
const int Nitems = 1 << 24;
typedef cplx0 fcomplex0;
typedef cplx1 fcomplex1;
{
fcomplex0* in;
fcomplex0* out;
cudaMalloc((void **)&in, Nitems * sizeof(fcomplex0));
cudaMalloc((void **)&out, Nitems * sizeof(fcomplex1));
for(int i=0; i<10; i++) {
memset<fcomplex0>(in, fcomplex0(1.0f,1.0f), Nitems);
memcpy<fcomplex0>(in, out, Nitems);
}
cudaFree(in);
cudaFree(out);
}
{
fcomplex1* in;
fcomplex1* out;
cudaMalloc((void **)&in, Nitems * sizeof(fcomplex1));
cudaMalloc((void **)&out, Nitems * sizeof(fcomplex1));
for(int i=0; i<10; i++) {
memset<fcomplex1>(in, fcomplex1(1.0f,1.0f), Nitems);
memcpy<fcomplex1>(in, out, Nitems);
cudaDeviceSynchronize();
}
cudaFree(in);
cudaFree(out);
}
cudaDeviceReset();
return 0;
}
这里我们有两种自制的复杂类型，一种具有显式对齐规范，一种没有。否则它们是相同的。在此测试工具中让它们通过一个简单的 mempcy 和 memset 内核，使我们能够检查每种类型的工具链的代码生成行为并对性能进行基准测试。
首先是代码。对于具有显式 8 字节对齐的 cplx0 class，编译器在两个内核中发出矢量化加载和存储：
memcpykernel
ld.global.nc.v2.f32 {%f5, %f6}, [%rd17];
st.global.v2.f32 [%rd18], {%f5, %f6};
memsetkernel
st.global.v2.f32 [%rd11], {%f1, %f2};
而对于 cplx1 情况，它不会：
memcpykernel
ld.global.nc.f32 %f1, [%rd16];
ld.global.nc.f32 %f2, [%rd16+4];
st.global.f32 [%rd15+4], %f2;
st.global.f32 [%rd15], %f1;
memsetkernel
st.global.f32 [%rd11+4], %f2;
st.global.f32 [%rd11], %f1;
在性能方面，memset 案例（CUDA 8 发布工具包，带有 Linux 367.48 驱动程序的 GTX 970）在性能上存在重要差异：
$ nvprof ./complex_types
==29074== NVPROF is profiling process 29074, command: ./complex_types
==29074== Profiling application: ./complex_types
==29074== Profiling result:
Time(%) Time Calls Avg Min Max Name
33.04% 19.264ms 10 1.9264ms 1.9238ms 1.9303ms void memcpykernel<cplx1>(cplx1 const *, cplx1*, int)
32.72% 19.080ms 10 1.9080ms 1.9055ms 1.9106ms void memcpykernel<cplx0>(cplx0 const *, cplx0*, int)
19.15% 11.165ms 10 1.1165ms 1.1120ms 1.1217ms void memsetkernel<cplx1>(cplx1*, cplx1, int)
15.09% 8.7985ms 10 879.85us 877.67us 884.13us void memsetkernel<cplx0>(cplx0*, cplx0, int)
Thrust 模板化复杂类型没有明确的对齐定义（尽管它可能可以通过专门化来实现，尽管这会在某种程度上破坏目的）。所以你在这里唯一的选择是要么制作你自己的具有显式对齐的 Thrust 类型版本，要么使用另一个复杂类型（如 CUBLAS 和 CUFFT 使用的 cuComplex 类型）。