cuda::std::complex division is slower than expected

@jrhemstad edit: This was originally from the Thrust repo about `thrust::complex`. I pointed @NickKarpowicz at `cuda::std::complex` as it will eventually replace `thrust::complex`. @NickKarpowicz reported that `cuda::std::complex` was even slower than `thrust::complex` :upside_down_face:

Hi, I noticed that the division of a (real) double precision number by a thust::complex isn't as fast as it could be. Maybe this is a case of the compiler not optimizing something it should, but there is an easy workaround. I posted about this on the nVidia forum, and they suggested I create an issue here.

It seems to do the operation "literally": first turn the double into a complex number, then divide complex/complex. This can be done more efficiently, in a way that saves a division. I pasted a simple program below that I made to isolate and test this.

Long story short: If I write a function by hand to do the operation without the additional divide, kernel calls just doing this division ~64 million times average 6.22 ms according to the profiler. Doing it with Thrust’s division operator, they take 10.97 ms on average, on a 2080 Super. On a 3060, the number are similar, 13.23 ms vs. 23.46 ms. This is compiling on Windows in Visual Studio 2022, CUDA 11.7.

So one can simply overload the / operator for a bit of a speedup, as:

__device__ thrust::complex operator/(double a, thrust::complex b) {
double divByDenominator = a / (b.real() * b.real() + b.imag() * b.imag());
return thrust::complex(b.real() * divByDenominator, -b.imag() * divByDenominator);
}

I tried this with floats instead of doubles, and it doesn’t seem to matter there. If it turns out to be the case in general and not just a me thing, maybe it’s worth putting something like that explicitly in the library, or maybe the compiler is just currently missing something it shouldn’t? I'm not sure if this is a compiler issue or thrust issue, but one can most easily work around it when interacting with thrust so I put it here...

The code I used for testing is here:

#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include
#include

#define TESTSIZE 64*1048576
#define THREADS_PER_BLOCK 128
#define NLAUNCHES 5

//divide the arrays using thrust standard operator
__global__ void divideWithThrust(double* x, thrust::complex* y, thrust::complex* z) {
unsigned int i = threadIdx.x + blockIdx.x * blockDim.x;
z[i] = x[i] / y[i];
}

//divide the arrays by hand
__global__ void divideDZ(double* x, thrust::complex* y, thrust::complex* z) {
unsigned int i = threadIdx.x + blockIdx.x * blockDim.x;
double divByDenominator = x[i] / (y[i].real() * y[i].real() + y[i].imag() * y[i].imag());
z[i] = thrust::complex(y[i].real() * divByDenominator, -y[i].imag() * divByDenominator);
}

//divide the arrays by explicitly turning the double into a complex double
__global__ void divideDZupcast(double* x, thrust::complex* y, thrust::complex* z) {
unsigned int i = threadIdx.x + blockIdx.x * blockDim.x;
z[i] = thrust::complex(x[i], 0) / y[i];
}

//float math for comparison
__global__ void divideWithThrustFloat(float* x, thrust::complex* y, thrust::complex* z) {
unsigned int i = threadIdx.x + blockIdx.x * blockDim.x;
z[i] = x[i] / y[i];
}

//float by hand for comparison
__global__ void divideFC(float* x, thrust::complex* y, thrust::complex* z) {
unsigned int i = threadIdx.x + blockIdx.x * blockDim.x;
float divByDenominator = x[i] / (y[i].real() * y[i].real() + y[i].imag() * y[i].imag());
z[i] = thrust::complex(y[i].real() * divByDenominator, -y[i].imag() * divByDenominator);
}

//fill arrays
__global__ void initArrays(double* x, thrust::complex* y) {
unsigned int i = threadIdx.x + blockIdx.x * blockDim.x;
x[i] = sin(0.1 * i);
y[i] = thrust::complex(cos(0.2 * i), sin(0.5 * i));
}
__global__ void initArraysFloat(float* x, thrust::complex* y) {
unsigned int i = threadIdx.x + blockIdx.x * blockDim.x;
x[i] = sin(0.1 * i);
y[i] = thrust::complex(cos(0.2 * i), sin(0.5 * i));
}

int main()
{
//first check with doubles
double *x;
thrust::complex *y, *z;
cudaMalloc(&x, TESTSIZE * sizeof(double));
cudaMalloc(&y, TESTSIZE * sizeof(thrust::complex));
cudaMalloc(&z, TESTSIZE * sizeof(thrust::complex));

//divide by hand
initArrays<<>>(x, y);
for (int i = 0; i < NLAUNCHES; i++) {
divideDZ<<>>(x, y, z);
}

//divide with thrust
initArrays<<>>(x, y);
for (int i = 0; i < NLAUNCHES; i++) {
divideWithThrust<<>>(x, y, z);
}

//divide by turning double into complex explicitly
initArrays << > > (x, y);
for (int i = 0; i < NLAUNCHES; i++) {
divideDZupcast << > > (x, y, z);
}

cudaFree(x);
cudaFree(y);
cudaFree(z);

//compare float division
float *xf;
thrust::complex *yf, * zf;
cudaMalloc(&xf, TESTSIZE * sizeof(float));
cudaMalloc(&yf, TESTSIZE * sizeof(thrust::complex));
cudaMalloc(&zf, TESTSIZE * sizeof(thrust::complex));

initArraysFloat<<>>(xf, yf);
for (int i = 0; i < NLAUNCHES; i++) {
divideFC<<>>(xf, yf, zf);
}

initArraysFloat<<>>(xf, yf);
for (int i = 0; i < NLAUNCHES; i++) {
divideWithThrustFloat<<>>(xf, yf, zf);
}

cudaFree(xf);
cudaFree(yf);
cudaFree(zf);

return 0;
}