为 1 个 GPU 创建 16 位输入到 cufftXtMakePlanMany 和 workSize
creating 16-bit input to cufftXtMakePlanMany and workSize for 1 GPU
我需要对 unsigned int 8 位数据计算 FFT。以前，我使用的是 cufftPlanMany，我的输入是 cufftReal，输出是 cufftComplex，我在 FFT 之前和之后使用转换从无符号 8 位转换为 cufftReal，然后从 cufftComplex 转换为有符号 8 位。
我注意到 cuFFT 有一个很好的选项 运行 FFT 用于半精度数据，我希望它能缩短 运行 宁时间。根据文档，它目前不支持所有的 cudaDataType（如果将来可以的话那就太好了），但至少我可以 运行 它具有 16 位浮点数（半精度）和以下 signature:
cufftResult
cufftXtMakePlanMany(cufftHandle plan, int rank, long long int *n, long long int *inembed,
long long int istride, long long int idist, cudaDataType inputtype,
long long int *onembed, long long int ostride, long long int odist,
cudaDataType outputtype, long long int batch, size_t *workSize,
cudaDataType executiontype);
输入、输出和执行的数据类型分别为： CUDA_R_16F, CUDA_C_16F and CUDA_C_16F。我可以用我的 U8 数据提供这个 cuFFT 是最理想的，有什么办法吗？否则，如果第一次从 U8 转换为 cufftReal 是必要的，我如何将我的数据从 cufftReal 转换为 CUDA_R_16F 然后从 CUDA_C_16F ？ cuda 是否足够聪明，可以将输入从 float 转换为半精度数据，因为 cufftExecR2C ultimaltey 是相同的并且没有其他函数可以调用半精度？
另一个问题是关于 workSize 的，它是为多 GPU 情况设计的。知道如何计算这个尺寸吗？ （我只有 1 个 GPU）。我负责管理那个缓冲区吗？
TL;DR：我可以在这里看到两种可能的方法，一种使用 half-precision 转换，另一种使用 single-precision 转换（可能使用 CUFFT 回调）。选择一个或另一个的原因可能取决于许多因素，例如转换的大小、输入数据范围的控制、您 运行 使用的 GPU，以及其他因素。
我不会尝试解决您在此处指出的输出数据的处理问题：
then from cufftComplex to signed 8bit.
因为在没有更多信息的情况下我不知道该怎么做。但是，每种情况下输入数据的处理都应该说明如何处理输出数据。
- 使用half-precision 转换
这里需要注意的几件事是您 cannot (currently) use callbacks with half-precision transforms, and half-precision transforms can be more sensitive 输入数据特征（例如 DC 偏移、变换大小等）而不是单个或 double-precision 变换。此外，half-precision 大部分转换需要 Pascal 或更新的 GPU（忽略 Jetson 系列）。
因为half-precision转换不支持回调，我们将使用“普通”主机代码来处理输入数据；您也可以在转换之前在设备上执行此处理，提供的代码概述了这两种可能性。我这里的“预处理”主要是为了防止 16 位转换溢出。如果您尝试使用此代码，您将很快看到溢出的样子（输出中的 inf and/or nan）。
$ cat t1961.cu
#include <cufft.h>
#include <stdio.h>
#include <stdlib.h>
#include <cufftXt.h>
#include <cuda_fp16.h>
#include <assert.h>
#include <iostream>
typedef half2 ctype;
typedef half rtype;
typedef unsigned char dtype;
long long sig_size = 1<<18;
const int amplitude = 127;
const float ramplitude = 1/(float)(4*amplitude);
__host__ __device__ half convert(int val){
return __float2half_rn((val - amplitude)*ramplitude);
}
__global__ void dev_convert(rtype *out, dtype *in, int sz){
int idx = threadIdx.x+blockDim.x*blockIdx.x;
if (idx < sz)
out[idx] = convert(in[idx]);
}
int main(){
//put 4x sine waves into a U8 array
dtype *my_data = (dtype *)malloc(sig_size*sizeof(dtype));
for (int i = 0; i < sig_size; i++) my_data[i] = amplitude*(sin((i*8*3.141592654f)/sig_size)+1.0);
rtype *d_idata;
ctype *d_odata;
cudaMalloc(&d_idata, sizeof(rtype)*sig_size);
#ifdef USE_HOST
rtype *h_idata = (rtype *)malloc(sig_size*sizeof(rtype));
// convert to 16 bit float non-offset suitable for cufft
for (int i = 0; i < sig_size; i++) h_idata[i] = convert(my_data[i]);
cudaMemcpy(d_idata, h_idata, sig_size*sizeof(rtype), cudaMemcpyHostToDevice);
#else
const int bs = 256;
dtype *d_mydata;
cudaMalloc(&d_mydata, sig_size*sizeof(dtype));
cudaMemcpy(d_mydata, my_data, sig_size*sizeof(dtype), cudaMemcpyHostToDevice);
dev_convert<<<(sig_size+bs-1)/bs, bs>>>(d_idata, d_mydata, sig_size);
#endif
cudaMalloc(&d_odata, sizeof(ctype)*(sig_size/2+1));
cufftHandle plan;
cufftResult r;
r = cufftCreate(&plan);
assert(r == CUFFT_SUCCESS);
size_t ws = 0;
r = cufftXtMakePlanMany(plan, 1, &sig_size, NULL, 1, 1, CUDA_R_16F, NULL, 1, 1, CUDA_C_16F, 1, &ws, CUDA_C_16F);
assert(r == CUFFT_SUCCESS);
r = cufftXtExec(plan, d_idata, d_odata, CUFFT_FORWARD); // warm-up
assert(r == CUFFT_SUCCESS);
cudaEvent_t start, stop;
cudaEventCreate(&start); cudaEventCreate(&stop);
cudaEventRecord(start);
r = cufftXtExec(plan, d_idata, d_odata, CUFFT_FORWARD);
assert(r == CUFFT_SUCCESS);
cudaEventRecord(stop);
cudaEventSynchronize(stop);
float et;
cudaEventElapsedTime(&et, start, stop);
printf("forward FFT time for %ld samples: %fms\n", sig_size, et);
ctype *h_odata = (ctype *)malloc((sig_size/2+1)*sizeof(ctype));
cudaMemcpy(h_odata, d_odata, (sig_size/2+1)*sizeof(ctype), cudaMemcpyDeviceToHost);
for (int i = 0; i < 8; i++)
std::cout << __half2float(h_odata[i].x) << " + " << __half2float(h_odata[i].y) << "i" << std::endl;
return 0;
}
$ nvcc -o t1961 t1961.cu -lcufft
$ ./t1961
forward FFT time for 262144 samples: 0.027520ms
-258 + 0i
0.00349998 + 0.00127506i
-0.000146866 + -0.000833511i
0.00140095 + -0.00501251i
-1.57031 + -32752i
-0.00198174 + 0.00856018i
0.00474548 + 0.00359917i
-0.00226784 + 0.00987244i
$
- 对负载使用单精度转换 callback
我认为这有一些好处。它不像半精度变换那样容易出现溢出现象，the (load) callback routine allows us to still operate on U8 input data.
$ cat t1962.cu
#include <cufft.h>
#include <stdio.h>
#include <stdlib.h>
#include <cufftXt.h>
#include <cuda_fp16.h>
#include <assert.h>
#include <iostream>
typedef cufftComplex ctype;
typedef cufftReal rtype;
typedef unsigned char dtype;
long long sig_size = 1<<18;
const int amplitude = 127;
const cufftReal ramplitude = 1/(float)(4*amplitude);
__device__ rtype convert(int val){
return (val - amplitude)*ramplitude;
}
__device__ rtype myOwnCallback(void *dataIn,
size_t offset,
void *callerInfo,
void *sharedPtr) {
rtype ret;
ret = convert(((dtype *)dataIn)[offset]);
return ret;
}
__device__ cufftCallbackLoadR myOwnCallbackPtr = myOwnCallback;
int main(){
cufftCallbackLoadR hostCopyOfCallbackPtr;
cudaMemcpyFromSymbol(&hostCopyOfCallbackPtr,
myOwnCallbackPtr,
sizeof(hostCopyOfCallbackPtr));
//put 4x sine waves into a U8 array
dtype *my_data = (dtype *)malloc(sig_size*sizeof(dtype));
for (int i = 0; i < sig_size; i++) my_data[i] = amplitude*(sin((i*8*3.141592654f)/sig_size)+1.0);
ctype *d_odata;
dtype *d_mydata;
cudaMalloc(&d_mydata, sig_size*sizeof(dtype));
cudaMemcpy(d_mydata, my_data, sig_size*sizeof(dtype), cudaMemcpyHostToDevice);
cudaMalloc(&d_odata, sizeof(ctype)*(sig_size/2+1));
cufftHandle plan;
cufftResult r;
r = cufftCreate(&plan);
assert(r == CUFFT_SUCCESS);
size_t ws = 0;
r = cufftXtMakePlanMany(plan, 1, &sig_size, NULL, 1, 1, CUDA_R_32F, NULL, 1, 1, CUDA_C_32F, 1, &ws, CUDA_C_32F);
assert(r == CUFFT_SUCCESS);
void *rps[] = {(void *)hostCopyOfCallbackPtr};
r = cufftXtSetCallback(plan, rps, CUFFT_CB_LD_REAL, NULL);
assert(r == CUFFT_SUCCESS);
r = cufftXtExec(plan, (cufftReal *)d_mydata, d_odata, CUFFT_FORWARD); // warm-up
assert(r == CUFFT_SUCCESS);
cudaEvent_t start, stop;
cudaEventCreate(&start); cudaEventCreate(&stop);
cudaEventRecord(start);
r = cufftXtExec(plan, (cufftReal *)d_mydata, d_odata, CUFFT_FORWARD);
assert(r == CUFFT_SUCCESS);
cudaEventRecord(stop);
cudaEventSynchronize(stop);
float et;
cudaEventElapsedTime(&et, start, stop);
printf("forward FFT time for %ld samples: %fms\n", sig_size, et);
ctype *h_odata = (ctype *)malloc((sig_size/2+1)*sizeof(ctype));
cudaMemcpy(h_odata, d_odata, (sig_size/2+1)*sizeof(ctype), cudaMemcpyDeviceToHost);
for (int i = 0; i < 8; i++)
std::cout << h_odata[i].x << " + " << h_odata[i].y << "i" << std::endl;
return 0;
}
$ nvcc -o t1962 t1962.cu -rdc=true -lcufft_static -lculibos
$ ./t1962
forward FFT time for 262144 samples: 0.031488ms
-257.969 + 0i
0.00344251 + 0.00137726i
-3.96543e-05 + -0.00106905i
0.0013994 + -0.00490045i
0.0331312 + -32759.4i
-0.00190887 + 0.00865401i
0.00454092 + 0.00368094i
-0.00219025 + 0.00983646i
$
是的，两种转换类型的结果在数值上并不相同。期望 16 位浮点计算和 32 位浮点计算相同是不合理的。 32 位计算很可能“更准确”。对于这种正弦波情况，我认为最重要的项是直流项以及基波处的幅度尖峰。这些在数字上彼此接近。其他术语是“在噪音中”。时序结果也不完全可比，因为 16 位计算案例忽略了将数据从 U8 转换为 F16 的内核调用成本。您可以使用探查器或重构代码以获得更具可比性的时间。
workSize对于单GPU情况下使用cufftXtMakePlanMany可以忽略，否则使用provided routines确定workSize。
我需要对 unsigned int 8 位数据计算 FFT。以前，我使用的是 cufftPlanMany，我的输入是 cufftReal，输出是 cufftComplex，我在 FFT 之前和之后使用转换从无符号 8 位转换为 cufftReal，然后从 cufftComplex 转换为有符号 8 位。
我注意到 cuFFT 有一个很好的选项 运行 FFT 用于半精度数据，我希望它能缩短 运行 宁时间。根据文档，它目前不支持所有的 cudaDataType（如果将来可以的话那就太好了），但至少我可以 运行 它具有 16 位浮点数（半精度）和以下 signature:
cufftResult
cufftXtMakePlanMany(cufftHandle plan, int rank, long long int *n, long long int *inembed,
long long int istride, long long int idist, cudaDataType inputtype,
long long int *onembed, long long int ostride, long long int odist,
cudaDataType outputtype, long long int batch, size_t *workSize,
cudaDataType executiontype);
输入、输出和执行的数据类型分别为： CUDA_R_16F, CUDA_C_16F and CUDA_C_16F。我可以用我的 U8 数据提供这个 cuFFT 是最理想的，有什么办法吗？否则，如果第一次从 U8 转换为 cufftReal 是必要的，我如何将我的数据从 cufftReal 转换为 CUDA_R_16F 然后从 CUDA_C_16F ？ cuda 是否足够聪明，可以将输入从 float 转换为半精度数据，因为 cufftExecR2C ultimaltey 是相同的并且没有其他函数可以调用半精度？
另一个问题是关于 workSize 的，它是为多 GPU 情况设计的。知道如何计算这个尺寸吗？ （我只有 1 个 GPU）。我负责管理那个缓冲区吗？
TL;DR：我可以在这里看到两种可能的方法，一种使用 half-precision 转换，另一种使用 single-precision 转换（可能使用 CUFFT 回调）。选择一个或另一个的原因可能取决于许多因素，例如转换的大小、输入数据范围的控制、您 运行 使用的 GPU，以及其他因素。
我不会尝试解决您在此处指出的输出数据的处理问题：
then from cufftComplex to signed 8bit.
因为在没有更多信息的情况下我不知道该怎么做。但是，每种情况下输入数据的处理都应该说明如何处理输出数据。
- 使用half-precision 转换
这里需要注意的几件事是您 cannot (currently) use callbacks with half-precision transforms, and half-precision transforms can be more sensitive 输入数据特征（例如 DC 偏移、变换大小等）而不是单个或 double-precision 变换。此外，half-precision 大部分转换需要 Pascal 或更新的 GPU（忽略 Jetson 系列）。
因为half-precision转换不支持回调，我们将使用“普通”主机代码来处理输入数据；您也可以在转换之前在设备上执行此处理，提供的代码概述了这两种可能性。我这里的“预处理”主要是为了防止 16 位转换溢出。如果您尝试使用此代码，您将很快看到溢出的样子（输出中的 inf and/or nan）。
$ cat t1961.cu
#include <cufft.h>
#include <stdio.h>
#include <stdlib.h>
#include <cufftXt.h>
#include <cuda_fp16.h>
#include <assert.h>
#include <iostream>
typedef half2 ctype;
typedef half rtype;
typedef unsigned char dtype;
long long sig_size = 1<<18;
const int amplitude = 127;
const float ramplitude = 1/(float)(4*amplitude);
__host__ __device__ half convert(int val){
return __float2half_rn((val - amplitude)*ramplitude);
}
__global__ void dev_convert(rtype *out, dtype *in, int sz){
int idx = threadIdx.x+blockDim.x*blockIdx.x;
if (idx < sz)
out[idx] = convert(in[idx]);
}
int main(){
//put 4x sine waves into a U8 array
dtype *my_data = (dtype *)malloc(sig_size*sizeof(dtype));
for (int i = 0; i < sig_size; i++) my_data[i] = amplitude*(sin((i*8*3.141592654f)/sig_size)+1.0);
rtype *d_idata;
ctype *d_odata;
cudaMalloc(&d_idata, sizeof(rtype)*sig_size);
#ifdef USE_HOST
rtype *h_idata = (rtype *)malloc(sig_size*sizeof(rtype));
// convert to 16 bit float non-offset suitable for cufft
for (int i = 0; i < sig_size; i++) h_idata[i] = convert(my_data[i]);
cudaMemcpy(d_idata, h_idata, sig_size*sizeof(rtype), cudaMemcpyHostToDevice);
#else
const int bs = 256;
dtype *d_mydata;
cudaMalloc(&d_mydata, sig_size*sizeof(dtype));
cudaMemcpy(d_mydata, my_data, sig_size*sizeof(dtype), cudaMemcpyHostToDevice);
dev_convert<<<(sig_size+bs-1)/bs, bs>>>(d_idata, d_mydata, sig_size);
#endif
cudaMalloc(&d_odata, sizeof(ctype)*(sig_size/2+1));
cufftHandle plan;
cufftResult r;
r = cufftCreate(&plan);
assert(r == CUFFT_SUCCESS);
size_t ws = 0;
r = cufftXtMakePlanMany(plan, 1, &sig_size, NULL, 1, 1, CUDA_R_16F, NULL, 1, 1, CUDA_C_16F, 1, &ws, CUDA_C_16F);
assert(r == CUFFT_SUCCESS);
r = cufftXtExec(plan, d_idata, d_odata, CUFFT_FORWARD); // warm-up
assert(r == CUFFT_SUCCESS);
cudaEvent_t start, stop;
cudaEventCreate(&start); cudaEventCreate(&stop);
cudaEventRecord(start);
r = cufftXtExec(plan, d_idata, d_odata, CUFFT_FORWARD);
assert(r == CUFFT_SUCCESS);
cudaEventRecord(stop);
cudaEventSynchronize(stop);
float et;
cudaEventElapsedTime(&et, start, stop);
printf("forward FFT time for %ld samples: %fms\n", sig_size, et);
ctype *h_odata = (ctype *)malloc((sig_size/2+1)*sizeof(ctype));
cudaMemcpy(h_odata, d_odata, (sig_size/2+1)*sizeof(ctype), cudaMemcpyDeviceToHost);
for (int i = 0; i < 8; i++)
std::cout << __half2float(h_odata[i].x) << " + " << __half2float(h_odata[i].y) << "i" << std::endl;
return 0;
}
$ nvcc -o t1961 t1961.cu -lcufft
$ ./t1961
forward FFT time for 262144 samples: 0.027520ms
-258 + 0i
0.00349998 + 0.00127506i
-0.000146866 + -0.000833511i
0.00140095 + -0.00501251i
-1.57031 + -32752i
-0.00198174 + 0.00856018i
0.00474548 + 0.00359917i
-0.00226784 + 0.00987244i
$
- 对负载使用单精度转换 callback
我认为这有一些好处。它不像半精度变换那样容易出现溢出现象，the (load) callback routine allows us to still operate on U8 input data.
$ cat t1962.cu
#include <cufft.h>
#include <stdio.h>
#include <stdlib.h>
#include <cufftXt.h>
#include <cuda_fp16.h>
#include <assert.h>
#include <iostream>
typedef cufftComplex ctype;
typedef cufftReal rtype;
typedef unsigned char dtype;
long long sig_size = 1<<18;
const int amplitude = 127;
const cufftReal ramplitude = 1/(float)(4*amplitude);
__device__ rtype convert(int val){
return (val - amplitude)*ramplitude;
}
__device__ rtype myOwnCallback(void *dataIn,
size_t offset,
void *callerInfo,
void *sharedPtr) {
rtype ret;
ret = convert(((dtype *)dataIn)[offset]);
return ret;
}
__device__ cufftCallbackLoadR myOwnCallbackPtr = myOwnCallback;
int main(){
cufftCallbackLoadR hostCopyOfCallbackPtr;
cudaMemcpyFromSymbol(&hostCopyOfCallbackPtr,
myOwnCallbackPtr,
sizeof(hostCopyOfCallbackPtr));
//put 4x sine waves into a U8 array
dtype *my_data = (dtype *)malloc(sig_size*sizeof(dtype));
for (int i = 0; i < sig_size; i++) my_data[i] = amplitude*(sin((i*8*3.141592654f)/sig_size)+1.0);
ctype *d_odata;
dtype *d_mydata;
cudaMalloc(&d_mydata, sig_size*sizeof(dtype));
cudaMemcpy(d_mydata, my_data, sig_size*sizeof(dtype), cudaMemcpyHostToDevice);
cudaMalloc(&d_odata, sizeof(ctype)*(sig_size/2+1));
cufftHandle plan;
cufftResult r;
r = cufftCreate(&plan);
assert(r == CUFFT_SUCCESS);
size_t ws = 0;
r = cufftXtMakePlanMany(plan, 1, &sig_size, NULL, 1, 1, CUDA_R_32F, NULL, 1, 1, CUDA_C_32F, 1, &ws, CUDA_C_32F);
assert(r == CUFFT_SUCCESS);
void *rps[] = {(void *)hostCopyOfCallbackPtr};
r = cufftXtSetCallback(plan, rps, CUFFT_CB_LD_REAL, NULL);
assert(r == CUFFT_SUCCESS);
r = cufftXtExec(plan, (cufftReal *)d_mydata, d_odata, CUFFT_FORWARD); // warm-up
assert(r == CUFFT_SUCCESS);
cudaEvent_t start, stop;
cudaEventCreate(&start); cudaEventCreate(&stop);
cudaEventRecord(start);
r = cufftXtExec(plan, (cufftReal *)d_mydata, d_odata, CUFFT_FORWARD);
assert(r == CUFFT_SUCCESS);
cudaEventRecord(stop);
cudaEventSynchronize(stop);
float et;
cudaEventElapsedTime(&et, start, stop);
printf("forward FFT time for %ld samples: %fms\n", sig_size, et);
ctype *h_odata = (ctype *)malloc((sig_size/2+1)*sizeof(ctype));
cudaMemcpy(h_odata, d_odata, (sig_size/2+1)*sizeof(ctype), cudaMemcpyDeviceToHost);
for (int i = 0; i < 8; i++)
std::cout << h_odata[i].x << " + " << h_odata[i].y << "i" << std::endl;
return 0;
}
$ nvcc -o t1962 t1962.cu -rdc=true -lcufft_static -lculibos
$ ./t1962
forward FFT time for 262144 samples: 0.031488ms
-257.969 + 0i
0.00344251 + 0.00137726i
-3.96543e-05 + -0.00106905i
0.0013994 + -0.00490045i
0.0331312 + -32759.4i
-0.00190887 + 0.00865401i
0.00454092 + 0.00368094i
-0.00219025 + 0.00983646i
$
是的，两种转换类型的结果在数值上并不相同。期望 16 位浮点计算和 32 位浮点计算相同是不合理的。 32 位计算很可能“更准确”。对于这种正弦波情况，我认为最重要的项是直流项以及基波处的幅度尖峰。这些在数字上彼此接近。其他术语是“在噪音中”。时序结果也不完全可比，因为 16 位计算案例忽略了将数据从 U8 转换为 F16 的内核调用成本。您可以使用探查器或重构代码以获得更具可比性的时间。
workSize对于单GPU情况下使用cufftXtMakePlanMany可以忽略，否则使用provided routines确定workSize。