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