Révision 271
| ETSN/MyDFT_1.py (revision 271) | ||
|---|---|---|
| 3 | 3 |
import numpy as np |
| 4 | 4 |
import pyopencl as cl |
| 5 | 5 |
|
| 6 |
# piling 16 arithmetical functions |
|
| 7 |
def MySillyFunction(x): |
|
| 8 |
return(np.power(np.sqrt(np.log(np.exp(np.arctanh(np.tanh(np.arcsinh(np.sinh(np.arccosh(np.cosh(np.arctan(np.tan(np.arcsin(np.sin(np.arccos(np.cos(x))))))))))))))),2)) |
|
| 9 |
|
|
| 10 |
# Native Operation under Numpy (for prototyping & tests |
|
| 11 |
def NativeAddition(a_np,b_np): |
|
| 12 |
return(a_np+b_np) |
|
| 13 |
|
|
| 14 |
# Native Operation with MySillyFunction under Numpy (for prototyping & tests |
|
| 15 |
def NativeSillyAddition(a_np,b_np): |
|
| 16 |
return(MySillyFunction(a_np)+MySillyFunction(b_np)) |
|
| 17 |
|
|
| 18 | 6 |
# Naive Discrete Fourier Transform |
| 19 | 7 |
def MyDFT(x,y): |
| 20 | 8 |
from numpy import pi,cos,sin |
| ... | ... | |
| 27 | 15 |
Y[i]=Y[i]+x[j]*sin(2.*pi*i*j/size)+y[j]*cos(2.*pi*i*j/size) |
| 28 | 16 |
return(X,Y) |
| 29 | 17 |
|
| 30 |
# CUDA complete operation |
|
| 31 |
def CUDAAddition(a_np,b_np): |
|
| 32 |
import pycuda.autoinit |
|
| 33 |
import pycuda.driver as drv |
|
| 34 |
import numpy |
|
| 35 |
|
|
| 36 |
from pycuda.compiler import SourceModule |
|
| 37 |
mod = SourceModule("""
|
|
| 38 |
__global__ void sum(float *dest, float *a, float *b) |
|
| 39 |
{
|
|
| 40 |
// const int i = threadIdx.x; |
|
| 41 |
const int i = blockIdx.x; |
|
| 42 |
dest[i] = a[i] + b[i]; |
|
| 43 |
} |
|
| 44 |
""") |
|
| 45 |
|
|
| 46 |
# sum = mod.get_function("sum")
|
|
| 47 |
sum = mod.get_function("sum")
|
|
| 48 |
|
|
| 49 |
res_np = numpy.zeros_like(a_np) |
|
| 50 |
sum(drv.Out(res_np), drv.In(a_np), drv.In(b_np), |
|
| 51 |
block=(1,1,1), grid=(a_np.size,1)) |
|
| 52 |
return(res_np) |
|
| 53 |
|
|
| 54 |
# CUDA Silly complete operation |
|
| 55 |
def CUDASillyAddition(a_np,b_np): |
|
| 56 |
import pycuda.autoinit |
|
| 57 |
import pycuda.driver as drv |
|
| 58 |
import numpy |
|
| 59 |
|
|
| 60 |
from pycuda.compiler import SourceModule |
|
| 61 |
TimeIn=time.time() |
|
| 62 |
mod = SourceModule("""
|
|
| 63 |
__device__ float MySillyFunction(float x) |
|
| 64 |
{
|
|
| 65 |
return(pow(sqrt(log(exp(atanh(tanh(asinh(sinh(acosh(cosh(atan(tan(asin(sin(acos(cos(x))))))))))))))),2)); |
|
| 66 |
} |
|
| 67 |
|
|
| 68 |
__global__ void sillysum(float *dest, float *a, float *b) |
|
| 69 |
{
|
|
| 70 |
const int i = blockIdx.x; |
|
| 71 |
dest[i] = MySillyFunction(a[i]) + MySillyFunction(b[i]); |
|
| 72 |
} |
|
| 73 |
""") |
|
| 74 |
Elapsed=time.time()-TimeIn |
|
| 75 |
print("Definition of kernel : %.3f" % Elapsed)
|
|
| 76 |
|
|
| 77 |
TimeIn=time.time() |
|
| 78 |
# sum = mod.get_function("sum")
|
|
| 79 |
sillysum = mod.get_function("sillysum")
|
|
| 80 |
Elapsed=time.time()-TimeIn |
|
| 81 |
print("Synthesis of kernel : %.3f" % Elapsed)
|
|
| 82 |
|
|
| 83 |
TimeIn=time.time() |
|
| 84 |
res_np = numpy.zeros_like(a_np) |
|
| 85 |
Elapsed=time.time()-TimeIn |
|
| 86 |
print("Allocation on Host for results : %.3f" % Elapsed)
|
|
| 87 |
|
|
| 88 |
TimeIn=time.time() |
|
| 89 |
sillysum(drv.Out(res_np), drv.In(a_np), drv.In(b_np), |
|
| 90 |
block=(1,1,1), grid=(a_np.size,1)) |
|
| 91 |
Elapsed=time.time()-TimeIn |
|
| 92 |
print("Execution of kernel : %.3f" % Elapsed)
|
|
| 93 |
return(res_np) |
|
| 94 |
|
|
| 95 |
# OpenCL complete operation |
|
| 96 |
def OpenCLAddition(a_np,b_np): |
|
| 97 |
|
|
| 98 |
# Context creation |
|
| 99 |
ctx = cl.create_some_context() |
|
| 100 |
# Every process is stored in a queue |
|
| 101 |
queue = cl.CommandQueue(ctx) |
|
| 102 |
|
|
| 103 |
TimeIn=time.time() |
|
| 104 |
# Copy from Host to Device using pointers |
|
| 105 |
mf = cl.mem_flags |
|
| 106 |
a_g = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=a_np) |
|
| 107 |
b_g = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=b_np) |
|
| 108 |
Elapsed=time.time()-TimeIn |
|
| 109 |
print("Copy from Host 2 Device : %.3f" % Elapsed)
|
|
| 110 |
|
|
| 111 |
TimeIn=time.time() |
|
| 112 |
# Definition of kernel under OpenCL |
|
| 113 |
prg = cl.Program(ctx, """ |
|
| 114 |
__kernel void sum( |
|
| 115 |
__global const float *a_g, __global const float *b_g, __global float *res_g) |
|
| 116 |
{
|
|
| 117 |
int gid = get_global_id(0); |
|
| 118 |
res_g[gid] = a_g[gid] + b_g[gid]; |
|
| 119 |
} |
|
| 120 |
""").build() |
|
| 121 |
Elapsed=time.time()-TimeIn |
|
| 122 |
print("Building kernels : %.3f" % Elapsed)
|
|
| 123 |
|
|
| 124 |
TimeIn=time.time() |
|
| 125 |
# Memory allocation on Device for result |
|
| 126 |
res_g = cl.Buffer(ctx, mf.WRITE_ONLY, a_np.nbytes) |
|
| 127 |
Elapsed=time.time()-TimeIn |
|
| 128 |
print("Allocation on Device for results : %.3f" % Elapsed)
|
|
| 129 |
|
|
| 130 |
TimeIn=time.time() |
|
| 131 |
# Synthesis of function "sum" inside Kernel Sources |
|
| 132 |
knl = prg.sum # Use this Kernel object for repeated calls |
|
| 133 |
Elapsed=time.time()-TimeIn |
|
| 134 |
print("Synthesis of kernel : %.3f" % Elapsed)
|
|
| 135 |
|
|
| 136 |
TimeIn=time.time() |
|
| 137 |
# Call of kernel previously defined |
|
| 138 |
knl(queue, a_np.shape, None, a_g, b_g, res_g) |
|
| 139 |
Elapsed=time.time()-TimeIn |
|
| 140 |
print("Execution of kernel : %.3f" % Elapsed)
|
|
| 141 |
|
|
| 142 |
TimeIn=time.time() |
|
| 143 |
# Creation of vector for result with same size as input vectors |
|
| 144 |
res_np = np.empty_like(a_np) |
|
| 145 |
Elapsed=time.time()-TimeIn |
|
| 146 |
print("Allocation on Host for results: %.3f" % Elapsed)
|
|
| 147 |
|
|
| 148 |
TimeIn=time.time() |
|
| 149 |
# Copy from Device to Host |
|
| 150 |
cl.enqueue_copy(queue, res_np, res_g) |
|
| 151 |
Elapsed=time.time()-TimeIn |
|
| 152 |
print("Copy from Device 2 Host : %.3f" % Elapsed)
|
|
| 153 |
|
|
| 154 |
return(res_np) |
|
| 155 |
|
|
| 156 |
# OpenCL complete operation |
|
| 157 |
def OpenCLSillyAddition(a_np,b_np): |
|
| 158 |
|
|
| 159 |
# Context creation |
|
| 160 |
ctx = cl.create_some_context() |
|
| 161 |
# Every process is stored in a queue |
|
| 162 |
queue = cl.CommandQueue(ctx) |
|
| 163 |
|
|
| 164 |
TimeIn=time.time() |
|
| 165 |
# Copy from Host to Device using pointers |
|
| 166 |
mf = cl.mem_flags |
|
| 167 |
a_g = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=a_np) |
|
| 168 |
b_g = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=b_np) |
|
| 169 |
Elapsed=time.time()-TimeIn |
|
| 170 |
print("Copy from Host 2 Device : %.3f" % Elapsed)
|
|
| 171 |
|
|
| 172 |
TimeIn=time.time() |
|
| 173 |
# Definition of kernel under OpenCL |
|
| 174 |
prg = cl.Program(ctx, """ |
|
| 175 |
|
|
| 176 |
float MySillyFunction(float x) |
|
| 177 |
{
|
|
| 178 |
return(pow(sqrt(log(exp(atanh(tanh(asinh(sinh(acosh(cosh(atan(tan(asin(sin(acos(cos(x))))))))))))))),2)); |
|
| 179 |
} |
|
| 180 |
|
|
| 181 |
__kernel void sillysum( |
|
| 182 |
__global const float *a_g, __global const float *b_g, __global float *res_g) |
|
| 183 |
{
|
|
| 184 |
int gid = get_global_id(0); |
|
| 185 |
res_g[gid] = MySillyFunction(a_g[gid]) + MySillyFunction(b_g[gid]); |
|
| 186 |
} |
|
| 187 |
|
|
| 188 |
__kernel void sum( |
|
| 189 |
__global const float *a_g, __global const float *b_g, __global float *res_g) |
|
| 190 |
{
|
|
| 191 |
int gid = get_global_id(0); |
|
| 192 |
res_g[gid] = a_g[gid] + b_g[gid]; |
|
| 193 |
} |
|
| 194 |
""").build() |
|
| 195 |
Elapsed=time.time()-TimeIn |
|
| 196 |
print("Building kernels : %.3f" % Elapsed)
|
|
| 197 |
|
|
| 198 |
TimeIn=time.time() |
|
| 199 |
# Memory allocation on Device for result |
|
| 200 |
res_g = cl.Buffer(ctx, mf.WRITE_ONLY, a_np.nbytes) |
|
| 201 |
Elapsed=time.time()-TimeIn |
|
| 202 |
print("Allocation on Device for results : %.3f" % Elapsed)
|
|
| 203 |
|
|
| 204 |
TimeIn=time.time() |
|
| 205 |
# Synthesis of function "sillysum" inside Kernel Sources |
|
| 206 |
knl = prg.sillysum # Use this Kernel object for repeated calls |
|
| 207 |
Elapsed=time.time()-TimeIn |
|
| 208 |
print("Synthesis of kernel : %.3f" % Elapsed)
|
|
| 209 |
|
|
| 210 |
TimeIn=time.time() |
|
| 211 |
# Call of kernel previously defined |
|
| 212 |
CallCL=knl(queue, a_np.shape, None, a_g, b_g, res_g) |
|
| 213 |
# |
|
| 214 |
CallCL.wait() |
|
| 215 |
Elapsed=time.time()-TimeIn |
|
| 216 |
print("Execution of kernel : %.3f" % Elapsed)
|
|
| 217 |
|
|
| 218 |
TimeIn=time.time() |
|
| 219 |
# Creation of vector for result with same size as input vectors |
|
| 220 |
res_np = np.empty_like(a_np) |
|
| 221 |
Elapsed=time.time()-TimeIn |
|
| 222 |
print("Allocation on Host for results: %.3f" % Elapsed)
|
|
| 223 |
|
|
| 224 |
TimeIn=time.time() |
|
| 225 |
# Copy from Device to Host |
|
| 226 |
cl.enqueue_copy(queue, res_np, res_g) |
|
| 227 |
Elapsed=time.time()-TimeIn |
|
| 228 |
print("Copy from Device 2 Host : %.3f" % Elapsed)
|
|
| 229 |
|
|
| 230 |
return(res_np) |
|
| 231 |
|
|
| 232 | 18 |
import sys |
| 233 | 19 |
import time |
| 234 | 20 |
|
| ... | ... | |
| 240 | 26 |
SIZE=int(sys.argv[1]) |
| 241 | 27 |
print("Size of vectors set to %i" % SIZE)
|
| 242 | 28 |
except: |
| 243 |
SIZE=50000
|
|
| 29 |
SIZE=256
|
|
| 244 | 30 |
print("Size of vectors set to default size %i" % SIZE)
|
| 245 | 31 |
|
| 246 |
# a_np = np.random.rand(SIZE).astype(np.float32) |
|
| 247 |
# b_np = np.random.rand(SIZE).astype(np.float32) |
|
| 248 |
|
|
| 249 | 32 |
a_np = np.ones(SIZE).astype(np.float32) |
| 250 | 33 |
b_np = np.ones(SIZE).astype(np.float32) |
| 251 | 34 |
|
| 35 |
C_np = np.zeros(SIZE).astype(np.float32) |
|
| 36 |
D_np = np.zeros(SIZE).astype(np.float32) |
|
| 37 |
C_np[0] = np.float32(SIZE) |
|
| 38 |
D_np[0] = np.float32(SIZE) |
|
| 39 |
|
|
| 252 | 40 |
# Native Implementation |
| 253 | 41 |
TimeIn=time.time() |
| 254 |
# res_np=NativeAddition(a_np,b_np) |
|
| 255 |
# res_np=NativeSillyAddition(a_np,b_np) |
|
| 256 | 42 |
c_np,d_np=MyDFT(a_np,b_np) |
| 257 | 43 |
NativeElapsed=time.time()-TimeIn |
| 258 | 44 |
NativeRate=int(SIZE/NativeElapsed) |
| 259 | 45 |
print("NativeRate: %i" % NativeRate)
|
| 260 |
|
|
| 261 |
print(c_np,d_np) |
|
| 262 |
|
|
| 263 |
# # OpenCL Implementation |
|
| 264 |
# TimeIn=time.time() |
|
| 265 |
# # res_cl=OpenCLAddition(a_np,b_np) |
|
| 266 |
# res_cl=OpenCLSillyAddition(a_np,b_np) |
|
| 267 |
# OpenCLElapsed=time.time()-TimeIn |
|
| 268 |
# OpenCLRate=int(SIZE/OpenCLElapsed) |
|
| 269 |
# print("OpenCLRate: %i" % OpenCLRate)
|
|
| 270 |
|
|
| 271 |
# # CUDA Implementation |
|
| 272 |
# TimeIn=time.time() |
|
| 273 |
# # res_cuda=CUDAAddition(a_np,b_np) |
|
| 274 |
# res_cuda=CUDASillyAddition(a_np,b_np) |
|
| 275 |
# CUDAElapsed=time.time()-TimeIn |
|
| 276 |
# CUDARate=int(SIZE/CUDAElapsed) |
|
| 277 |
# print("CUDARate: %i" % CUDARate)
|
|
| 278 |
|
|
| 279 |
# print("OpenCLvsNative ratio: %f" % (OpenCLRate/NativeRate))
|
|
| 280 |
# print("CUDAvsNative ratio: %f" % (CUDARate/NativeRate))
|
|
| 281 |
|
|
| 282 |
# # Check on OpenCL with Numpy: |
|
| 283 |
# print(res_cl - res_np) |
|
| 284 |
# print(np.linalg.norm(res_cl - res_np)) |
|
| 285 |
# try: |
|
| 286 |
# assert np.allclose(res_np, res_cl) |
|
| 287 |
# except: |
|
| 288 |
# print("Results between Native & OpenCL seem to be too different!")
|
|
| 289 |
|
|
| 290 |
# # Check on CUDA with Numpy: |
|
| 291 |
# print(res_cuda - res_np) |
|
| 292 |
# print(np.linalg.norm(res_cuda - res_np)) |
|
| 293 |
# try: |
|
| 294 |
# assert np.allclose(res_np, res_cuda) |
|
| 295 |
# except: |
|
| 296 |
# print("Results between Native & CUDA seem to be too different!")
|
|
| 297 |
|
|
| 298 |
|
|
| 46 |
print("Precision: ",np.linalg.norm(c_np-C_np),np.linalg.norm(d_np-D_np))
|
|
| ETSN/MyDFT_2.py (revision 271) | ||
|---|---|---|
| 3 | 3 |
import numpy as np |
| 4 | 4 |
import pyopencl as cl |
| 5 | 5 |
|
| 6 |
# piling 16 arithmetical functions |
|
| 7 |
def MySillyFunction(x): |
|
| 8 |
return(np.power(np.sqrt(np.log(np.exp(np.arctanh(np.tanh(np.arcsinh(np.sinh(np.arccosh(np.cosh(np.arctan(np.tan(np.arcsin(np.sin(np.arccos(np.cos(x))))))))))))))),2)) |
|
| 9 |
|
|
| 10 |
# Native Operation under Numpy (for prototyping & tests |
|
| 11 |
def NativeAddition(a_np,b_np): |
|
| 12 |
return(a_np+b_np) |
|
| 13 |
|
|
| 14 |
# Native Operation with MySillyFunction under Numpy (for prototyping & tests |
|
| 15 |
def NativeSillyAddition(a_np,b_np): |
|
| 16 |
return(MySillyFunction(a_np)+MySillyFunction(b_np)) |
|
| 17 |
|
|
| 18 | 6 |
# Naive Discrete Fourier Transform |
| 19 | 7 |
def MyDFT(x,y): |
| 20 | 8 |
from numpy import pi,cos,sin |
| ... | ... | |
| 39 | 27 |
Y[i]=np.sum(np.add(np.multiply(np.sin(i*nj),x),np.multiply(np.cos(i*nj),y))) |
| 40 | 28 |
return(X,Y) |
| 41 | 29 |
|
| 42 |
# CUDA complete operation |
|
| 43 |
def CUDAAddition(a_np,b_np): |
|
| 44 |
import pycuda.autoinit |
|
| 45 |
import pycuda.driver as drv |
|
| 46 |
import numpy |
|
| 47 |
|
|
| 48 |
from pycuda.compiler import SourceModule |
|
| 49 |
mod = SourceModule("""
|
|
| 50 |
__global__ void sum(float *dest, float *a, float *b) |
|
| 51 |
{
|
|
| 52 |
// const int i = threadIdx.x; |
|
| 53 |
const int i = blockIdx.x; |
|
| 54 |
dest[i] = a[i] + b[i]; |
|
| 55 |
} |
|
| 56 |
""") |
|
| 57 |
|
|
| 58 |
# sum = mod.get_function("sum")
|
|
| 59 |
sum = mod.get_function("sum")
|
|
| 60 |
|
|
| 61 |
res_np = numpy.zeros_like(a_np) |
|
| 62 |
sum(drv.Out(res_np), drv.In(a_np), drv.In(b_np), |
|
| 63 |
block=(1,1,1), grid=(a_np.size,1)) |
|
| 64 |
return(res_np) |
|
| 65 |
|
|
| 66 |
# CUDA Silly complete operation |
|
| 67 |
def CUDASillyAddition(a_np,b_np): |
|
| 68 |
import pycuda.autoinit |
|
| 69 |
import pycuda.driver as drv |
|
| 70 |
import numpy |
|
| 71 |
|
|
| 72 |
from pycuda.compiler import SourceModule |
|
| 73 |
TimeIn=time.time() |
|
| 74 |
mod = SourceModule("""
|
|
| 75 |
__device__ float MySillyFunction(float x) |
|
| 76 |
{
|
|
| 77 |
return(pow(sqrt(log(exp(atanh(tanh(asinh(sinh(acosh(cosh(atan(tan(asin(sin(acos(cos(x))))))))))))))),2)); |
|
| 78 |
} |
|
| 79 |
|
|
| 80 |
__global__ void sillysum(float *dest, float *a, float *b) |
|
| 81 |
{
|
|
| 82 |
const int i = blockIdx.x; |
|
| 83 |
dest[i] = MySillyFunction(a[i]) + MySillyFunction(b[i]); |
|
| 84 |
} |
|
| 85 |
""") |
|
| 86 |
Elapsed=time.time()-TimeIn |
|
| 87 |
print("Definition of kernel : %.3f" % Elapsed)
|
|
| 88 |
|
|
| 89 |
TimeIn=time.time() |
|
| 90 |
# sum = mod.get_function("sum")
|
|
| 91 |
sillysum = mod.get_function("sillysum")
|
|
| 92 |
Elapsed=time.time()-TimeIn |
|
| 93 |
print("Synthesis of kernel : %.3f" % Elapsed)
|
|
| 94 |
|
|
| 95 |
TimeIn=time.time() |
|
| 96 |
res_np = numpy.zeros_like(a_np) |
|
| 97 |
Elapsed=time.time()-TimeIn |
|
| 98 |
print("Allocation on Host for results : %.3f" % Elapsed)
|
|
| 99 |
|
|
| 100 |
TimeIn=time.time() |
|
| 101 |
sillysum(drv.Out(res_np), drv.In(a_np), drv.In(b_np), |
|
| 102 |
block=(1,1,1), grid=(a_np.size,1)) |
|
| 103 |
Elapsed=time.time()-TimeIn |
|
| 104 |
print("Execution of kernel : %.3f" % Elapsed)
|
|
| 105 |
return(res_np) |
|
| 106 |
|
|
| 107 |
# OpenCL complete operation |
|
| 108 |
def OpenCLAddition(a_np,b_np): |
|
| 109 |
|
|
| 110 |
# Context creation |
|
| 111 |
ctx = cl.create_some_context() |
|
| 112 |
# Every process is stored in a queue |
|
| 113 |
queue = cl.CommandQueue(ctx) |
|
| 114 |
|
|
| 115 |
TimeIn=time.time() |
|
| 116 |
# Copy from Host to Device using pointers |
|
| 117 |
mf = cl.mem_flags |
|
| 118 |
a_g = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=a_np) |
|
| 119 |
b_g = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=b_np) |
|
| 120 |
Elapsed=time.time()-TimeIn |
|
| 121 |
print("Copy from Host 2 Device : %.3f" % Elapsed)
|
|
| 122 |
|
|
| 123 |
TimeIn=time.time() |
|
| 124 |
# Definition of kernel under OpenCL |
|
| 125 |
prg = cl.Program(ctx, """ |
|
| 126 |
__kernel void sum( |
|
| 127 |
__global const float *a_g, __global const float *b_g, __global float *res_g) |
|
| 128 |
{
|
|
| 129 |
int gid = get_global_id(0); |
|
| 130 |
res_g[gid] = a_g[gid] + b_g[gid]; |
|
| 131 |
} |
|
| 132 |
""").build() |
|
| 133 |
Elapsed=time.time()-TimeIn |
|
| 134 |
print("Building kernels : %.3f" % Elapsed)
|
|
| 135 |
|
|
| 136 |
TimeIn=time.time() |
|
| 137 |
# Memory allocation on Device for result |
|
| 138 |
res_g = cl.Buffer(ctx, mf.WRITE_ONLY, a_np.nbytes) |
|
| 139 |
Elapsed=time.time()-TimeIn |
|
| 140 |
print("Allocation on Device for results : %.3f" % Elapsed)
|
|
| 141 |
|
|
| 142 |
TimeIn=time.time() |
|
| 143 |
# Synthesis of function "sum" inside Kernel Sources |
|
| 144 |
knl = prg.sum # Use this Kernel object for repeated calls |
|
| 145 |
Elapsed=time.time()-TimeIn |
|
| 146 |
print("Synthesis of kernel : %.3f" % Elapsed)
|
|
| 147 |
|
|
| 148 |
TimeIn=time.time() |
|
| 149 |
# Call of kernel previously defined |
|
| 150 |
knl(queue, a_np.shape, None, a_g, b_g, res_g) |
|
| 151 |
Elapsed=time.time()-TimeIn |
|
| 152 |
print("Execution of kernel : %.3f" % Elapsed)
|
|
| 153 |
|
|
| 154 |
TimeIn=time.time() |
|
| 155 |
# Creation of vector for result with same size as input vectors |
|
| 156 |
res_np = np.empty_like(a_np) |
|
| 157 |
Elapsed=time.time()-TimeIn |
|
| 158 |
print("Allocation on Host for results: %.3f" % Elapsed)
|
|
| 159 |
|
|
| 160 |
TimeIn=time.time() |
|
| 161 |
# Copy from Device to Host |
|
| 162 |
cl.enqueue_copy(queue, res_np, res_g) |
|
| 163 |
Elapsed=time.time()-TimeIn |
|
| 164 |
print("Copy from Device 2 Host : %.3f" % Elapsed)
|
|
| 165 |
|
|
| 166 |
return(res_np) |
|
| 167 |
|
|
| 168 |
# OpenCL complete operation |
|
| 169 |
def OpenCLSillyAddition(a_np,b_np): |
|
| 170 |
|
|
| 171 |
# Context creation |
|
| 172 |
ctx = cl.create_some_context() |
|
| 173 |
# Every process is stored in a queue |
|
| 174 |
queue = cl.CommandQueue(ctx) |
|
| 175 |
|
|
| 176 |
TimeIn=time.time() |
|
| 177 |
# Copy from Host to Device using pointers |
|
| 178 |
mf = cl.mem_flags |
|
| 179 |
a_g = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=a_np) |
|
| 180 |
b_g = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=b_np) |
|
| 181 |
Elapsed=time.time()-TimeIn |
|
| 182 |
print("Copy from Host 2 Device : %.3f" % Elapsed)
|
|
| 183 |
|
|
| 184 |
TimeIn=time.time() |
|
| 185 |
# Definition of kernel under OpenCL |
|
| 186 |
prg = cl.Program(ctx, """ |
|
| 187 |
|
|
| 188 |
float MySillyFunction(float x) |
|
| 189 |
{
|
|
| 190 |
return(pow(sqrt(log(exp(atanh(tanh(asinh(sinh(acosh(cosh(atan(tan(asin(sin(acos(cos(x))))))))))))))),2)); |
|
| 191 |
} |
|
| 192 |
|
|
| 193 |
__kernel void sillysum( |
|
| 194 |
__global const float *a_g, __global const float *b_g, __global float *res_g) |
|
| 195 |
{
|
|
| 196 |
int gid = get_global_id(0); |
|
| 197 |
res_g[gid] = MySillyFunction(a_g[gid]) + MySillyFunction(b_g[gid]); |
|
| 198 |
} |
|
| 199 |
|
|
| 200 |
__kernel void sum( |
|
| 201 |
__global const float *a_g, __global const float *b_g, __global float *res_g) |
|
| 202 |
{
|
|
| 203 |
int gid = get_global_id(0); |
|
| 204 |
res_g[gid] = a_g[gid] + b_g[gid]; |
|
| 205 |
} |
|
| 206 |
""").build() |
|
| 207 |
Elapsed=time.time()-TimeIn |
|
| 208 |
print("Building kernels : %.3f" % Elapsed)
|
|
| 209 |
|
|
| 210 |
TimeIn=time.time() |
|
| 211 |
# Memory allocation on Device for result |
|
| 212 |
res_g = cl.Buffer(ctx, mf.WRITE_ONLY, a_np.nbytes) |
|
| 213 |
Elapsed=time.time()-TimeIn |
|
| 214 |
print("Allocation on Device for results : %.3f" % Elapsed)
|
|
| 215 |
|
|
| 216 |
TimeIn=time.time() |
|
| 217 |
# Synthesis of function "sillysum" inside Kernel Sources |
|
| 218 |
knl = prg.sillysum # Use this Kernel object for repeated calls |
|
| 219 |
Elapsed=time.time()-TimeIn |
|
| 220 |
print("Synthesis of kernel : %.3f" % Elapsed)
|
|
| 221 |
|
|
| 222 |
TimeIn=time.time() |
|
| 223 |
# Call of kernel previously defined |
|
| 224 |
CallCL=knl(queue, a_np.shape, None, a_g, b_g, res_g) |
|
| 225 |
# |
|
| 226 |
CallCL.wait() |
|
| 227 |
Elapsed=time.time()-TimeIn |
|
| 228 |
print("Execution of kernel : %.3f" % Elapsed)
|
|
| 229 |
|
|
| 230 |
TimeIn=time.time() |
|
| 231 |
# Creation of vector for result with same size as input vectors |
|
| 232 |
res_np = np.empty_like(a_np) |
|
| 233 |
Elapsed=time.time()-TimeIn |
|
| 234 |
print("Allocation on Host for results: %.3f" % Elapsed)
|
|
| 235 |
|
|
| 236 |
TimeIn=time.time() |
|
| 237 |
# Copy from Device to Host |
|
| 238 |
cl.enqueue_copy(queue, res_np, res_g) |
|
| 239 |
Elapsed=time.time()-TimeIn |
|
| 240 |
print("Copy from Device 2 Host : %.3f" % Elapsed)
|
|
| 241 |
|
|
| 242 |
return(res_np) |
|
| 243 |
|
|
| 244 | 30 |
import sys |
| 245 | 31 |
import time |
| 246 | 32 |
|
| ... | ... | |
| 252 | 38 |
SIZE=int(sys.argv[1]) |
| 253 | 39 |
print("Size of vectors set to %i" % SIZE)
|
| 254 | 40 |
except: |
| 255 |
SIZE=50000
|
|
| 41 |
SIZE=256
|
|
| 256 | 42 |
print("Size of vectors set to default size %i" % SIZE)
|
| 257 | 43 |
|
| 258 |
# a_np = np.random.rand(SIZE).astype(np.float32) |
|
| 259 |
# b_np = np.random.rand(SIZE).astype(np.float32) |
|
| 260 |
|
|
| 261 | 44 |
a_np = np.ones(SIZE).astype(np.float32) |
| 262 | 45 |
b_np = np.ones(SIZE).astype(np.float32) |
| 263 | 46 |
|
| 47 |
C_np = np.zeros(SIZE).astype(np.float32) |
|
| 48 |
D_np = np.zeros(SIZE).astype(np.float32) |
|
| 49 |
C_np[0] = np.float32(SIZE) |
|
| 50 |
D_np[0] = np.float32(SIZE) |
|
| 51 |
|
|
| 264 | 52 |
# Native & Naive Implementation |
| 265 | 53 |
print("Performing naive implementation")
|
| 266 | 54 |
TimeIn=time.time() |
| ... | ... | |
| 268 | 56 |
NativeElapsed=time.time()-TimeIn |
| 269 | 57 |
NativeRate=int(SIZE/NativeElapsed) |
| 270 | 58 |
print("NativeRate: %i" % NativeRate)
|
| 271 |
print(c_np,d_np)
|
|
| 59 |
print("Precision: ",np.linalg.norm(c_np-C_np),np.linalg.norm(d_np-D_np))
|
|
| 272 | 60 |
|
| 273 | 61 |
# Native & Numpy Implementation |
| 274 | 62 |
print("Performing Numpy implementation")
|
| ... | ... | |
| 277 | 65 |
NumpyElapsed=time.time()-TimeIn |
| 278 | 66 |
NumpyRate=int(SIZE/NumpyElapsed) |
| 279 | 67 |
print("NumpyRate: %i" % NumpyRate)
|
| 280 |
print(e_np,f_np) |
|
| 281 |
|
|
| 282 |
print(np.linalg.norm(c_np-e_np)) |
|
| 283 |
print(np.linalg.norm(d_np-f_np)) |
|
| 284 |
|
|
| 285 |
# # OpenCL Implementation |
|
| 286 |
# TimeIn=time.time() |
|
| 287 |
# # res_cl=OpenCLAddition(a_np,b_np) |
|
| 288 |
# res_cl=OpenCLSillyAddition(a_np,b_np) |
|
| 289 |
# OpenCLElapsed=time.time()-TimeIn |
|
| 290 |
# OpenCLRate=int(SIZE/OpenCLElapsed) |
|
| 291 |
# print("OpenCLRate: %i" % OpenCLRate)
|
|
| 292 |
|
|
| 293 |
# # CUDA Implementation |
|
| 294 |
# TimeIn=time.time() |
|
| 295 |
# # res_cuda=CUDAAddition(a_np,b_np) |
|
| 296 |
# res_cuda=CUDASillyAddition(a_np,b_np) |
|
| 297 |
# CUDAElapsed=time.time()-TimeIn |
|
| 298 |
# CUDARate=int(SIZE/CUDAElapsed) |
|
| 299 |
# print("CUDARate: %i" % CUDARate)
|
|
| 300 |
|
|
| 301 |
# print("OpenCLvsNative ratio: %f" % (OpenCLRate/NativeRate))
|
|
| 302 |
# print("CUDAvsNative ratio: %f" % (CUDARate/NativeRate))
|
|
| 303 |
|
|
| 304 |
# # Check on OpenCL with Numpy: |
|
| 305 |
# print(res_cl - res_np) |
|
| 306 |
# print(np.linalg.norm(res_cl - res_np)) |
|
| 307 |
# try: |
|
| 308 |
# assert np.allclose(res_np, res_cl) |
|
| 309 |
# except: |
|
| 310 |
# print("Results between Native & OpenCL seem to be too different!")
|
|
| 311 |
|
|
| 312 |
# # Check on CUDA with Numpy: |
|
| 313 |
# print(res_cuda - res_np) |
|
| 314 |
# print(np.linalg.norm(res_cuda - res_np)) |
|
| 315 |
# try: |
|
| 316 |
# assert np.allclose(res_np, res_cuda) |
|
| 317 |
# except: |
|
| 318 |
# print("Results between Native & CUDA seem to be too different!")
|
|
| 319 |
|
|
| 320 |
|
|
| 68 |
print("Precision: ",np.linalg.norm(e_np-C_np),np.linalg.norm(f_np-D_np))
|
|
| ETSN/MyDFT_3.py (revision 271) | ||
|---|---|---|
| 4 | 4 |
import pyopencl as cl |
| 5 | 5 |
from numpy import pi,cos,sin |
| 6 | 6 |
|
| 7 |
# piling 16 arithmetical functions |
|
| 8 |
def MySillyFunction(x): |
|
| 9 |
return(np.power(np.sqrt(np.log(np.exp(np.arctanh(np.tanh(np.arcsinh(np.sinh(np.arccosh(np.cosh(np.arctan(np.tan(np.arcsin(np.sin(np.arccos(np.cos(x))))))))))))))),2)) |
|
| 10 |
|
|
| 11 |
# Native Operation under Numpy (for prototyping & tests |
|
| 12 |
def NativeAddition(a_np,b_np): |
|
| 13 |
return(a_np+b_np) |
|
| 14 |
|
|
| 15 |
# Native Operation with MySillyFunction under Numpy (for prototyping & tests |
|
| 16 |
def NativeSillyAddition(a_np,b_np): |
|
| 17 |
return(MySillyFunction(a_np)+MySillyFunction(b_np)) |
|
| 18 |
|
|
| 19 | 7 |
# Naive Discrete Fourier Transform |
| 20 | 8 |
def MyDFT(x,y): |
| 21 |
from numpy import pi,cos,sin |
|
| 22 | 9 |
size=x.shape[0] |
| 23 | 10 |
X=np.zeros(size).astype(np.float32) |
| 24 | 11 |
Y=np.zeros(size).astype(np.float32) |
| ... | ... | |
| 44 | 31 |
@numba.njit(parallel=True) |
| 45 | 32 |
def NumbaDFT(x,y): |
| 46 | 33 |
size=x.shape[0] |
| 47 |
X=np.zeros(size) |
|
| 48 |
Y=np.zeros(size) |
|
| 34 |
X=np.zeros(size).astype(np.float32)
|
|
| 35 |
Y=np.zeros(size).astype(np.float32)
|
|
| 49 | 36 |
nj=np.multiply(2.0*np.pi/size,np.arange(size)).astype(np.float32) |
| 50 | 37 |
for i in numba.prange(size): |
| 51 | 38 |
X[i]=np.sum(np.subtract(np.multiply(np.cos(i*nj),x),np.multiply(np.sin(i*nj),y))) |
| 52 | 39 |
Y[i]=np.sum(np.add(np.multiply(np.sin(i*nj),x),np.multiply(np.cos(i*nj),y))) |
| 53 | 40 |
return(X,Y) |
| 54 | 41 |
|
| 55 |
# CUDA complete operation |
|
| 56 |
def CUDAAddition(a_np,b_np): |
|
| 57 |
import pycuda.autoinit |
|
| 58 |
import pycuda.driver as drv |
|
| 59 |
import numpy |
|
| 60 |
|
|
| 61 |
from pycuda.compiler import SourceModule |
|
| 62 |
mod = SourceModule("""
|
|
| 63 |
__global__ void sum(float *dest, float *a, float *b) |
|
| 64 |
{
|
|
| 65 |
// const int i = threadIdx.x; |
|
| 66 |
const int i = blockIdx.x; |
|
| 67 |
dest[i] = a[i] + b[i]; |
|
| 68 |
} |
|
| 69 |
""") |
|
| 70 |
|
|
| 71 |
# sum = mod.get_function("sum")
|
|
| 72 |
sum = mod.get_function("sum")
|
|
| 73 |
|
|
| 74 |
res_np = numpy.zeros_like(a_np) |
|
| 75 |
sum(drv.Out(res_np), drv.In(a_np), drv.In(b_np), |
|
| 76 |
block=(1,1,1), grid=(a_np.size,1)) |
|
| 77 |
return(res_np) |
|
| 78 |
|
|
| 79 |
# CUDA Silly complete operation |
|
| 80 |
def CUDASillyAddition(a_np,b_np): |
|
| 81 |
import pycuda.autoinit |
|
| 82 |
import pycuda.driver as drv |
|
| 83 |
import numpy |
|
| 84 |
|
|
| 85 |
from pycuda.compiler import SourceModule |
|
| 86 |
TimeIn=time.time() |
|
| 87 |
mod = SourceModule("""
|
|
| 88 |
__device__ float MySillyFunction(float x) |
|
| 89 |
{
|
|
| 90 |
return(pow(sqrt(log(exp(atanh(tanh(asinh(sinh(acosh(cosh(atan(tan(asin(sin(acos(cos(x))))))))))))))),2)); |
|
| 91 |
} |
|
| 92 |
|
|
| 93 |
__global__ void sillysum(float *dest, float *a, float *b) |
|
| 94 |
{
|
|
| 95 |
const int i = blockIdx.x; |
|
| 96 |
dest[i] = MySillyFunction(a[i]) + MySillyFunction(b[i]); |
|
| 97 |
} |
|
| 98 |
""") |
|
| 99 |
Elapsed=time.time()-TimeIn |
|
| 100 |
print("Definition of kernel : %.3f" % Elapsed)
|
|
| 101 |
|
|
| 102 |
TimeIn=time.time() |
|
| 103 |
# sum = mod.get_function("sum")
|
|
| 104 |
sillysum = mod.get_function("sillysum")
|
|
| 105 |
Elapsed=time.time()-TimeIn |
|
| 106 |
print("Synthesis of kernel : %.3f" % Elapsed)
|
|
| 107 |
|
|
| 108 |
TimeIn=time.time() |
|
| 109 |
res_np = numpy.zeros_like(a_np) |
|
| 110 |
Elapsed=time.time()-TimeIn |
|
| 111 |
print("Allocation on Host for results : %.3f" % Elapsed)
|
|
| 112 |
|
|
| 113 |
TimeIn=time.time() |
|
| 114 |
sillysum(drv.Out(res_np), drv.In(a_np), drv.In(b_np), |
|
| 115 |
block=(1,1,1), grid=(a_np.size,1)) |
|
| 116 |
Elapsed=time.time()-TimeIn |
|
| 117 |
print("Execution of kernel : %.3f" % Elapsed)
|
|
| 118 |
return(res_np) |
|
| 119 |
|
|
| 120 |
# OpenCL complete operation |
|
| 121 |
def OpenCLAddition(a_np,b_np): |
|
| 122 |
|
|
| 123 |
# Context creation |
|
| 124 |
ctx = cl.create_some_context() |
|
| 125 |
# Every process is stored in a queue |
|
| 126 |
queue = cl.CommandQueue(ctx) |
|
| 127 |
|
|
| 128 |
TimeIn=time.time() |
|
| 129 |
# Copy from Host to Device using pointers |
|
| 130 |
mf = cl.mem_flags |
|
| 131 |
a_g = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=a_np) |
|
| 132 |
b_g = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=b_np) |
|
| 133 |
Elapsed=time.time()-TimeIn |
|
| 134 |
print("Copy from Host 2 Device : %.3f" % Elapsed)
|
|
| 135 |
|
|
| 136 |
TimeIn=time.time() |
|
| 137 |
# Definition of kernel under OpenCL |
|
| 138 |
prg = cl.Program(ctx, """ |
|
| 139 |
__kernel void sum( |
|
| 140 |
__global const float *a_g, __global const float *b_g, __global float *res_g) |
|
| 141 |
{
|
|
| 142 |
int gid = get_global_id(0); |
|
| 143 |
res_g[gid] = a_g[gid] + b_g[gid]; |
|
| 144 |
} |
|
| 145 |
""").build() |
|
| 146 |
Elapsed=time.time()-TimeIn |
|
| 147 |
print("Building kernels : %.3f" % Elapsed)
|
|
| 148 |
|
|
| 149 |
TimeIn=time.time() |
|
| 150 |
# Memory allocation on Device for result |
|
| 151 |
res_g = cl.Buffer(ctx, mf.WRITE_ONLY, a_np.nbytes) |
|
| 152 |
Elapsed=time.time()-TimeIn |
|
| 153 |
print("Allocation on Device for results : %.3f" % Elapsed)
|
|
| 154 |
|
|
| 155 |
TimeIn=time.time() |
|
| 156 |
# Synthesis of function "sum" inside Kernel Sources |
|
| 157 |
knl = prg.sum # Use this Kernel object for repeated calls |
|
| 158 |
Elapsed=time.time()-TimeIn |
|
| 159 |
print("Synthesis of kernel : %.3f" % Elapsed)
|
|
| 160 |
|
|
| 161 |
TimeIn=time.time() |
|
| 162 |
# Call of kernel previously defined |
|
| 163 |
knl(queue, a_np.shape, None, a_g, b_g, res_g) |
|
| 164 |
Elapsed=time.time()-TimeIn |
|
| 165 |
print("Execution of kernel : %.3f" % Elapsed)
|
|
| 166 |
|
|
| 167 |
TimeIn=time.time() |
|
| 168 |
# Creation of vector for result with same size as input vectors |
|
| 169 |
res_np = np.empty_like(a_np) |
|
| 170 |
Elapsed=time.time()-TimeIn |
|
| 171 |
print("Allocation on Host for results: %.3f" % Elapsed)
|
|
| 172 |
|
|
| 173 |
TimeIn=time.time() |
|
| 174 |
# Copy from Device to Host |
|
| 175 |
cl.enqueue_copy(queue, res_np, res_g) |
|
| 176 |
Elapsed=time.time()-TimeIn |
|
| 177 |
print("Copy from Device 2 Host : %.3f" % Elapsed)
|
|
| 178 |
|
|
| 179 |
return(res_np) |
|
| 180 |
|
|
| 181 |
# OpenCL complete operation |
|
| 182 |
def OpenCLSillyAddition(a_np,b_np): |
|
| 183 |
|
|
| 184 |
# Context creation |
|
| 185 |
ctx = cl.create_some_context() |
|
| 186 |
# Every process is stored in a queue |
|
| 187 |
queue = cl.CommandQueue(ctx) |
|
| 188 |
|
|
| 189 |
TimeIn=time.time() |
|
| 190 |
# Copy from Host to Device using pointers |
|
| 191 |
mf = cl.mem_flags |
|
| 192 |
a_g = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=a_np) |
|
| 193 |
b_g = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=b_np) |
|
| 194 |
Elapsed=time.time()-TimeIn |
|
| 195 |
print("Copy from Host 2 Device : %.3f" % Elapsed)
|
|
| 196 |
|
|
| 197 |
TimeIn=time.time() |
|
| 198 |
# Definition of kernel under OpenCL |
|
| 199 |
prg = cl.Program(ctx, """ |
|
| 200 |
|
|
| 201 |
float MySillyFunction(float x) |
|
| 202 |
{
|
|
| 203 |
return(pow(sqrt(log(exp(atanh(tanh(asinh(sinh(acosh(cosh(atan(tan(asin(sin(acos(cos(x))))))))))))))),2)); |
|
| 204 |
} |
|
| 205 |
|
|
| 206 |
__kernel void sillysum( |
|
| 207 |
__global const float *a_g, __global const float *b_g, __global float *res_g) |
|
| 208 |
{
|
|
| 209 |
int gid = get_global_id(0); |
|
| 210 |
res_g[gid] = MySillyFunction(a_g[gid]) + MySillyFunction(b_g[gid]); |
|
| 211 |
} |
|
| 212 |
|
|
| 213 |
__kernel void sum( |
|
| 214 |
__global const float *a_g, __global const float *b_g, __global float *res_g) |
|
| 215 |
{
|
|
| 216 |
int gid = get_global_id(0); |
|
| 217 |
res_g[gid] = a_g[gid] + b_g[gid]; |
|
| 218 |
} |
|
| 219 |
""").build() |
|
| 220 |
Elapsed=time.time()-TimeIn |
|
| 221 |
print("Building kernels : %.3f" % Elapsed)
|
|
| 222 |
|
|
| 223 |
TimeIn=time.time() |
|
| 224 |
# Memory allocation on Device for result |
|
| 225 |
res_g = cl.Buffer(ctx, mf.WRITE_ONLY, a_np.nbytes) |
|
| 226 |
Elapsed=time.time()-TimeIn |
|
| 227 |
print("Allocation on Device for results : %.3f" % Elapsed)
|
|
| 228 |
|
|
| 229 |
TimeIn=time.time() |
|
| 230 |
# Synthesis of function "sillysum" inside Kernel Sources |
|
| 231 |
knl = prg.sillysum # Use this Kernel object for repeated calls |
|
| 232 |
Elapsed=time.time()-TimeIn |
|
| 233 |
print("Synthesis of kernel : %.3f" % Elapsed)
|
|
| 234 |
|
|
| 235 |
TimeIn=time.time() |
|
| 236 |
# Call of kernel previously defined |
|
| 237 |
CallCL=knl(queue, a_np.shape, None, a_g, b_g, res_g) |
|
| 238 |
# |
|
| 239 |
CallCL.wait() |
|
| 240 |
Elapsed=time.time()-TimeIn |
|
| 241 |
print("Execution of kernel : %.3f" % Elapsed)
|
|
| 242 |
|
|
| 243 |
TimeIn=time.time() |
|
| 244 |
# Creation of vector for result with same size as input vectors |
|
| 245 |
res_np = np.empty_like(a_np) |
|
| 246 |
Elapsed=time.time()-TimeIn |
|
| 247 |
print("Allocation on Host for results: %.3f" % Elapsed)
|
|
| 248 |
|
|
| 249 |
TimeIn=time.time() |
|
| 250 |
# Copy from Device to Host |
|
| 251 |
cl.enqueue_copy(queue, res_np, res_g) |
|
| 252 |
Elapsed=time.time()-TimeIn |
|
| 253 |
print("Copy from Device 2 Host : %.3f" % Elapsed)
|
|
| 254 |
|
|
| 255 |
return(res_np) |
|
| 256 |
|
|
| 257 | 42 |
import sys |
| 258 | 43 |
import time |
| 259 | 44 |
|
| ... | ... | |
| 265 | 50 |
SIZE=int(sys.argv[1]) |
| 266 | 51 |
print("Size of vectors set to %i" % SIZE)
|
| 267 | 52 |
except: |
| 268 |
SIZE=50000
|
|
| 53 |
SIZE=256
|
|
| 269 | 54 |
print("Size of vectors set to default size %i" % SIZE)
|
| 270 | 55 |
|
| 271 |
# a_np = np.random.rand(SIZE).astype(np.float32) |
|
| 272 |
# b_np = np.random.rand(SIZE).astype(np.float32) |
|
| 273 |
|
|
| 274 | 56 |
a_np = np.ones(SIZE).astype(np.float32) |
| 275 | 57 |
b_np = np.ones(SIZE).astype(np.float32) |
| 276 | 58 |
|
| 59 |
C_np = np.zeros(SIZE).astype(np.float32) |
|
| 60 |
D_np = np.zeros(SIZE).astype(np.float32) |
|
| 61 |
C_np[0] = np.float32(SIZE) |
|
| 62 |
D_np[0] = np.float32(SIZE) |
|
| 63 |
|
|
| 277 | 64 |
# Native & Naive Implementation |
| 278 | 65 |
print("Performing naive implementation")
|
| 279 | 66 |
TimeIn=time.time() |
| ... | ... | |
| 281 | 68 |
NativeElapsed=time.time()-TimeIn |
| 282 | 69 |
NativeRate=int(SIZE/NativeElapsed) |
| 283 | 70 |
print("NativeRate: %i" % NativeRate)
|
| 284 |
|
|
| 71 |
print("Precision: ",np.linalg.norm(c_np-C_np),np.linalg.norm(d_np-D_np))
|
|
| 72 |
|
|
| 285 | 73 |
# Native & Numpy Implementation |
| 286 | 74 |
print("Performing Numpy implementation")
|
| 287 | 75 |
TimeIn=time.time() |
| ... | ... | |
| 289 | 77 |
NumpyElapsed=time.time()-TimeIn |
| 290 | 78 |
NumpyRate=int(SIZE/NumpyElapsed) |
| 291 | 79 |
print("NumpyRate: %i" % NumpyRate)
|
| 292 |
|
|
| 293 |
print(np.linalg.norm(c_np-e_np)) |
|
| 294 |
print(np.linalg.norm(d_np-f_np)) |
|
| 80 |
print("Precision: ",np.linalg.norm(e_np-C_np),np.linalg.norm(f_np-D_np))
|
|
| 295 | 81 |
|
| 296 |
# Native & Numpy Implementation
|
|
| 82 |
# Native & Numba Implementation
|
|
| 297 | 83 |
print("Performing Numba implementation")
|
| 298 | 84 |
TimeIn=time.time() |
| 299 | 85 |
g_np,h_np=NumbaDFT(a_np,b_np) |
| 300 |
NumpyElapsed=time.time()-TimeIn |
|
| 301 |
NumpyRate=int(SIZE/NumpyElapsed) |
|
| 302 |
print("NumpyRate: %i" % NumpyRate)
|
|
| 303 |
|
|
| 304 |
print(np.linalg.norm(c_np-g_np)) |
|
| 305 |
print(np.linalg.norm(d_np-h_np)) |
|
| 306 |
|
|
| 307 |
# # OpenCL Implementation |
|
| 308 |
# TimeIn=time.time() |
|
| 309 |
# # res_cl=OpenCLAddition(a_np,b_np) |
|
| 310 |
# res_cl=OpenCLSillyAddition(a_np,b_np) |
|
| 311 |
# OpenCLElapsed=time.time()-TimeIn |
|
| 312 |
# OpenCLRate=int(SIZE/OpenCLElapsed) |
|
| 313 |
# print("OpenCLRate: %i" % OpenCLRate)
|
|
| 314 |
|
|
| 315 |
# # CUDA Implementation |
|
| 316 |
# TimeIn=time.time() |
|
| 317 |
# # res_cuda=CUDAAddition(a_np,b_np) |
|
| 318 |
# res_cuda=CUDASillyAddition(a_np,b_np) |
|
| 319 |
# CUDAElapsed=time.time()-TimeIn |
|
| 320 |
# CUDARate=int(SIZE/CUDAElapsed) |
|
| 321 |
# print("CUDARate: %i" % CUDARate)
|
|
| 322 |
|
|
| 323 |
# print("OpenCLvsNative ratio: %f" % (OpenCLRate/NativeRate))
|
|
| 324 |
# print("CUDAvsNative ratio: %f" % (CUDARate/NativeRate))
|
|
| 325 |
|
|
| 326 |
# # Check on OpenCL with Numpy: |
|
| 327 |
# print(res_cl - res_np) |
|
| 328 |
# print(np.linalg.norm(res_cl - res_np)) |
|
| 329 |
# try: |
|
| 330 |
# assert np.allclose(res_np, res_cl) |
|
| 331 |
# except: |
|
| 332 |
# print("Results between Native & OpenCL seem to be too different!")
|
|
| 333 |
|
|
| 334 |
# # Check on CUDA with Numpy: |
|
| 335 |
# print(res_cuda - res_np) |
|
| 336 |
# print(np.linalg.norm(res_cuda - res_np)) |
|
| 337 |
# try: |
|
| 338 |
# assert np.allclose(res_np, res_cuda) |
|
| 339 |
# except: |
|
| 340 |
# print("Results between Native & CUDA seem to be too different!")
|
|
| 341 |
|
|
| 342 |
|
|
| 86 |
NumbaElapsed=time.time()-TimeIn |
|
| 87 |
NumbaRate=int(SIZE/NumbaElapsed) |
|
| 88 |
print("NumbaRate: %i" % NumbaRate)
|
|
| 89 |
print("Precision: ",np.linalg.norm(g_np-C_np),np.linalg.norm(h_np-D_np))
|
|
| ETSN/MyDFT_4.py (revision 271) | ||
|---|---|---|
| 1 |
#!/usr/bin/env python3 |
|
| 2 |
|
|
| 3 |
import numpy as np |
|
| 4 |
import pyopencl as cl |
|
| 5 |
from numpy import pi,cos,sin |
|
| 6 |
|
|
| 7 |
# Naive Discrete Fourier Transform |
|
| 8 |
def MyDFT(x,y): |
|
| 9 |
size=x.shape[0] |
|
| 10 |
X=np.zeros(size).astype(np.float32) |
|
| 11 |
Y=np.zeros(size).astype(np.float32) |
|
| 12 |
for i in range(size): |
|
| 13 |
for j in range(size): |
|
| 14 |
X[i]=X[i]+x[j]*cos(2.*pi*i*j/size)-y[j]*sin(2.*pi*i*j/size) |
|
| 15 |
Y[i]=Y[i]+x[j]*sin(2.*pi*i*j/size)+y[j]*cos(2.*pi*i*j/size) |
|
| 16 |
return(X,Y) |
|
| 17 |
|
|
| 18 |
# Numpy Discrete Fourier Transform |
|
| 19 |
def NumpyDFT(x,y): |
|
| 20 |
size=x.shape[0] |
|
| 21 |
X=np.zeros(size).astype(np.float32) |
|
| 22 |
Y=np.zeros(size).astype(np.float32) |
|
| 23 |
nj=np.multiply(2.0*np.pi/size,np.arange(size)).astype(np.float32) |
|
| 24 |
for i in range(size): |
|
| 25 |
X[i]=np.sum(np.subtract(np.multiply(np.cos(i*nj),x),np.multiply(np.sin(i*nj),y))) |
|
| 26 |
Y[i]=np.sum(np.add(np.multiply(np.sin(i*nj),x),np.multiply(np.cos(i*nj),y))) |
|
| 27 |
return(X,Y) |
|
| 28 |
|
|
| 29 |
# Numba Discrete Fourier Transform |
|
| 30 |
import numba |
|
| 31 |
@numba.njit(parallel=True) |
|
| 32 |
def NumbaDFT(x,y): |
|
| 33 |
size=x.shape[0] |
|
| 34 |
X=np.zeros(size).astype(np.float32) |
|
| 35 |
Y=np.zeros(size).astype(np.float32) |
|
| 36 |
nj=np.multiply(2.0*np.pi/size,np.arange(size)).astype(np.float32) |
|
| 37 |
for i in numba.prange(size): |
|
| 38 |
X[i]=np.sum(np.subtract(np.multiply(np.cos(i*nj),x),np.multiply(np.sin(i*nj),y))) |
|
| 39 |
Y[i]=np.sum(np.add(np.multiply(np.sin(i*nj),x),np.multiply(np.cos(i*nj),y))) |
|
| 40 |
return(X,Y) |
|
| 41 |
|
|
| 42 |
# OpenCL complete operation |
|
| 43 |
def OpenCLDFT(a_np,b_np): |
|
| 44 |
|
|
| 45 |
# Context creation |
|
| 46 |
ctx = cl.create_some_context() |
|
| 47 |
# Every process is stored in a queue |
|
| 48 |
queue = cl.CommandQueue(ctx) |
|
| 49 |
|
|
| 50 |
TimeIn=time.time() |
|
| 51 |
# Copy from Host to Device using pointers |
|
| 52 |
mf = cl.mem_flags |
|
| 53 |
a_g = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=a_np) |
|
| 54 |
b_g = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=b_np) |
|
| 55 |
Elapsed=time.time()-TimeIn |
|
| 56 |
print("Copy from Host 2 Device : %.3f" % Elapsed)
|
|
| 57 |
|
|
| 58 |
TimeIn=time.time() |
|
| 59 |
# Definition of kernel under OpenCL |
|
| 60 |
prg = cl.Program(ctx, """ |
|
| 61 |
|
|
| 62 |
#define PI 3.141592653589793 |
|
| 63 |
|
|
| 64 |
__kernel void MyDFT( |
|
| 65 |
__global const float *a_g, __global const float *b_g, __global float *A_g, __global float *B_g) |
|
| 66 |
{
|
|
| 67 |
int gid = get_global_id(0); |
|
| 68 |
uint size = get_global_size(0); |
|
| 69 |
float A=0.,B=0.; |
|
| 70 |
for (uint i=0; i<size;i++) |
|
| 71 |
{
|
|
| 72 |
A+=a_g[i]*cos(2.*PI*(float)(gid*i)/(float)size)-b_g[i]*sin(2.*PI*(float)(gid*i)/(float)size); |
|
| 73 |
B+=a_g[i]*sin(2.*PI*(float)(gid*i)/(float)size)+b_g[i]*cos(2.*PI*(float)(gid*i)/(float)size); |
|
| 74 |
} |
|
| 75 |
A_g[gid]=A; |
|
| 76 |
B_g[gid]=B; |
|
| 77 |
} |
|
| 78 |
""").build() |
|
| 79 |
Elapsed=time.time()-TimeIn |
|
| 80 |
print("Building kernels : %.3f" % Elapsed)
|
|
| 81 |
|
|
| 82 |
TimeIn=time.time() |
|
| 83 |
# Memory allocation on Device for result |
|
| 84 |
A_ocl = np.empty_like(a_np) |
|
| 85 |
B_ocl = np.empty_like(a_np) |
|
| 86 |
Elapsed=time.time()-TimeIn |
|
| 87 |
print("Allocation on Host for results : %.3f" % Elapsed)
|
|
| 88 |
|
|
| 89 |
A_g = cl.Buffer(ctx, mf.WRITE_ONLY, A_ocl.nbytes) |
|
| 90 |
B_g = cl.Buffer(ctx, mf.WRITE_ONLY, B_ocl.nbytes) |
|
| 91 |
Elapsed=time.time()-TimeIn |
|
| 92 |
print("Allocation on Device for results : %.3f" % Elapsed)
|
|
| 93 |
|
|
| 94 |
TimeIn=time.time() |
|
| 95 |
# Synthesis of function "sillysum" inside Kernel Sources |
|
| 96 |
knl = prg.MyDFT # Use this Kernel object for repeated calls |
|
| 97 |
Elapsed=time.time()-TimeIn |
|
| 98 |
print("Synthesis of kernel : %.3f" % Elapsed)
|
|
| 99 |
|
|
| 100 |
TimeIn=time.time() |
|
| 101 |
# Call of kernel previously defined |
|
| 102 |
CallCL=knl(queue, a_np.shape, None, a_g, b_g, A_g, B_g) |
|
| 103 |
# |
|
| 104 |
CallCL.wait() |
|
| 105 |
Elapsed=time.time()-TimeIn |
|
| 106 |
print("Execution of kernel : %.3f" % Elapsed)
|
|
| 107 |
|
|
| 108 |
TimeIn=time.time() |
|
| 109 |
# Copy from Device to Host |
|
| 110 |
cl.enqueue_copy(queue, A_ocl, A_g) |
|
| 111 |
cl.enqueue_copy(queue, B_ocl, B_g) |
|
| 112 |
Elapsed=time.time()-TimeIn |
|
| 113 |
print("Copy from Device 2 Host : %.3f" % Elapsed)
|
|
| 114 |
|
|
| 115 |
return(A_ocl,B_ocl) |
|
| 116 |
|
|
| 117 |
import sys |
|
| 118 |
import time |
|
| 119 |
|
|
| 120 |
if __name__=='__main__': |
|
| 121 |
|
|
| 122 |
# Size of input vectors definition based on stdin |
|
| 123 |
import sys |
|
| 124 |
try: |
|
| 125 |
SIZE=int(sys.argv[1]) |
|
| 126 |
print("Size of vectors set to %i" % SIZE)
|
|
| 127 |
except: |
|
| 128 |
SIZE=256 |
|
| 129 |
print("Size of vectors set to default size %i" % SIZE)
|
|
| 130 |
|
|
| 131 |
a_np = np.ones(SIZE).astype(np.float32) |
|
| 132 |
b_np = np.ones(SIZE).astype(np.float32) |
|
| 133 |
|
|
| 134 |
C_np = np.zeros(SIZE).astype(np.float32) |
|
| 135 |
D_np = np.zeros(SIZE).astype(np.float32) |
|
| 136 |
C_np[0] = np.float32(SIZE) |
|
| 137 |
D_np[0] = np.float32(SIZE) |
|
| 138 |
|
|
| 139 |
# # Native & Naive Implementation |
|
| 140 |
# print("Performing naive implementation")
|
|
| 141 |
# TimeIn=time.time() |
|
| 142 |
# c_np,d_np=MyDFT(a_np,b_np) |
|
| 143 |
# NativeElapsed=time.time()-TimeIn |
|
| 144 |
# NativeRate=int(SIZE/NativeElapsed) |
|
| 145 |
# print("NativeRate: %i" % NativeRate)
|
|
| 146 |
# print("Precision: ",np.linalg.norm(c_np-C_np),np.linalg.norm(d_np-D_np))
|
|
| 147 |
|
|
| 148 |
# Native & Numpy Implementation |
|
| 149 |
print("Performing Numpy implementation")
|
|
| 150 |
TimeIn=time.time() |
|
| 151 |
e_np,f_np=NumpyDFT(a_np,b_np) |
|
| 152 |
NumpyElapsed=time.time()-TimeIn |
|
| 153 |
NumpyRate=int(SIZE/NumpyElapsed) |
|
| 154 |
print("NumpyRate: %i" % NumpyRate)
|
|
| 155 |
print("Precision: ",np.linalg.norm(e_np-C_np),np.linalg.norm(f_np-D_np))
|
|
| 156 |
|
|
| 157 |
# Native & Numba Implementation |
|
| 158 |
print("Performing Numba implementation")
|
|
| 159 |
TimeIn=time.time() |
|
| 160 |
g_np,h_np=NumbaDFT(a_np,b_np) |
|
| 161 |
NumbaElapsed=time.time()-TimeIn |
|
| 162 |
NumbaRate=int(SIZE/NumbaElapsed) |
|
| 163 |
print("NumbaRate: %i" % NumbaRate)
|
|
| 164 |
print("Precision: ",np.linalg.norm(g_np-C_np),np.linalg.norm(h_np-D_np))
|
|
| 165 |
|
|
| 166 |
# OpenCL Implementation |
|
| 167 |
TimeIn=time.time() |
|
| 168 |
i_np,j_np=OpenCLDFT(a_np,b_np) |
|
| 169 |
OpenCLElapsed=time.time()-TimeIn |
|
| 170 |
OpenCLRate=int(SIZE/OpenCLElapsed) |
|
| 171 |
print("OpenCLRate: %i" % OpenCLRate)
|
|
| 172 |
print("Precision: ",np.linalg.norm(i_np-C_np),np.linalg.norm(j_np-D_np))
|
|
| 173 |
|
|
| 0 | 174 | |
| ETSN/MyDFT_5.py (revision 271) | ||
|---|---|---|
| 1 |
#!/usr/bin/env python3 |
|
| 2 |
|
|
| 3 |
import numpy as np |
|
| 4 |
import pyopencl as cl |
|
| 5 |
from numpy import pi,cos,sin |
|
| 6 |
|
|
| 7 |
# Naive Discrete Fourier Transform |
|
| 8 |
def MyDFT(x,y): |
|
| 9 |
size=x.shape[0] |
|
| 10 |
X=np.zeros(size).astype(np.float32) |
|
| 11 |
Y=np.zeros(size).astype(np.float32) |
|
| 12 |
for i in range(size): |
|
| 13 |
for j in range(size): |
|
| 14 |
X[i]=X[i]+x[j]*cos(2.*pi*i*j/size)-y[j]*sin(2.*pi*i*j/size) |
|
| 15 |
Y[i]=Y[i]+x[j]*sin(2.*pi*i*j/size)+y[j]*cos(2.*pi*i*j/size) |
|
| 16 |
return(X,Y) |
|
| 17 |
|
|
| 18 |
# Numpy Discrete Fourier Transform |
|
| 19 |
def NumpyDFT(x,y): |
|
| 20 |
size=x.shape[0] |
|
| 21 |
X=np.zeros(size).astype(np.float32) |
|
| 22 |
Y=np.zeros(size).astype(np.float32) |
|
| 23 |
nj=np.multiply(2.0*np.pi/size,np.arange(size)).astype(np.float32) |
|
| 24 |
for i in range(size): |
|
| 25 |
X[i]=np.sum(np.subtract(np.multiply(np.cos(i*nj),x),np.multiply(np.sin(i*nj),y))) |
|
| 26 |
Y[i]=np.sum(np.add(np.multiply(np.sin(i*nj),x),np.multiply(np.cos(i*nj),y))) |
|
| 27 |
return(X,Y) |
|
| 28 |
|
|
| 29 |
# Numba Discrete Fourier Transform |
|
| 30 |
import numba |
|
| 31 |
@numba.njit(parallel=True) |
|
| 32 |
def NumbaDFT(x,y): |
|
| 33 |
size=x.shape[0] |
|
| 34 |
X=np.zeros(size).astype(np.float32) |
|
| 35 |
Y=np.zeros(size).astype(np.float32) |
|
| 36 |
nj=np.multiply(2.0*np.pi/size,np.arange(size)).astype(np.float32) |
|
| 37 |
for i in numba.prange(size): |
|
| 38 |
X[i]=np.sum(np.subtract(np.multiply(np.cos(i*nj),x),np.multiply(np.sin(i*nj),y))) |
|
| 39 |
Y[i]=np.sum(np.add(np.multiply(np.sin(i*nj),x),np.multiply(np.cos(i*nj),y))) |
|
| 40 |
return(X,Y) |
|
| 41 |
|
|
| 42 |
# OpenCL complete operation |
|
| 43 |
def OpenCLDFT(a_np,b_np): |
|
| 44 |
|
|
| 45 |
# Context creation |
|
| 46 |
ctx = cl.create_some_context() |
|
| 47 |
# Every process is stored in a queue |
|
| 48 |
queue = cl.CommandQueue(ctx) |
|
| 49 |
|
|
| 50 |
TimeIn=time.time() |
|
| 51 |
# Copy from Host to Device using pointers |
|
| 52 |
mf = cl.mem_flags |
|
| 53 |
a_g = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=a_np) |
|
| 54 |
b_g = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=b_np) |
|
| 55 |
Elapsed=time.time()-TimeIn |
|
| 56 |
print("Copy from Host 2 Device : %.3f" % Elapsed)
|
|
| 57 |
|
|
| 58 |
TimeIn=time.time() |
|
| 59 |
# Definition of kernel under OpenCL |
|
| 60 |
prg = cl.Program(ctx, """ |
|
| 61 |
|
|
| 62 |
#define PI 3.141592653589793 |
|
| 63 |
|
|
| 64 |
__kernel void MyDFT( |
|
| 65 |
__global const float *a_g, __global const float *b_g, __global float *A_g, __global float *B_g) |
|
| 66 |
{
|
|
| 67 |
int gid = get_global_id(0); |
|
| 68 |
uint size = get_global_size(0); |
|
| 69 |
float A=0.,B=0.; |
|
| 70 |
for (uint i=0; i<size;i++) |
|
| 71 |
{
|
|
| 72 |
A+=a_g[i]*cos(2.*PI*(float)(gid*i)/(float)size)-b_g[i]*sin(2.*PI*(float)(gid*i)/(float)size); |
|
| 73 |
B+=a_g[i]*sin(2.*PI*(float)(gid*i)/(float)size)+b_g[i]*cos(2.*PI*(float)(gid*i)/(float)size); |
|
| 74 |
} |
|
| 75 |
A_g[gid]=A; |
|
| 76 |
B_g[gid]=B; |
|
| 77 |
} |
|
| 78 |
""").build() |
|
| 79 |
Elapsed=time.time()-TimeIn |
|
| 80 |
print("Building kernels : %.3f" % Elapsed)
|
|
| 81 |
|
|
| 82 |
TimeIn=time.time() |
|
| 83 |
# Memory allocation on Device for result |
|
| 84 |
A_ocl = np.empty_like(a_np) |
|
| 85 |
B_ocl = np.empty_like(a_np) |
|
| 86 |
Elapsed=time.time()-TimeIn |
|
| 87 |
print("Allocation on Host for results : %.3f" % Elapsed)
|
|
| 88 |
|
|
| 89 |
A_g = cl.Buffer(ctx, mf.WRITE_ONLY, A_ocl.nbytes) |
|
| 90 |
B_g = cl.Buffer(ctx, mf.WRITE_ONLY, B_ocl.nbytes) |
|
| 91 |
Elapsed=time.time()-TimeIn |
|
| 92 |
print("Allocation on Device for results : %.3f" % Elapsed)
|
|
| 93 |
|
|
| 94 |
TimeIn=time.time() |
|
| 95 |
# Synthesis of function "sillysum" inside Kernel Sources |
|
| 96 |
knl = prg.MyDFT # Use this Kernel object for repeated calls |
|
| 97 |
Elapsed=time.time()-TimeIn |
|
| 98 |
print("Synthesis of kernel : %.3f" % Elapsed)
|
|
| 99 |
|
|
| 100 |
TimeIn=time.time() |
|
| 101 |
# Call of kernel previously defined |
|
| 102 |
CallCL=knl(queue, a_np.shape, None, a_g, b_g, A_g, B_g) |
|
| 103 |
# |
|
| 104 |
CallCL.wait() |
|
| 105 |
Elapsed=time.time()-TimeIn |
|
| 106 |
print("Execution of kernel : %.3f" % Elapsed)
|
|
| 107 |
|
|
| 108 |
TimeIn=time.time() |
|
| 109 |
# Copy from Device to Host |
|
| 110 |
cl.enqueue_copy(queue, A_ocl, A_g) |
|
| 111 |
cl.enqueue_copy(queue, B_ocl, B_g) |
|
| 112 |
Elapsed=time.time()-TimeIn |
|
| 113 |
print("Copy from Device 2 Host : %.3f" % Elapsed)
|
|
| 114 |
|
|
| 115 |
return(A_ocl,B_ocl) |
|
| 116 |
|
|
| 117 |
# CUDA Silly complete operation |
|
| 118 |
def CUDADFT(a_np,b_np): |
|
| 119 |
import pycuda.autoinit |
|
| 120 |
import pycuda.driver as drv |
|
| 121 |
import numpy |
|
| 122 |
|
|
| 123 |
from pycuda.compiler import SourceModule |
|
| 124 |
TimeIn=time.time() |
|
| 125 |
mod = SourceModule("""
|
|
| 126 |
|
|
| 127 |
#define PI 3.141592653589793 |
|
| 128 |
|
|
| 129 |
__global__ void MyDFT(float *A_g, float *B_g, const float *a_g,const float *b_g) |
|
| 130 |
{
|
|
| 131 |
const int gid = blockIdx.x; |
|
| 132 |
uint size = gridDim.x; |
|
| 133 |
float A=0.,B=0.; |
|
| 134 |
for (uint i=0; i<size;i++) |
|
| 135 |
{
|
|
| 136 |
A+=a_g[i]*cos(2.*PI*(float)(gid*i)/(float)size)-b_g[i]*sin(2.*PI*(float)(gid*i)/(float)size); |
|
| 137 |
B+=a_g[i]*sin(2.*PI*(float)(gid*i)/(float)size)+b_g[i]*cos(2.*PI*(float)(gid*i)/(float)size); |
|
| 138 |
} |
|
| 139 |
A_g[gid]=A; |
|
| 140 |
B_g[gid]=B; |
|
| 141 |
} |
|
| 142 |
|
|
| 143 |
""") |
|
| 144 |
Elapsed=time.time()-TimeIn |
|
| 145 |
print("Definition of kernel : %.3f" % Elapsed)
|
|
| 146 |
|
|
| 147 |
TimeIn=time.time() |
|
| 148 |
MyDFT = mod.get_function("MyDFT")
|
|
| 149 |
Elapsed=time.time()-TimeIn |
|
| 150 |
print("Synthesis of kernel : %.3f" % Elapsed)
|
|
| 151 |
|
|
| 152 |
TimeIn=time.time() |
|
| 153 |
A_np = numpy.zeros_like(a_np) |
|
| 154 |
B_np = numpy.zeros_like(a_np) |
|
| 155 |
Elapsed=time.time()-TimeIn |
|
| 156 |
print("Allocation on Host for results : %.3f" % Elapsed)
|
|
| 157 |
|
|
| 158 |
TimeIn=time.time() |
|
| 159 |
MyDFT(drv.Out(A_np), drv.Out(B_np), drv.In(a_np), drv.In(b_np), |
|
| 160 |
block=(1,1,1), grid=(a_np.size,1)) |
|
| 161 |
Elapsed=time.time()-TimeIn |
|
| 162 |
print("Execution of kernel : %.3f" % Elapsed)
|
|
| 163 |
return(A_np,B_np) |
|
| 164 |
|
|
| 165 |
import sys |
|
| 166 |
import time |
|
| 167 |
|
|
| 168 |
if __name__=='__main__': |
|
| 169 |
|
|
| 170 |
# Size of input vectors definition based on stdin |
|
| 171 |
import sys |
|
| 172 |
try: |
|
| 173 |
SIZE=int(sys.argv[1]) |
|
| 174 |
print("Size of vectors set to %i" % SIZE)
|
|
| 175 |
except: |
|
| 176 |
SIZE=256 |
|
| 177 |
print("Size of vectors set to default size %i" % SIZE)
|
|
| 178 |
|
|
| 179 |
a_np = np.ones(SIZE).astype(np.float32) |
|
| 180 |
b_np = np.ones(SIZE).astype(np.float32) |
|
| 181 |
|
|
| 182 |
C_np = np.zeros(SIZE).astype(np.float32) |
|
| 183 |
D_np = np.zeros(SIZE).astype(np.float32) |
|
| 184 |
C_np[0] = np.float32(SIZE) |
|
| 185 |
D_np[0] = np.float32(SIZE) |
|
| 186 |
|
|
| 187 |
# Native & Naive Implementation |
|
| 188 |
print("Performing naive implementation")
|
|
| 189 |
TimeIn=time.time() |
|
| 190 |
c_np,d_np=MyDFT(a_np,b_np) |
|
| 191 |
NativeElapsed=time.time()-TimeIn |
|
| 192 |
NativeRate=int(SIZE/NativeElapsed) |
|
| 193 |
print("NativeRate: %i" % NativeRate)
|
|
| 194 |
print("Precision: ",np.linalg.norm(c_np-C_np),np.linalg.norm(d_np-D_np))
|
|
| 195 |
|
|
| 196 |
# Native & Numpy Implementation |
|
| 197 |
print("Performing Numpy implementation")
|
|
| 198 |
TimeIn=time.time() |
|
| 199 |
e_np,f_np=NumpyDFT(a_np,b_np) |
|
| 200 |
NumpyElapsed=time.time()-TimeIn |
|
| 201 |
NumpyRate=int(SIZE/NumpyElapsed) |
|
| 202 |
print("NumpyRate: %i" % NumpyRate)
|
|
| 203 |
print("Precision: ",np.linalg.norm(e_np-C_np),np.linalg.norm(f_np-D_np))
|
|
| 204 |
|
|
| 205 |
# Native & Numba Implementation |
|
| 206 |
print("Performing Numba implementation")
|
|
| 207 |
TimeIn=time.time() |
|
| 208 |
g_np,h_np=NumbaDFT(a_np,b_np) |
|
| 209 |
NumbaElapsed=time.time()-TimeIn |
|
| 210 |
NumbaRate=int(SIZE/NumbaElapsed) |
|
| 211 |
print("NumbaRate: %i" % NumbaRate)
|
|
| 212 |
print("Precision: ",np.linalg.norm(g_np-C_np),np.linalg.norm(h_np-D_np))
|
|
| 213 |
|
|
| 214 |
# OpenCL Implementation |
|
| 215 |
TimeIn=time.time() |
|
| 216 |
i_np,j_np=OpenCLDFT(a_np,b_np) |
|
| 217 |
OpenCLElapsed=time.time()-TimeIn |
|
| 218 |
OpenCLRate=int(SIZE/OpenCLElapsed) |
|
| 219 |
print("OpenCLRate: %i" % OpenCLRate)
|
|
| 220 |
print("Precision: ",np.linalg.norm(i_np-C_np),np.linalg.norm(j_np-D_np))
|
|
| 221 |
|
|
| 222 |
# CUDA Implementation |
|
| 223 |
TimeIn=time.time() |
|
| 224 |
k_np,l_np=CUDADFT(a_np,b_np) |
|
| 225 |
CUDAElapsed=time.time()-TimeIn |
|
| 226 |
CUDARate=int(SIZE/CUDAElapsed) |
|
| 227 |
print("CUDARate: %i" % CUDARate)
|
|
| 228 |
print("Precision: ",np.linalg.norm(k_np-C_np),np.linalg.norm(l_np-D_np))
|
|
| 229 |
|
|
| 0 | 230 | |
| ETSN/MyDFT_5b.py (revision 271) | ||
|---|---|---|
| 1 |
#!/usr/bin/env python3 |
|
| 2 |
|
|
| 3 |
import numpy as np |
|
| 4 |
import pyopencl as cl |
|
| 5 |
from numpy import pi,cos,sin |
|
| 6 |
|
|
| 7 |
# Naive Discrete Fourier Transform |
|
| 8 |
def MyDFT(x,y): |
|
| 9 |
size=x.shape[0] |
|
| 10 |
X=np.zeros(size).astype(np.float32) |
|
| 11 |
Y=np.zeros(size).astype(np.float32) |
|
| 12 |
for i in range(size): |
|
| 13 |
for j in range(size): |
|
| 14 |
X[i]=X[i]+x[j]*cos(2.*pi*i*j/size)-y[j]*sin(2.*pi*i*j/size) |
|
| 15 |
Y[i]=Y[i]+x[j]*sin(2.*pi*i*j/size)+y[j]*cos(2.*pi*i*j/size) |
|
| 16 |
return(X,Y) |
|
| 17 |
|
|
| 18 |
# Numpy Discrete Fourier Transform |
|
| 19 |
def NumpyDFT(x,y): |
|
| 20 |
size=x.shape[0] |
|
| 21 |
X=np.zeros(size).astype(np.float32) |
|
| 22 |
Y=np.zeros(size).astype(np.float32) |
|
| 23 |
nj=np.multiply(2.0*np.pi/size,np.arange(size)).astype(np.float32) |
|
| 24 |
for i in range(size): |
|
| 25 |
X[i]=np.sum(np.subtract(np.multiply(np.cos(i*nj),x),np.multiply(np.sin(i*nj),y))) |
|
| 26 |
Y[i]=np.sum(np.add(np.multiply(np.sin(i*nj),x),np.multiply(np.cos(i*nj),y))) |
|
| 27 |
return(X,Y) |
|
| 28 |
|
|
| 29 |
# Numba Discrete Fourier Transform |
|
| 30 |
import numba |
|
| 31 |
@numba.njit(parallel=True) |
|
| 32 |
def NumbaDFT(x,y): |
|
| 33 |
size=x.shape[0] |
|
| 34 |
X=np.zeros(size).astype(np.float32) |
|
| 35 |
Y=np.zeros(size).astype(np.float32) |
|
| 36 |
nj=np.multiply(2.0*np.pi/size,np.arange(size)).astype(np.float32) |
|
| 37 |
for i in numba.prange(size): |
|
| 38 |
X[i]=np.sum(np.subtract(np.multiply(np.cos(i*nj),x),np.multiply(np.sin(i*nj),y))) |
|
| 39 |
Y[i]=np.sum(np.add(np.multiply(np.sin(i*nj),x),np.multiply(np.cos(i*nj),y))) |
|
| 40 |
return(X,Y) |
|
| 41 |
|
|
| 42 |
# OpenCL complete operation |
|
| 43 |
def OpenCLDFT(a_np,b_np): |
|
| 44 |
|
|
| 45 |
# Context creation |
|
| 46 |
ctx = cl.create_some_context() |
|
| 47 |
# Every process is stored in a queue |
|
| 48 |
queue = cl.CommandQueue(ctx) |
|
| 49 |
|
|
| 50 |
TimeIn=time.time() |
|
| 51 |
# Copy from Host to Device using pointers |
|
| 52 |
mf = cl.mem_flags |
|
| 53 |
a_g = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=a_np) |
|
| 54 |
b_g = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=b_np) |
|
| 55 |
Elapsed=time.time()-TimeIn |
|
| 56 |
print("Copy from Host 2 Device : %.3f" % Elapsed)
|
|
| 57 |
|
|
| 58 |
TimeIn=time.time() |
|
| 59 |
# Definition of kernel under OpenCL |
|
| 60 |
prg = cl.Program(ctx, """ |
|
| 61 |
|
|
| 62 |
#define PI 3.141592653589793 |
|
| 63 |
|
|
| 64 |
__kernel void MyDFT( |
|
| 65 |
__global const float *a_g, __global const float *b_g, __global float *A_g, __global float *B_g) |
|
| 66 |
{
|
|
| 67 |
int gid = get_global_id(0); |
|
| 68 |
uint size = get_global_size(0); |
|
| 69 |
float A=0.,B=0.; |
|
| 70 |
for (uint i=0; i<size;i++) |
|
| 71 |
{
|
|
| 72 |
A+=a_g[i]*cos(2.*PI*(float)(gid*i)/(float)size)-b_g[i]*sin(2.*PI*(float)(gid*i)/(float)size); |
|
| 73 |
B+=a_g[i]*sin(2.*PI*(float)(gid*i)/(float)size)+b_g[i]*cos(2.*PI*(float)(gid*i)/(float)size); |
|
| 74 |
} |
|
| 75 |
A_g[gid]=A; |
|
| 76 |
B_g[gid]=B; |
|
| 77 |
} |
|
| 78 |
""").build() |
|
| 79 |
Elapsed=time.time()-TimeIn |
|
| 80 |
print("Building kernels : %.3f" % Elapsed)
|
|
| 81 |
|
|
| 82 |
TimeIn=time.time() |
|
| 83 |
# Memory allocation on Device for result |
|
| 84 |
A_ocl = np.empty_like(a_np) |
|
| 85 |
B_ocl = np.empty_like(a_np) |
|
| 86 |
Elapsed=time.time()-TimeIn |
|
| 87 |
print("Allocation on Host for results : %.3f" % Elapsed)
|
|
| 88 |
|
|
| 89 |
A_g = cl.Buffer(ctx, mf.WRITE_ONLY, A_ocl.nbytes) |
|
| 90 |
B_g = cl.Buffer(ctx, mf.WRITE_ONLY, B_ocl.nbytes) |
|
| 91 |
Elapsed=time.time()-TimeIn |
|
| 92 |
print("Allocation on Device for results : %.3f" % Elapsed)
|
|
| 93 |
|
|
| 94 |
TimeIn=time.time() |
|
| 95 |
# Synthesis of function "sillysum" inside Kernel Sources |
|
| 96 |
knl = prg.MyDFT # Use this Kernel object for repeated calls |
|
| 97 |
Elapsed=time.time()-TimeIn |
|
| 98 |
print("Synthesis of kernel : %.3f" % Elapsed)
|
|
| 99 |
|
|
| 100 |
TimeIn=time.time() |
|
| 101 |
# Call of kernel previously defined |
|
| 102 |
CallCL=knl(queue, a_np.shape, None, a_g, b_g, A_g, B_g) |
|
| 103 |
# |
|
| 104 |
CallCL.wait() |
|
| 105 |
Elapsed=time.time()-TimeIn |
|
| 106 |
print("Execution of kernel : %.3f" % Elapsed)
|
|
| 107 |
|
|
| 108 |
TimeIn=time.time() |
|
| 109 |
# Copy from Device to Host |
|
| 110 |
cl.enqueue_copy(queue, A_ocl, A_g) |
|
| 111 |
cl.enqueue_copy(queue, B_ocl, B_g) |
|
| 112 |
Elapsed=time.time()-TimeIn |
|
| 113 |
print("Copy from Device 2 Host : %.3f" % Elapsed)
|
|
| 114 |
|
|
| 115 |
return(A_ocl,B_ocl) |
|
| 116 |
|
|
| 117 |
# CUDA Silly complete operation |
|
| 118 |
def CUDADFT(a_np,b_np): |
|
| 119 |
import pycuda.autoinit |
|
| 120 |
import pycuda.driver as drv |
|
| 121 |
import numpy |
|
| 122 |
|
|
| 123 |
from pycuda.compiler import SourceModule |
|
| 124 |
TimeIn=time.time() |
|
| 125 |
mod = SourceModule("""
|
|
| 126 |
|
|
Formats disponibles : Unified diff