Révision 109
Epidevomath/vector.py (revision 109) | ||
---|---|---|
1 |
import numpy as np |
|
2 |
import pyopencl as cl |
|
3 |
import pyopencl.array as cl_array |
|
4 |
from numpy.random import randint as nprnd |
|
5 |
|
|
6 |
deviceID = 0 |
|
7 |
platformID = 0 |
|
8 |
workGroup=(1,1) |
|
9 |
|
|
10 |
N = 10 |
|
11 |
MyData = np.zeros(N, dtype=cl_array.vec.float4) |
|
12 |
|
|
13 |
dev = cl.get_platforms()[platformID].get_devices()[deviceID] |
|
14 |
|
|
15 |
ctx = cl.Context([dev]) |
|
16 |
queue = cl.CommandQueue(ctx) |
|
17 |
mf = cl.mem_flags |
|
18 |
clData = cl.Buffer(ctx, mf.READ_WRITE, MyData.nbytes) |
|
19 |
|
|
20 |
|
|
21 |
prg = cl.Program(ctx, """ |
|
22 |
#define znew ((z=36969*(z&65535)+(z>>16))<<16) |
|
23 |
#define wnew ((w=18000*(w&65535)+(w>>16))&65535) |
|
24 |
#define MWC (znew+wnew) |
|
25 |
#define SHR3 (jsr=(jsr=(jsr=jsr^(jsr<<17))^(jsr>>13))^(jsr<<5)) |
|
26 |
#define CONG (jcong=69069*jcong+1234567) |
|
27 |
#define KISS ((MWC^CONG)+SHR3) |
|
28 |
|
|
29 |
#define MWCfp MWC * 2.328306435454494e-10f |
|
30 |
#define KISSfp KISS * 2.328306435454494e-10f |
|
31 |
#define SHR3fp SHR3 * 2.328306435454494e-10f |
|
32 |
#define CONGfp CONG * 2.328306435454494e-10f |
|
33 |
|
|
34 |
__kernel void SplutterSpace(__global float4* clData, |
|
35 |
uint seed_z,uint seed_w) |
|
36 |
{ |
|
37 |
int gid = get_global_id(0); |
|
38 |
uint z=seed_z+(uint)gid; |
|
39 |
uint w=seed_w-(uint)gid; |
|
40 |
|
|
41 |
clData[gid].xyzw = (float4) (MWCfp,MWCfp,MWCfp,0.); |
|
42 |
} |
|
43 |
|
|
44 |
|
|
45 |
""").build() |
|
46 |
|
|
47 |
#prg.Pack_Cmplx(queue, (N,1), workGroup, Data_In, np.int32(N)) |
|
48 |
|
|
49 |
prg.SplutterSpace(queue, (N,1), None, clData, |
|
50 |
numpy.uint32(nprnd(2**32)),numpy.uint32(nprnd(2**32))) |
|
51 |
cl.enqueue_copy(queue, MyData, clData) |
|
52 |
|
|
53 |
|
|
54 |
print MyData |
Epidevomath/float.py (revision 109) | ||
---|---|---|
1 |
import numpy as np |
|
2 |
import pyopencl as cl |
|
3 |
import pyopencl.array as cl_array |
|
4 |
|
|
5 |
|
|
6 |
deviceID = 0 |
|
7 |
platformID = 0 |
|
8 |
workGroup=(1,1) |
|
9 |
|
|
10 |
N = 10 |
|
11 |
testData = np.zeros(N, dtype=cl_array.vec.float4) |
|
12 |
|
|
13 |
dev = cl.get_platforms()[platformID].get_devices()[deviceID] |
|
14 |
|
|
15 |
ctx = cl.Context([dev]) |
|
16 |
queue = cl.CommandQueue(ctx) |
|
17 |
mf = cl.mem_flags |
|
18 |
Data_In = cl.Buffer(ctx, mf.READ_WRITE, testData.nbytes) |
|
19 |
|
|
20 |
|
|
21 |
prg = cl.Program(ctx, """ |
|
22 |
|
|
23 |
__kernel void Pack_Cmplx( __global float4* Data_In, int N) |
|
24 |
{ |
|
25 |
int gid = get_global_id(0); |
|
26 |
|
|
27 |
Data_In[gid] = 1; |
|
28 |
} |
|
29 |
""").build() |
|
30 |
|
|
31 |
prg.Pack_Cmplx(queue, (N,1), workGroup, Data_In, np.int32(N)) |
|
32 |
cl.enqueue_copy(queue, testData, Data_In) |
|
33 |
|
|
34 |
|
|
35 |
print testData |
Epidevomath/vector4.py (revision 109) | ||
---|---|---|
1 |
import numpy as np |
|
2 |
import pyopencl as cl |
|
3 |
import pyopencl.array as cl_array |
|
4 |
from numpy.random import randint as nprnd |
|
5 |
|
|
6 |
deviceID = 0 |
|
7 |
platformID = 0 |
|
8 |
workGroup=(1,1) |
|
9 |
|
|
10 |
N = 10 |
|
11 |
MyData = np.zeros(N, dtype=cl_array.vec.float4) |
|
12 |
|
|
13 |
dev = cl.get_platforms()[platformID].get_devices()[deviceID] |
|
14 |
|
|
15 |
ctx = cl.Context([dev]) |
|
16 |
queue = cl.CommandQueue(ctx) |
|
17 |
mf = cl.mem_flags |
|
18 |
clData = cl.Buffer(ctx, mf.READ_WRITE, MyData.nbytes) |
|
19 |
|
|
20 |
|
|
21 |
prg = cl.Program(ctx, """ |
|
22 |
#define znew ((z=36969*(z&65535)+(z>>16))<<16) |
|
23 |
#define wnew ((w=18000*(w&65535)+(w>>16))&65535) |
|
24 |
#define MWC (znew+wnew) |
|
25 |
#define SHR3 (jsr=(jsr=(jsr=jsr^(jsr<<17))^(jsr>>13))^(jsr<<5)) |
|
26 |
#define CONG (jcong=69069*jcong+1234567) |
|
27 |
#define KISS ((MWC^CONG)+SHR3) |
|
28 |
|
|
29 |
#define MWCfp MWC * 2.328306435454494e-10f |
|
30 |
#define KISSfp KISS * 2.328306435454494e-10f |
|
31 |
#define SHR3fp SHR3 * 2.328306435454494e-10f |
|
32 |
#define CONGfp CONG * 2.328306435454494e-10f |
|
33 |
|
|
34 |
__kernel void SplutterSpace(__global float4* clData, |
|
35 |
uint seed_z,uint seed_w) |
|
36 |
{ |
|
37 |
int gid = get_global_id(0); |
|
38 |
uint z=seed_z+(uint)gid; |
|
39 |
uint w=seed_w-(uint)gid; |
|
40 |
|
|
41 |
clData[gid].xyzw = (float4) (MWCfp,MWCfp,MWCfp,0.); |
|
42 |
} |
|
43 |
|
|
44 |
|
|
45 |
""").build() |
|
46 |
|
|
47 |
#prg.Pack_Cmplx(queue, (N,1), workGroup, Data_In, np.int32(N)) |
|
48 |
|
|
49 |
prg.SplutterSpace(queue, (N,1), None, clData, |
|
50 |
numpy.uint32(nprnd(2**32)),numpy.uint32(nprnd(2**32))) |
|
51 |
cl.enqueue_copy(queue, MyData, clData) |
|
52 |
|
|
53 |
|
|
54 |
print MyData |
Epidevomath/vector8.py (revision 109) | ||
---|---|---|
1 |
import numpy as np |
|
2 |
import pyopencl as cl |
|
3 |
import pyopencl.array as cl_array |
|
4 |
from numpy.random import randint as nprnd |
|
5 |
|
|
6 |
deviceID = 0 |
|
7 |
platformID = 0 |
|
8 |
workGroup=(1,1) |
|
9 |
|
|
10 |
N = 32768 |
|
11 |
MyData = np.zeros(N, dtype=cl_array.vec.float8) |
|
12 |
|
|
13 |
dev = cl.get_platforms()[platformID].get_devices()[deviceID] |
|
14 |
|
|
15 |
ctx = cl.Context([dev]) |
|
16 |
queue = cl.CommandQueue(ctx) |
|
17 |
mf = cl.mem_flags |
|
18 |
clData = cl.Buffer(ctx, mf.READ_WRITE, MyData.nbytes) |
|
19 |
|
|
20 |
MyRoutines = cl.Program(ctx, """ |
|
21 |
#define znew ((z=36969*(z&65535)+(z>>16))<<16) |
|
22 |
#define wnew ((w=18000*(w&65535)+(w>>16))&65535) |
|
23 |
#define MWC (znew+wnew) |
|
24 |
#define SHR3 (jsr=(jsr=(jsr=jsr^(jsr<<17))^(jsr>>13))^(jsr<<5)) |
|
25 |
#define CONG (jcong=69069*jcong+1234567) |
|
26 |
#define KISS ((MWC^CONG)+SHR3) |
|
27 |
|
|
28 |
#define MWCfp MWC * 2.328306435454494e-10f |
|
29 |
#define KISSfp KISS * 2.328306435454494e-10f |
|
30 |
#define SHR3fp SHR3 * 2.328306435454494e-10f |
|
31 |
#define CONGfp CONG * 2.328306435454494e-10f |
|
32 |
|
|
33 |
#define LENGTH 1. |
|
34 |
|
|
35 |
#define PI 3.14159265359 |
|
36 |
|
|
37 |
#define SMALL_NUM 0.0000001 |
|
38 |
|
|
39 |
__kernel void SplutterSpace(__global float8* clData, |
|
40 |
uint seed_z,uint seed_w) |
|
41 |
{ |
|
42 |
int gid = get_global_id(0); |
|
43 |
uint z=seed_z+(uint)gid; |
|
44 |
uint w=seed_w-(uint)gid; |
|
45 |
|
|
46 |
clData[gid].s01234567 = (float8) (MWCfp,MWCfp,MWCfp,0.,0.,0.,0.,0.); |
|
47 |
} |
|
48 |
|
|
49 |
__kernel void ExtendSegment(__global float8* clData, |
|
50 |
uint seed_z,uint seed_w) |
|
51 |
{ |
|
52 |
int gid = get_global_id(0); |
|
53 |
uint z=seed_z+(uint)gid; |
|
54 |
uint w=seed_w-(uint)gid; |
|
55 |
|
|
56 |
float theta=MWCfp*PI; |
|
57 |
float phi=MWCfp*PI*2.; |
|
58 |
float sinTheta=sin(theta); |
|
59 |
clData[gid].s4=clData[gid].s0+LENGTH*sinTheta*cos(phi); |
|
60 |
clData[gid].s5=clData[gid].s1+LENGTH*sinTheta*sin(phi); |
|
61 |
clData[gid].s6=clData[gid].s2+LENGTH*cos(theta); |
|
62 |
|
|
63 |
} |
|
64 |
|
|
65 |
__kernel void EstimateLength(__global float8* clData,__global float* clSize) |
|
66 |
{ |
|
67 |
int gid = get_global_id(0); |
|
68 |
|
|
69 |
clSize[gid]=distance(clData[gid].lo,clData[gid].hi); |
|
70 |
} |
|
71 |
|
|
72 |
// Get from http://geomalgorithms.com/a07-_distance.html |
|
73 |
|
|
74 |
__kernel void ShortestDistance(__global float8* clData,__global float* clDistance) |
|
75 |
{ |
|
76 |
int gidx = get_global_id(0); |
|
77 |
int ggsz = get_global_size(0); |
|
78 |
int gidy = get_global_id(1); |
|
79 |
|
|
80 |
float4 u = clData[gidx].hi - clData[gidx].lo; |
|
81 |
float4 v = clData[gidy].hi - clData[gidy].lo; |
|
82 |
float4 w = clData[gidx].lo - clData[gidy].lo; |
|
83 |
|
|
84 |
float a = dot(u,u); // always >= 0 |
|
85 |
float b = dot(u,v); |
|
86 |
float c = dot(v,v); // always >= 0 |
|
87 |
float d = dot(u,w); |
|
88 |
float e = dot(v,w); |
|
89 |
|
|
90 |
float D = a*c - b*b; // always >= 0 |
|
91 |
float sc, sN, sD = D; // sc = sN / sD, default sD = D >= 0 |
|
92 |
float tc, tN, tD = D; // tc = tN / tD, default tD = D >= 0 |
|
93 |
|
|
94 |
// compute the line parameters of the two closest points |
|
95 |
if (D < SMALL_NUM) { // the lines are almost parallel |
|
96 |
sN = 0.0; // force using point P0 on segment S1 |
|
97 |
sD = 1.0; // to prevent possible division by 0.0 later |
|
98 |
tN = e; |
|
99 |
tD = c; |
|
100 |
} |
|
101 |
else { // get the closest points on the infinite lines |
|
102 |
sN = (b*e - c*d); |
|
103 |
tN = (a*e - b*d); |
|
104 |
if (sN < 0.0) { // sc < 0 => the s=0 edge is visible |
|
105 |
sN = 0.0; |
|
106 |
tN = e; |
|
107 |
tD = c; |
|
108 |
} |
|
109 |
else if (sN > sD) { // sc > 1 => the s=1 edge is visible |
|
110 |
sN = sD; |
|
111 |
tN = e + b; |
|
112 |
tD = c; |
|
113 |
} |
|
114 |
} |
|
115 |
|
|
116 |
if (tN < 0.0) { // tc < 0 => the t=0 edge is visible |
|
117 |
tN = 0.0; |
|
118 |
// recompute sc for this edge |
|
119 |
if (-d < 0.0) |
|
120 |
sN = 0.0; |
|
121 |
else if (-d > a) |
|
122 |
sN = sD; |
|
123 |
else { |
|
124 |
sN = -d; |
|
125 |
sD = a; |
|
126 |
} |
|
127 |
} |
|
128 |
else if (tN > tD) { // tc > 1 => the t=1 edge is visible |
|
129 |
tN = tD; |
|
130 |
// recompute sc for this edge |
|
131 |
if ((-d + b) < 0.0) |
|
132 |
sN = 0; |
|
133 |
else if ((-d + b) > a) |
|
134 |
sN = sD; |
|
135 |
else { |
|
136 |
sN = (-d + b); |
|
137 |
sD = a; |
|
138 |
} |
|
139 |
} |
|
140 |
// finally do the division to get sc and tc |
|
141 |
sc = (fabs(sN) < SMALL_NUM ? 0.0 : sN / sD); |
|
142 |
tc = (fabs(tN) < SMALL_NUM ? 0.0 : tN / tD); |
|
143 |
|
|
144 |
// get the difference of the two closest points |
|
145 |
float4 dP = w + (sc * u) - (tc * v); // = S1(sc) - S2(tc) |
|
146 |
|
|
147 |
clDistance[ggsz*gidy+gidx]=length(dP); // return the closest distance |
|
148 |
} |
|
149 |
|
|
150 |
""").build() |
|
151 |
|
|
152 |
print 'Tous au meme endroit',MyData |
|
153 |
|
|
154 |
MyRoutines.SplutterSpace(queue, (N,1), None, clData, numpy.uint32(nprnd(2**32)),numpy.uint32(nprnd(2**32))) |
|
155 |
|
|
156 |
cl.enqueue_copy(queue, MyData, clData) |
|
157 |
|
|
158 |
print 'Tous distribues',MyData |
|
159 |
|
|
160 |
MyRoutines.ExtendSegment(queue, (N,1), None, clData,numpy.uint32(nprnd(2**32)),numpy.uint32(nprnd(2**32))) |
|
161 |
|
|
162 |
cl.enqueue_copy(queue, MyData, clData) |
|
163 |
|
|
164 |
print 'Tous avec leur extremite',MyData |
|
165 |
|
|
166 |
MySize = np.zeros(len(MyData), dtype=numpy.float32) |
|
167 |
clSize = cl.Buffer(ctx, mf.READ_WRITE, MySize.nbytes) |
|
168 |
|
|
169 |
MyRoutines.EstimateLength(queue, (N,1), None, clData, clSize) |
|
170 |
cl.enqueue_copy(queue, MySize, clSize) |
|
171 |
|
|
172 |
print 'La distance de chacun avec son extremite',MySize |
|
173 |
|
|
174 |
MyDistance = np.zeros(len(MyData)*len(MyData), dtype=numpy.float32) |
|
175 |
clDistance = cl.Buffer(ctx, mf.READ_WRITE, MyDistance.nbytes) |
|
176 |
|
|
177 |
MyRoutines.ShortestDistance(queue, (N,N), None, clData, clDistance) |
|
178 |
cl.enqueue_copy(queue, MyDistance, clDistance) |
|
179 |
|
|
180 |
MyDistance=numpy.reshape(MyDistance,(N,N)) |
|
181 |
|
|
182 |
print 'La distance de chacun',MyDistance |
|
183 |
|
Formats disponibles : Unified diff