how to do power of complex number in CUBLAS?

Question

I am porting my c++ code to CUDA & CUBLAS. I use stl::complex for complex computation (i.e. pow, log, exp, etc.) but I didn't see the same functions defined in CuComplex library. I don't know how to create those functions but I found some codes online

#include <iostream>
#include <cublas_v2.h>
#include <cuComplex.h>
using namespace std;

typedef cuDoubleComplex Complex;

#define complex(x, y) make_cuDoubleComplex(x, y)

__host__ __device__ double cabs(const Complex& z) {return cuCabs(z);}
__host__ __device__ double carg(const Complex& z) {return atan2(cuCreal(z),   cuCimag(z));}
__host__ __device__ Complex polar(const double &magnitude, const double &angle) {return complex(magnitude*cos(angle), magnitude*sin(angle));}
__host__ __device__ Complex cexp(const Complex& z) {return polar( exp(cuCreal(z)), cuCimag(z));}
__host__ __device__ Complex czlog(const Complex& z) {return complex( ::log(cabs(z)), carg(z) );}
__host__ __device__ Complex cpow(const Complex& z, const int &exponent) {return cexp( cuCmul(czlog(z), complex((double )exponent, 0)) );}

void main(void)
{
  Complex z=complex(0.34, 0.56);
  cout << cuCreal(cpow(z, 2)) << "  " << cuCimag(cpow(z, 2)) << endl;
}

the above results didn't give right answer. Is that anything wrong with the cpow? Is that any better to do the power and other function on complex number?

Answer 1

This is not correct:

__host__ __device__ double carg(const Complex& z) {return atan2(cuCreal(z),   cuCimag(z));}

The polar angle of a complex number is given by the arctangent of the imaginary part divided by the real part of the complex number. This corresponds to the ratio of the first parameter divided by the second parameter of atan2

Therefore you should use:

__host__ __device__ double carg(const Complex& z) {return atan2(cuCimag(z),   cuCreal(z));}

I'm not sure about your power function (cpow) either. Have you tried DeMoivre's theorem? I don't know computationally the best method, but it seems like the first order of business is to get the right answer.

Additional notes:

1. I don't think this question really has anything to do with CUBLAS
2. When posting questions like this, it's helpful if you give the actual results you are observing along with the expected results.

Here's a worked example based on DeMoivre's theorem:

$ cat t233.cu
#include <iostream>
#include <math.h>
#include <cuComplex.h>
#include <complex>

typedef double     rtype;
typedef cuDoubleComplex ctype;
#define rpart(x)   (cuCreal(x))
#define ipart(x)   (cuCimag(x))
#define cmplx(x,y) (make_cuDoubleComplex(x,y))

__host__ __device__ rtype carg(const ctype& z) {return (rtype)atan2(ipart(z), rpart(z));} // polar angle
__host__ __device__ rtype cabs(const ctype& z) {return (rtype)cuCabs(z);}
__host__ __device__ ctype cp2c(const rtype d, const rtype a) {return cmplx(d*cos(a), d*sin(a));}
__host__ __device__ ctype cpow(const ctype& z, const int &n) {return cmplx((pow(cabs(z), n)*cos(n*carg(z))), (pow(cabs(z), n)*sin(n*carg(z))));}

int main(){

  double r = 0.34;
  double i = 0.56;
  int n = 2;

  std::complex<double> stl_num(r,i);
  std::complex<double> cn(n,0);
  ctype cu_num = cmplx(r,i);

  std::complex<double> stl_ans = std::pow(stl_num, cn);
  ctype cu_ans = cpow(cu_num, n);

  std::cout << "STL real: " << std::real(stl_ans) << " STL imag: " << std::imag(stl_ans) << std::endl;
  std::cout << "CU  real: " << rpart(cu_ans) << " CU  imag: " << ipart(cu_ans) << std::endl;
  return 0;
}
$ nvcc -arch=sm_20 -O3 -o t233 t233.cu
$ ./t233
STL real: -0.198 STL imag: 0.3808
CU  real: -0.198 CU  imag: 0.3808
$

I'm not suggesting this is thoroughly tested code, but it seems to be on the right track and gives the right answer for your test case.