如何在 Pyopencl 中使用复数函数？答案

【问题标题】：How to Use Functions of Complex Numbers In Pyopencl?如何在 Pyopencl 中使用复数函数？
【发布时间】：2022-09-27 19:55:21
【问题描述】：

我一直在尝试使用 pyopencl 的元素智能内核来实现牛顿分形，据我所知，它基本上编译了一段 C 代码，然后在 GPU 上为每个元素并行运行它。我遇到了让它处理复数的问题。 Pyopencl 文档（https://documen.tician.de/pyopencl/array.html）对我来说似乎不清楚，尽管我尝试将此处提到的定义和包含语句添加到序言中。但我仍然不确定如何做一些简单的事情，比如找到复数的绝对值。这是：https://github.com/inducer/pyopencl/blob/main/pyopencl/cl/pyopencl-complex.h 链接在 pyopencl 文档中。我认为它包含 pyopencl 的复杂内容的 .h 文件，但我不确定，因为我在 Python 之外的经验很少。它与我链接的方式意味着我应该能够通读代码以确定实现了哪些功能以及如何调用其中的内容。但是我无法弄清楚如何使用这些功能。我不知道是不是因为我没有正确导入/包含它们，或者他们的呼叫签名的某些部分我丢失了。

mapclstr=\"\"\"
int C = 0;
cdouble_t fval;
cdouble_t fpval;
fval.real=100;
fpval.real=100;
fval.imag=100;
fpval.imag=100;

while ((fval.real*fval.real+fval.imag*fval.imag)>precision && C<N) 
{
    
  fval=_add(_pow( X[i] , 3 ) + _pow( X[i] , 2 ) + X[i],- 1.0);
  fpval=_add(_mul(3,_pow( X[i] , 2 ))  + _mul(2,X[i]),1);
  
  X[i] =_add(X[i],-_divide(fval,fpval));
  C+=1;


} 
\"\"\"
ctx = cl.create_some_context()
queue = cl.CommandQueue(ctx)

A=np.linspace(3,6,1,dtype=np.complex128)#1val for easy testing
res_g = cl.array.to_device(queue, A)
mapcl = ElementwiseKernel(ctx,\"cdouble_t *X,int N,double precision\",mapclstr,\"mapcl\",preamble=\"#define PYOPENCL_DEFINE_CDOUBLE //#include <pyopencl-complex.h>  \")
mapcl(res_g,np.intc(500),np.float64(0.00001)) 
print(res_g.get())

我得到了所有 _ 函数的这个错误，我试过没有 \"_\" 和其他一些变体，但我得到了同样的错误。

Traceback (most recent call last):
  File \"c:/Users/kieran/Documents/Coding projects/Fractal-Viewer/pyopencltest6.py\", line 37, in <module>
    mapcl(res_g,np.intc(500),np.float64(0.00001))
  File \"C:\\ProgramData\\Anaconda3\\envs\\FractalStuff\\lib\\site-packages\\pyopencl\\elementwise.py\", line 249, in __call__
    kernel, arg_descrs = self.get_kernel(use_range)
  File \"C:\\ProgramData\\Anaconda3\\envs\\FractalStuff\\lib\\site-packages\\pytools\\__init__.py\", line 766, in wrapper
    result = function(obj, *args, **kwargs)
  File \"C:\\ProgramData\\Anaconda3\\envs\\FractalStuff\\lib\\site-packages\\pyopencl\\elementwise.py\", line 222, in get_kernel
    use_range=use_range, **self.kwargs)
  File \"C:\\ProgramData\\Anaconda3\\envs\\FractalStuff\\lib\\site-packages\\pyopencl\\elementwise.py\", line 159, in get_elwise_kernel_and_types
    use_range=use_range, loop_prep=loop_prep, **kwargs)
  File \"C:\\ProgramData\\Anaconda3\\envs\\FractalStuff\\lib\\site-packages\\pyopencl\\elementwise.py\", line 109, in get_elwise_program
    return Program(context, source).build(options)
  File \"C:\\ProgramData\\Anaconda3\\envs\\FractalStuff\\lib\\site-packages\\pyopencl\\__init__.py\", line 539, in build
    options_bytes=options_bytes, source=self._source)
  File \"C:\\ProgramData\\Anaconda3\\envs\\FractalStuff\\lib\\site-packages\\pyopencl\\__init__.py\", line 583, in _build_and_catch_errors
    raise err
pyopencl._cl.RuntimeError: clBuildProgram failed: BUILD_PROGRAM_FAILURE - clBuildProgram failed: BUILD_PROGRAM_FAILURE - clBuildProgram failed: BUILD_PROGRAM_FAILURE

Build on <pyopencl.Device \'gfx1010:xnack-\' on \'AMD Accelerated Parallel Processing\' at 0x271b4693a60>:

C:\\Users\\kieran\\AppData\\Local\\Temp\\comgr-26ac13\\input\\CompileSource:37:8: error: implicit declaration of function \'_add\' is invalid in OpenCL
  fval=_add(_pow( X[i] , 3 ) + _pow( X[i] , 2 ) + X[i],- 1.0);
       ^
C:\\Users\\kieran\\AppData\\Local\\Temp\\comgr-26ac13\\input\\CompileSource:37:8: note: did you mean \'hadd\'?
C:\\constructicon\\builds\\gfx\\two\\21.50\\stream\\opencl\\compiler\\llvm\\build\\wNow64a\\B_rel\\tools\\comgr\\<stdin>:9854:23: note: \'hadd\' declared here
ulong16 __ovld __cnfn hadd(ulong16 x, ulong16 y);
                      ^
C:\\Users\\kieran\\AppData\\Local\\Temp\\comgr-26ac13\\input\\CompileSource:37:13: error: implicit declaration of function \'_pow\' is invalid in OpenCL
  fval=_add(_pow( X[i] , 3 ) + _pow( X[i] , 2 ) + X[i],- 1.0);
            ^
C:\\Users\\kieran\\AppData\\Local\\Temp\\comgr-26ac13\\input\\CompileSource:37:13: note: did you mean \'pow\'?
C:\\constructicon\\builds\\gfx\\two\\21.50\\stream\\opencl\\compiler\\llvm\\build\\wNow64a\\B_rel\\tools\\comgr\\<stdin>:8733:22: note: \'pow\' declared here
half16 __ovld __cnfn pow(half16 x, half16 y);
                     ^
C:\\Users\\kieran\\AppData\\Local\\Temp\\comgr-26ac13\\input\\CompileSource:37:49: error: invalid operands to binary expression (\'int\' and \'__global cdouble_t\')
  fval=_add(_pow( X[i] , 3 ) + _pow( X[i] , 2 ) + X[i],- 1.0);
            ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ^ ~~~~
C:\\Users\\kieran\\AppData\\Local\\Temp\\comgr-26ac13\\input\\CompileSource:38:14: error: implicit declaration of function \'_mul\' is invalid in OpenCL
  fpval=_add(_mul(3,_pow( X[i] , 2 ))  + _mul(2,X[i]),1);
             ^
C:\\Users\\kieran\\AppData\\Local\\Temp\\comgr-26ac13\\input\\CompileSource:38:8: error: assigning to \'__private cdouble_t\' from incompatible type \'int\'
  fpval=_add(_mul(3,_pow( X[i] , 2 ))  + _mul(2,X[i]),1);
       ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C:\\Users\\kieran\\AppData\\Local\\Temp\\comgr-26ac13\\input\\CompileSource:40:20: error: implicit declaration of function \'_divide\' is invalid in OpenCL
  X[i] =_add(X[i],-_divide(fval,fpval));
                   ^
C:\\Users\\kieran\\AppData\\Local\\Temp\\comgr-26ac13\\input\\CompileSource:40:8: error: assigning to \'__global cdouble_t\' from incompatible type \'int\'
  X[i] =_add(X[i],-_divide(fval,fpval));
       ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
7 errors generated.
Error: Failed to compile source (from CL or HIP source to LLVM IR).

我在 PyopenCl 文档中苦苦挣扎，opencl 文档似乎假设我对 C 有很好的了解，而我没有。

我的实际要求是有人告诉我我在 pyopencl 的复杂模块上做错了什么，以及我如何做一些事情，比如找到一个模数，或者将一个复数提高到一个幂。在 opencl 内核中操作复数的另一种方法也将受到赞赏。解释这里发生了什么的加分点：https://github.com/inducer/pyopencl/blob/main/pyopencl/cl/pyopencl-complex.h为什么每个函数都是内联定义的并且有这个TP TPROOT##。

标签： python c opencl complex-numbers pyopencl

【解决方案1】：

我想我解决了 TPROOT## 就像一个通用数据类型，所以如果你想要复数的乘法：cdouble 你使用 cdouble_mul 替换 TPROOT## 为你想要的数据类型。

  { \
    return TPROOT##_new( \
        a.real*b.real - a.imag*b.imag, \
        a.real*b.imag + a.imag*b.real); \
  } \```

I spent ages trying to work this out and there was not help anywhere i could find on what `TPROOT##` meant. Gonna leave this here in hopes it saves someone time.

【讨论】：