root/src/testbed/testbed.cpp @ 174:d204809217c6

Revision 174:d204809217c6, 13.3 KB (checked in by Andreas Schaefer <gentryx@…>, 14 months ago)

testing vectorized struct of arrays prototype
Line 
1#include <cmath>
2#include <typeinfo>
3// #include <CL/cl.h>
4#include <iostream>
5#include <emmintrin.h>
6// #include <pmmintrin.h>
7#include <sys/time.h>
8#include <vector>
9#include <libgeodecomp/misc/grid.h>
10
11using namespace LibGeoDecomp;
12
13class Scalar2D
14{
15public:
16    static int coefficients()
17    {
18        return 0;
19    }
20
21    inline void step(double *src, double *dst, int offset, int startX, int endX)
22    {
23        for (int x = startX; x < endX; ++x) {
24            dst[x] = (src[x - offset] + src[x - 1] + src[x] + src[x + 1] + src[x + offset]) * 0.2;
25        }
26    }
27
28    int flops()
29    {
30        return 5;
31    }
32};
33
34class VectorizedSSEMelbourneShuffle2D
35{
36public:
37    static int coefficients()
38    {
39        return 0;
40    }
41
42    inline void step(double *src, double *dst, int offset, int startX, int endX)
43    {
44        int x = startX;
45        Scalar2D scalarUpdater;
46
47        if ((x & 1) == 1) {
48            scalarUpdater.step(src, dst, offset, x, x + 1);
49            x += 1;
50        }
51
52        __m128d oneFifth = _mm_set_pd(1.0/3.0, 1.0/3.0);
53        __m128d buff0 = _mm_loadu_pd(src + x - 1);
54        __m128d same0 = _mm_load_pd(src + x + 0);
55
56        int paddedEndX = endX - 7;
57        for (; x < paddedEndX; x += 8) {
58            // load center row
59            __m128d same1 = _mm_load_pd(src + x + 2);
60            __m128d same2 = _mm_load_pd(src + x + 4);
61            __m128d same3 = _mm_load_pd(src + x + 6);
62            __m128d same4 = _mm_load_pd(src + x + 8);
63           
64            // shuffle values obtain left/right neighbors
65            __m128d buff1 = _mm_shuffle_pd(same0, same1, (1 << 0) | (0 << 2));
66            __m128d buff2 = _mm_shuffle_pd(same1, same2, (1 << 0) | (0 << 2));
67            __m128d buff3 = _mm_shuffle_pd(same2, same3, (1 << 0) | (0 << 2));
68            __m128d buff4 = _mm_shuffle_pd(same3, same4, (1 << 0) | (0 << 2));
69
70            // load top row
71            __m128d temp0 = _mm_load_pd(src - offset + x + 0);
72            __m128d temp1 = _mm_load_pd(src - offset + x + 2);
73            __m128d temp2 = _mm_load_pd(src - offset + x + 4);
74            __m128d temp3 = _mm_load_pd(src - offset + x + 6);
75
76            // add center row with left...
77            same0 = _mm_add_pd(same0, buff0);
78            same1 = _mm_add_pd(same1, buff1);
79            same2 = _mm_add_pd(same2, buff2);
80            same3 = _mm_add_pd(same3, buff3);
81
82            // ...and right neighbors
83            same0 = _mm_add_pd(same0, buff1);
84            same1 = _mm_add_pd(same1, buff2);
85            same2 = _mm_add_pd(same2, buff3);
86            same3 = _mm_add_pd(same3, buff4);
87   
88            // load bottom row
89            buff0 = _mm_load_pd(src + offset + x + 0);
90            buff1 = _mm_load_pd(src + offset + x + 2);
91            buff2 = _mm_load_pd(src + offset + x + 4);
92            buff3 = _mm_load_pd(src + offset + x + 6);
93       
94            // add top row
95            same0 = _mm_add_pd(same0, temp0);
96            same1 = _mm_add_pd(same1, temp1);
97            same2 = _mm_add_pd(same2, temp2);
98            same3 = _mm_add_pd(same3, temp3);
99
100            // add bottom row
101            same0 = _mm_add_pd(same0, buff0);
102            same1 = _mm_add_pd(same1, buff1);
103            same2 = _mm_add_pd(same2, buff2);
104            same3 = _mm_add_pd(same3, buff3);
105
106            // scale down...
107            same0 = _mm_mul_pd(same0, oneFifth);
108            same1 = _mm_mul_pd(same1, oneFifth);
109            same2 = _mm_mul_pd(same2, oneFifth);
110            same3 = _mm_mul_pd(same3, oneFifth);
111
112            // ...and store
113            _mm_store_pd(dst + 0, same0);
114            _mm_store_pd(dst + 2, same1);
115            _mm_store_pd(dst + 4, same2);
116            _mm_store_pd(dst + 6, same3);
117
118            same0 = same4;
119            buff0 = buff4;
120        }
121
122        scalarUpdater.step(src, dst, offset, x, endX);
123    }
124
125    int flops()
126    {
127        return 5;
128    }
129};
130
131
132template<typename UPDATER, int DIM>
133class Benchmark
134{
135public:
136    typedef Grid<double, typename Topologies::Cube<DIM>::Topology> GridType;
137
138
139    void run(Coord<DIM> dim, int repeats)
140    {
141        Coord<DIM> coeffDim = dim;
142        coeffDim.c[DIM - 1] *= UPDATER::coefficients();
143        GridType coeff(coeffDim, 0.1);
144        GridType a(dim, 1.0);
145        GridType b(dim, 1.0);
146
147        GridType *oldGrid = &a;
148        GridType *newGrid = &b;
149
150        UPDATER updater;
151
152        long long tStart = getUTtime();
153        std::vector<double*> coefficients(UPDATER::coefficients());
154
155        for (int t = 0; t < repeats; ++t) {
156            for (int z = 1; z < dim.c[2] - 1; ++z) {
157                for (int y = 1; y < dim.c[1] - 1; ++y) {
158                    Coord<DIM> c(0, y, z);
159                    Coord<DIM> coeffCoord = c;
160                    for (int i = 0; i < UPDATER::coefficients(); ++i) {
161                        coefficients[i] = &coeff.at(coeffCoord);
162                        coeffCoord.c[DIM - 1] += dim.c[DIM - 1];
163                    }
164                    updater.step(&coefficients[0], &oldGrid->at(c), &newGrid->at(c), dim.c[1], dim.c[2], 1, dim.c[0] - 1);
165                }
166            }
167            std::swap(newGrid, oldGrid);
168        }
169
170        long long tEnd = getUTtime();
171        evaluate(dim, repeats, tEnd - tStart);
172    }
173
174    void exercise() 
175    {
176        std::cout << "# " << typeid(UPDATER).name() << "\n";
177        int lastDim = 0;
178        for (int i = 4; i <= 4096; i *= 2) {
179            int intermediateSteps = 8;
180            for (int j = 0; j < intermediateSteps; ++j) {
181                int d = i * std::pow(2, j * (1.0 / intermediateSteps));
182                if (d % 2) {
183                    d += 1;
184                }
185
186                if (d > lastDim) {
187                    lastDim = d;
188                    Coord<DIM> dim;
189                    for (int i = 0; i < DIM; ++i)
190                        dim.c[i] = d;
191                    int repeats = std::max(1, 10000000 / dim.prod());
192                    run(dim, repeats);
193                }
194            }
195        }
196        std::cout << "\n";
197    }
198
199private:
200    long long getUTtime()
201    {
202        timeval t;
203        gettimeofday(&t, 0);
204        return (long long)t.tv_sec * 1000000 + t.tv_usec;
205    }
206
207    void evaluate(Coord<DIM> dim, int repeats, long long uTime)
208    {
209        double seconds = 1.0 * uTime / 1000 / 1000;
210        Coord<DIM> inner = dim;
211        for (int i = 0; i < DIM; ++i)
212            inner.c[i] -= 2;
213        double gflops = 1.0 * UPDATER().flops() * inner.prod() * 
214            repeats / 1000 / 1000 / 1000 / seconds;
215        std::cout << dim.x() << " " << gflops << "\n";
216    }
217
218
219};
220
221class Scalar3D
222{
223public:
224    static int coefficients()
225    {
226        return 7;
227    }
228
229    inline void step(double *coeff[7], double *src, double *dst, int offsetY, int offsetZ, int startX, int endX)
230    {
231        for (int x = startX; x < endX; ++x) {
232            dst[x] = 
233                coeff[0][x] * src[x - offsetZ] +
234                coeff[1][x] * src[x - offsetY] +
235                coeff[2][x] * src[x - 1] +
236                coeff[3][x] * src[x] +
237                coeff[4][x] * src[x + 1] +
238                coeff[5][x] * src[x + offsetY] +
239                coeff[6][x] * src[x + offsetZ];
240        }
241    }
242
243    int flops()
244    {
245        return 13;
246    }
247};
248
249class Vectorized3D
250{
251public:
252    static int coefficients()
253    {
254        return 7;
255    }
256
257    inline void step(double **coeff, double *src, double *dst, int offsetY, int offsetZ, int startX, int endX)
258    {
259        int x = startX;
260        Scalar3D scalarUpdater;
261
262        if ((x & 1) == 1) {
263            scalarUpdater.step(coeff, src, dst, offsetY, offsetZ, x, x + 1);
264            x += 1;
265        }
266
267        __m128d same0 = _mm_load_pd(src + x + 0);
268        __m128d neig0 = _mm_loadu_pd(src + x + 1);
269       
270        int paddedEndX = endX - 7;
271        for (; x < paddedEndX; x += 8) {
272            __m128d same1 = _mm_load_pd(src + x + 2);
273            __m128d same2 = _mm_load_pd(src + x + 4);
274            __m128d same3 = _mm_load_pd(src + x + 6);
275            __m128d same4 = _mm_load_pd(src + x + 8);
276
277            __m128d neig1 = _mm_shuffle_pd(same0, same1, (1 << 0) | (0 << 2));
278            __m128d neig2 = _mm_shuffle_pd(same1, same2, (1 << 0) | (0 << 2));
279            __m128d neig3 = _mm_shuffle_pd(same2, same3, (1 << 0) | (0 << 2));
280            __m128d neig4 = _mm_shuffle_pd(same3, same4, (1 << 0) | (0 << 2));
281
282            same0 = _mm_mul_pd(same0, _mm_load_pd(&coeff[3][x + 0]));
283            same1 = _mm_mul_pd(same1, _mm_load_pd(&coeff[3][x + 2]));
284            same2 = _mm_mul_pd(same2, _mm_load_pd(&coeff[3][x + 4]));
285            same3 = _mm_mul_pd(same3, _mm_load_pd(&coeff[3][x + 6]));
286
287            __m128d temp1 = _mm_mul_pd(neig0, _mm_load_pd(&coeff[2][x + 0]));
288            __m128d temp2 = _mm_mul_pd(neig1, _mm_load_pd(&coeff[2][x + 2]));
289            __m128d temp3 = _mm_mul_pd(neig2, _mm_load_pd(&coeff[2][x + 4]));
290            __m128d temp4 = _mm_mul_pd(neig3, _mm_load_pd(&coeff[2][x + 6]));
291
292            same0 = _mm_add_pd(same0, temp1);
293            same1 = _mm_add_pd(same1, temp2);
294            same2 = _mm_add_pd(same2, temp3);
295            same3 = _mm_add_pd(same3, temp4);
296
297            temp1 = _mm_mul_pd(neig1, _mm_load_pd(&coeff[4][x + 0]));
298            temp2 = _mm_mul_pd(neig2, _mm_load_pd(&coeff[4][x + 2]));
299            temp3 = _mm_mul_pd(neig3, _mm_load_pd(&coeff[4][x + 4]));
300            temp4 = _mm_mul_pd(neig4, _mm_load_pd(&coeff[4][x + 6]));
301
302            same0 = _mm_add_pd(same0, temp1);
303            same1 = _mm_add_pd(same1, temp2);
304            same2 = _mm_add_pd(same2, temp3);
305            same3 = _mm_add_pd(same3, temp4);
306
307            neig0 = _mm_load_pd(src + x - offsetZ + 0);
308            neig1 = _mm_load_pd(src + x - offsetZ + 2);
309            neig2 = _mm_load_pd(src + x - offsetZ + 4);
310            neig3 = _mm_load_pd(src + x - offsetZ + 6);
311
312            temp1 = _mm_mul_pd(neig0, _mm_load_pd(&coeff[0][x + 0]));
313            temp2 = _mm_mul_pd(neig1, _mm_load_pd(&coeff[0][x + 2]));
314            temp3 = _mm_mul_pd(neig2, _mm_load_pd(&coeff[0][x + 4]));
315            temp4 = _mm_mul_pd(neig3, _mm_load_pd(&coeff[0][x + 6]));
316
317            same0 = _mm_add_pd(same0, temp1);
318            same1 = _mm_add_pd(same1, temp2);
319            same2 = _mm_add_pd(same2, temp3);
320            same3 = _mm_add_pd(same3, temp4);
321
322            neig0 = _mm_load_pd(src + x - offsetY + 0);
323            neig1 = _mm_load_pd(src + x - offsetY + 2);
324            neig2 = _mm_load_pd(src + x - offsetY + 4);
325            neig3 = _mm_load_pd(src + x - offsetY + 6);
326
327            temp1 = _mm_mul_pd(neig0, _mm_load_pd(&coeff[1][x + 0]));
328            temp2 = _mm_mul_pd(neig1, _mm_load_pd(&coeff[1][x + 2]));
329            temp3 = _mm_mul_pd(neig2, _mm_load_pd(&coeff[1][x + 4]));
330            temp4 = _mm_mul_pd(neig3, _mm_load_pd(&coeff[1][x + 6]));
331
332            same0 = _mm_add_pd(same0, temp1);
333            same1 = _mm_add_pd(same1, temp2);
334            same2 = _mm_add_pd(same2, temp3);
335            same3 = _mm_add_pd(same3, temp4);
336
337            neig0 = _mm_load_pd(src + x + offsetY + 0);
338            neig1 = _mm_load_pd(src + x + offsetY + 2);
339            neig2 = _mm_load_pd(src + x + offsetY + 4);
340            neig3 = _mm_load_pd(src + x + offsetY + 6);
341
342            temp1 = _mm_mul_pd(neig0, _mm_load_pd(&coeff[5][x + 0]));
343            temp2 = _mm_mul_pd(neig1, _mm_load_pd(&coeff[5][x + 2]));
344            temp3 = _mm_mul_pd(neig2, _mm_load_pd(&coeff[5][x + 4]));
345            temp4 = _mm_mul_pd(neig3, _mm_load_pd(&coeff[5][x + 6]));
346
347            same0 = _mm_add_pd(same0, temp1);
348            same1 = _mm_add_pd(same1, temp2);
349            same2 = _mm_add_pd(same2, temp3);
350            same3 = _mm_add_pd(same3, temp4);
351
352            neig0 = _mm_load_pd(src + x + offsetZ + 0);
353            neig1 = _mm_load_pd(src + x + offsetZ + 2);
354            neig2 = _mm_load_pd(src + x + offsetZ + 4);
355            neig3 = _mm_load_pd(src + x + offsetZ + 6);
356
357            temp1 = _mm_mul_pd(neig0, _mm_load_pd(&coeff[5][x + 0]));
358            temp2 = _mm_mul_pd(neig1, _mm_load_pd(&coeff[5][x + 2]));
359            temp3 = _mm_mul_pd(neig2, _mm_load_pd(&coeff[5][x + 4]));
360            temp4 = _mm_mul_pd(neig3, _mm_load_pd(&coeff[5][x + 6]));
361
362            same0 = _mm_add_pd(same0, temp1);
363            same1 = _mm_add_pd(same1, temp2);
364            same2 = _mm_add_pd(same2, temp3);
365            same3 = _mm_add_pd(same3, temp4);
366
367            _mm_store_pd(dst + 0, same0);
368            _mm_store_pd(dst + 2, same1);
369            _mm_store_pd(dst + 4, same2);
370            _mm_store_pd(dst + 6, same3);
371
372            same0 = same4;
373            neig0 = neig4;
374
375            // dst[x] =
376            //     coeff[0][x] * src[x - offsetZ] +
377            //     coeff[1][x] * src[x - offsetY] +
378            //     coeff[2][x] * src[x - 1] +
379            //     coeff[3][x] * src[x] +
380            //     coeff[4][x] * src[x + 1] +
381            //     coeff[5][x] * src[x + offsetY] +
382            //     coeff[6][x] * src[x + offsetZ];
383        }
384
385        scalarUpdater.step(coeff, src, dst, offsetY, offsetZ, x, endX);
386    }
387
388    int flops()
389    {
390        return 13;
391    }
392};
393
394
395int main(int argc, char *argv[])
396{
397    // Benchmark<Scalar3D, 3>().exercise();
398    Benchmark<Vectorized3D, 3>().exercise();
399    // Benchmark<VectorizedSSEMelbourneShuffle2D>().exercise();
400
401
402    // std::vector<cl::Platform> platforms;
403    // cl::Platform::get(&platforms);
404
405    // for (int i = 0; i < platforms.size(); ++i) {
406    //     std::string str;
407    //     platforms[i].getInfo(CL_PLATFORM_NAME, &str);
408    //     std::cout << "Platform[" << i << "] = " << str << std::endl;
409    // }
410
411    return 0;
412}
Note: See TracBrowser for help on using the browser.