# CUDA complex powers

**URL:** <https://numba.discourse.group/t/cuda-complex-powers/1158>\
**Category:** Community Support\
**Created:** [January 26, 2022, 11:28pm UTC](https://numba.discourse.group/t/cuda-complex-powers/1158 "2022-01-26T23:28:30Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![joej7](https://yyz2.discourse-cdn.com/free1/user_avatar/numba.discourse.group/joej7/32/514_2.png) [@joej7](https://numba.discourse.group/u/joej7)\
**Post date:** [January 26, 2022, 11:28pm UTC](https://numba.discourse.group/t/cuda-complex-powers/1158/1 "2022-01-26T23:28:30Z")

</div>

Migrating my question from [github](https://github.com/numba/numba/pull/6290)

I’m trying to find the complex powers using numba cuda, and it seems that it has been supported, but for some reason, my code doesn’t run. I’ve implemented my algorithm in CUDA and CPU, Ideally the plots should be the same. I’ve tried this with real numbers and things do work. Open to any suggestions

```python
import numpy as np
import matplotlib.pyplot as plt
import time, math
from matplotlib.colors import hsv_to_rgb
from numba import cuda, jit, prange

@cuda.jit
def get_grid_cuda(Z):
    i, j = cuda.grid(2)
    if i >=0 and i < Z.shape[0] and j >=0 and j < Z.shape[1]:     
        for i in range(200):
            fx=Z[i,j]**3+Z[i,j]**2+Z[i,j]+1
            dfx=3*Z[i,j]**2+2*Z[i,j]+1
            div=fx/dfx
            Z[i,j]=Z[i,j]-div
            if(abs(div)<1e-3):
                break

@jit(nopython=True)
def NewtonR(x):
    for i in range(100):
        fx=x **3+x** 2+x+1
        dfx=3*x**2+2*x+1
        div=fx/dfx
        x=x-div
        if(abs(div)<1e-3):
            break
    return x

@jit(nopython=True, parallel=True)
def get_grid(Z):
    for i in prange(len(Z)):
        for j in prange(len(Z[0])):
            Z[i,j]=NewtonR(Z[i,j])
    return Z

x=np.linspace(-1,1,512,dtype=complex)
xj=np.linspace(-1j,1j,512,dtype=complex)
xx,yy=np.meshgrid(x,xj)
X=xx+yy

Z=np.copy(X)
Z_cpu=np.copy(X)

threadsperblock = (16, 16)
blockspergrid_x = math.ceil(Z.shape[0] / threadsperblock[0])
blockspergrid_y = math.ceil(Z.shape[1] / threadsperblock[1])
blockspergrid = (blockspergrid_x, blockspergrid_y)

get_grid_cuda[blockspergrid, threadsperblock](Z)

Z_cpu=get_grid(Z_cpu)

plt.contourf(X.real,X.imag,Z,levels=30)
plt.show()
plt.close()

plt.contourf(X.real,X.imag,Z_cpu,levels=30)
plt.show()
plt.close()

```

---

<div class="post-metadata">

**Author:** ![stuartarchibald](https://yyz2.discourse-cdn.com/free1/user_avatar/numba.discourse.group/stuartarchibald/32/10_2.png) [@stuartarchibald](https://numba.discourse.group/u/stuartarchibald)\
**Post date:** [January 28, 2022, 12:00pm UTC](https://numba.discourse.group/t/cuda-complex-powers/1158/2 "2022-01-28T12:00:53Z")

</div>

Hi @joej7

I think the issue in the above code is that the induction variable (`i`) in the `for i in range(200)` loop in the CUDA implementation overwrites the CUDA indexing variable `i, j = cuda.grid(2)`, such that data ends up being read from and written to unexpected locations. This seems to work ok:

```auto
import numpy as np
import time, math
from numba import cuda, jit, prange

@cuda.jit
def get_grid_cuda(Z):
    i, j = cuda.grid(2)
    if i >=0 and i < Z.shape[0] and j >=0 and j < Z.shape[1]:
        # NOTE THE USE OF _ AS THE INDUCTION VARIABLE
        for _ in range(200):
            fx=Z[i,j]**3+Z[i,j]**2+Z[i,j]+1
            dfx=3*Z[i,j]**2+2*Z[i,j]+1
            div=fx/dfx
            Z[i,j]=Z[i,j]-div
            if(abs(div)<1e-3):
                break

@jit(nopython=True)
def NewtonR(x):
    for i in range(100):
        fx=x **3+x** 2+x+1
        dfx=3*x**2+2*x+1
        div=fx/dfx
        x=x-div
        if(abs(div)<1e-3):
            break
    return x

@jit(nopython=True, parallel=True)
def get_grid(Z):
    for i in prange(len(Z)):
        for j in prange(len(Z[0])):
            Z[i,j]=NewtonR(Z[i,j])
    return Z

x=np.linspace(-1,1,512,dtype=complex)
xj=np.linspace(-1j,1j,512,dtype=complex)
xx,yy=np.meshgrid(x,xj)
X=xx+yy

Z=np.copy(X)
Z_cpu=np.copy(X)

threadsperblock = (16, 16)
blockspergrid_x = math.ceil(Z.shape[0] / threadsperblock[0])
blockspergrid_y = math.ceil(Z.shape[1] / threadsperblock[1])
blockspergrid = (blockspergrid_x, blockspergrid_y)

get_grid_cuda[blockspergrid, threadsperblock](Z)

Z_cpu=get_grid(Z_cpu)

np.testing.assert_allclose(Z_cpu, Z)

```

Hope this helps?

---

<div class="post-metadata">

**Author:** ![joej7](https://yyz2.discourse-cdn.com/free1/user_avatar/numba.discourse.group/joej7/32/514_2.png) [@joej7](https://numba.discourse.group/u/joej7)\
**Post date:** [February 13, 2022, 9:44pm UTC](https://numba.discourse.group/t/cuda-complex-powers/1158/3 "2022-02-13T21:44:36Z")

</div>

OMG, I’m stupid, lol. Thank you so much!