Révision 271 ETSN/MyDFT_1.py
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)) |
Formats disponibles : Unified diff