root/src/testbed/testbed.cpp @ 183:c9f7e84ad779

Revision 183:c9f7e84ad779, 102.2 KB (checked in by Andreas Schaefer <gentryx@…>, 13 months ago)

slow evolutions towards a fast, convenient API
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// FC = FixedCoord
132template<int DIM_X=0, int DIM_Y=0, int DIM_Z=0>
133class FC
134{
135public:
136};
137
138template<int DIM_X, int DIM_Y, int DIM_Z, int INDEX=0>
139class Neighborhood
140{
141public:
142    Neighborhood(double *_coefficients, double *_source) :
143        coefficients(_coefficients),
144        source(_source)
145    {}
146
147    template<int X, int Y, int Z, int C> 
148    const double& coeff(const int& x) const
149    {
150        return coefficients[C * DIM_X * DIM_Y * DIM_Z + Z * DIM_X * DIM_Y + Y * DIM_X + X + x];
151    }
152
153    template<int X, int Y, int Z> 
154    const double& src(const int& x) const
155    {
156        return source[Z * DIM_X * DIM_Y + Y * DIM_X + X + x];
157    }
158
159    template<int X, int Y, int Z>
160    Neighborhood<DIM_X, DIM_Y, DIM_Z, X + Y * DIM_X + Z * DIM_X * DIM_Y>
161    operator[](FC<X, Y, Z> /*unused*/) const
162    {
163        return Neighborhood<DIM_X, DIM_Y, DIM_Z, X + Y * DIM_X + Z * DIM_X * DIM_Y>(coefficients, source);
164    }
165
166    const double& srcB(const int& x) const
167    {
168        return source[INDEX + x];
169    }
170
171    template<int C> 
172    const double& coeffB(const FC<C, 0, 0> /*unused*/, const int& x) const
173    { 
174        return coefficients[C * DIM_X * DIM_Y * DIM_Z + INDEX];
175    }
176
177private:
178    double *coefficients;
179    double *source;
180};
181
182template<typename UPDATER, int DIM>
183class Benchmark
184{
185public:
186    typedef Grid<double, typename Topologies::Cube<DIM>::Topology> GridType;
187
188
189    void run(Coord<DIM> dim, int repeats)
190    {
191        Coord<DIM> coeffDim = dim;
192        coeffDim.c[DIM - 1] *= UPDATER::coefficients();
193        GridType coeff(coeffDim, 0.1);
194        GridType a(dim, 1.0);
195        GridType b(dim, 1.0);
196
197        GridType *oldGrid = &a;
198        GridType *newGrid = &b;
199
200        UPDATER updater;
201
202        long long tStart = getUTtime();
203        std::vector<double*> coefficients(UPDATER::coefficients());
204
205        for (int t = 0; t < repeats; ++t) {
206            for (int z = 1; z < dim.c[2] - 1; ++z) {
207                for (int y = 1; y < dim.c[1] - 1; ++y) {
208                    Coord<DIM> c(0, y, z);
209                    Coord<DIM> coeffCoord = c;
210
211                    // Defiant
212                    updater.step(
213                        Neighborhood<4, 4, MY_SIZE>(&coeff.at(coeffCoord), &oldGrid->at(c)),
214                        &newGrid->at(c), 
215                        1, 
216                        dim.c[0] - 1);
217
218                    // updater.step(
219                    //     Neighborhood<4, 4, MY_SIZE>(&coeff.at(coeffCoord), &oldGrid->at(c)),
220                    //     &newGrid->at(c),
221                    //     dim.c[0],
222                    //     dim.c[0] * dim.c[1],
223                    //     1,
224                    //     dim.c[0] - 1);
225                   
226                    // for (int i = 0; i < UPDATER::coefficients(); ++i) {
227                    //     coefficients[i] = &coeff.at(coeffCoord);
228                    //     coeffCoord.c[DIM - 1] += dim.c[DIM - 1];
229                    // }
230                    // updater.step(
231                    //     &coefficients[0],
232                    //     &oldGrid->at(c),
233                    //     &newGrid->at(c),
234                    //     dim.c[0],
235                    //     dim.c[0] * dim.c[1],
236                    //     1,
237                    //     dim.c[0] - 1);
238
239                    // double *source[9] = {
240                    //     &oldGrid->at(Coord<DIM>(0, y - 1, z + 1)),
241                    //     &oldGrid->at(Coord<DIM>(0, y - 1, z    )),
242                    //     &oldGrid->at(Coord<DIM>(0, y - 1, z - 1)),
243                    //     &oldGrid->at(Coord<DIM>(0, y,     z + 1)),
244                    //     &oldGrid->at(Coord<DIM>(0, y,     z    )),
245                    //     &oldGrid->at(Coord<DIM>(0, y,     z - 1)),
246                    //     &oldGrid->at(Coord<DIM>(0, y + 1, z + 1)),
247                    //     &oldGrid->at(Coord<DIM>(0, y + 1, z    )),
248                    //     &oldGrid->at(Coord<DIM>(0, y + 1, z - 1))
249                    // };
250                    // updater.update(source, &newGrid->at(c), 1, dim.c[0] - 1);
251                }
252            }
253            std::swap(newGrid, oldGrid);
254        }
255
256        long long tEnd = getUTtime();
257        evaluate(dim, repeats, tEnd - tStart);
258    }
259
260    void exercise() 
261    {
262        std::cout << "# " << typeid(UPDATER).name() << "\n";
263        int lastDim = 0;
264        for (int i = 4; i <= 4096; i *= 2) {
265            int intermediateSteps = 8;
266            for (int j = 0; j < intermediateSteps; ++j) {
267                int d = i * std::pow(2, j * (1.0 / intermediateSteps));
268                if (d % 2) {
269                    d += 1;
270                }
271
272                if (d > lastDim) {
273                    lastDim = d;
274                    Coord<DIM> dim;
275                    dim.c[0] = d;
276                    dim.c[1] = 4;
277                    dim.c[2] = 4;
278                    // for (int i = 0; i < DIM; ++i)
279                    //     dim.c[i] = d;
280                    int repeats = std::max(1, 10000000 / dim.prod());
281                    run(dim, repeats);
282                }
283            }
284        }
285        std::cout << "\n";
286    }
287
288private:
289    long long getUTtime()
290    {
291        timeval t;
292        gettimeofday(&t, 0);
293        return (long long)t.tv_sec * 1000000 + t.tv_usec;
294    }
295
296    void evaluate(Coord<DIM> dim, int repeats, long long uTime)
297    {
298        double seconds = 1.0 * uTime / 1000 / 1000;
299        Coord<DIM> inner = dim;
300        for (int i = 0; i < DIM; ++i)
301            inner.c[i] -= 2;
302        double gflops = 1.0 * UPDATER().flops() * inner.prod() * 
303            repeats / 1000 / 1000 / 1000 / seconds;
304        std::cout << dim.x() << " " << gflops << "\n";
305    }
306
307
308};
309
310#define TN 0
311#define T  1
312#define TS 2
313#define N  3
314#define C  4
315#define S  5
316#define BN 6
317#define B  7
318#define BS 8
319
320class Jacobi3D
321{
322public:
323    static int coefficients()
324    {
325        return 1;
326    }
327
328    inline void update(double **src, double *dst, int startX, int endX)
329    {
330       int x = startX;
331
332       if ((x & 1) == 1) {
333           updateScalar(src, dst, x, x + 1);
334           x += 1;
335       }
336
337       __m128d oneSeventh = _mm_set_pd(1.0/7.0, 1.0/7.0);
338       __m128d buff0 = _mm_loadu_pd(src[C] + x - 1);
339       __m128d same0 = _mm_load_pd(src[C] + x + 0);
340
341       int paddedEndX = endX - 7;
342       for (; x < paddedEndX; x += 8) {
343           // load center row
344           __m128d same1 = _mm_load_pd(src[0] + x + 2);
345           __m128d same2 = _mm_load_pd(src[0] + x + 4);
346           __m128d same3 = _mm_load_pd(src[0] + x + 6);
347           __m128d same4 = _mm_load_pd(src[0] + x + 8);
348           
349           // shuffle values obtain left/right neighbors
350           __m128d buff1 = _mm_shuffle_pd(same0, same1, (1 << 0) | (0 << 2));
351           __m128d buff2 = _mm_shuffle_pd(same1, same2, (1 << 0) | (0 << 2));
352           __m128d buff3 = _mm_shuffle_pd(same2, same3, (1 << 0) | (0 << 2));
353           __m128d buff4 = _mm_shuffle_pd(same3, same4, (1 << 0) | (0 << 2));
354   
355           // load top row
356           __m128d temp0 = _mm_load_pd(src[T] + x + 0);
357           __m128d temp1 = _mm_load_pd(src[T] + x + 2);
358           __m128d temp2 = _mm_load_pd(src[T] + x + 4);
359           __m128d temp3 = _mm_load_pd(src[T] + x + 6);
360
361           // add center row with left...
362           same0 = _mm_add_pd(same0, buff0);
363           same1 = _mm_add_pd(same1, buff1);
364           same2 = _mm_add_pd(same2, buff2);
365           same3 = _mm_add_pd(same3, buff3);
366
367           // ...and right neighbors
368           same0 = _mm_add_pd(same0, buff1);
369           same1 = _mm_add_pd(same1, buff2);
370           same2 = _mm_add_pd(same2, buff3);
371           same3 = _mm_add_pd(same3, buff4);
372   
373           // load bottom row
374           buff0 = _mm_load_pd(src[B] + x + 0);
375           buff1 = _mm_load_pd(src[B] + x + 2);
376           buff2 = _mm_load_pd(src[B] + x + 4);
377           buff3 = _mm_load_pd(src[B] + x + 6);
378       
379           // add top row
380           same0 = _mm_add_pd(same0, temp0);
381           same1 = _mm_add_pd(same1, temp1);
382           same2 = _mm_add_pd(same2, temp2);
383           same3 = _mm_add_pd(same3, temp3);
384
385           // load north row
386           temp0 = _mm_load_pd(src[N] + x + 0);
387           temp1 = _mm_load_pd(src[N] + x + 2);
388           temp2 = _mm_load_pd(src[N] + x + 4);
389           temp3 = _mm_load_pd(src[N] + x + 6);
390
391           // add bottom row
392           same0 = _mm_add_pd(same0, buff0);
393           same1 = _mm_add_pd(same1, buff1);
394           same2 = _mm_add_pd(same2, buff2);
395           same3 = _mm_add_pd(same3, buff3);
396           
397           // load south row
398           buff0 = _mm_load_pd(src[S] + x + 0);
399           buff1 = _mm_load_pd(src[S] + x + 2);
400           buff2 = _mm_load_pd(src[S] + x + 4);
401           buff3 = _mm_load_pd(src[S] + x + 6);
402
403           // add north row
404           same0 = _mm_add_pd(same0, temp0);
405           same1 = _mm_add_pd(same1, temp1);
406           same2 = _mm_add_pd(same2, temp2);
407           same3 = _mm_add_pd(same3, temp3);
408
409           // add south row
410           same0 = _mm_add_pd(same0, buff0);
411           same1 = _mm_add_pd(same1, buff1);
412           same2 = _mm_add_pd(same2, buff2);
413           same3 = _mm_add_pd(same3, buff3);
414
415           // scale down...
416           same0 = _mm_mul_pd(same0, oneSeventh);
417           same1 = _mm_mul_pd(same1, oneSeventh);
418           same2 = _mm_mul_pd(same2, oneSeventh);
419           same3 = _mm_mul_pd(same3, oneSeventh);
420
421           // ...and store
422           _mm_store_pd(dst + x + 0, same0);
423           _mm_store_pd(dst + x + 2, same1);
424           _mm_store_pd(dst + x + 4, same2);
425           _mm_store_pd(dst + x + 6, same3);
426
427           same0 = same4;
428           buff0 = buff4;
429       }
430
431       updateScalar(src, dst, x, endX);
432    }
433
434    inline void updateScalar(double *src[9], double *dst, int startX, int endX)
435    {
436        for (int x = startX; x < endX; ++x) {
437            dst[x] = 
438                (src[S][x] +
439                 src[T][x] +
440                 src[C][x - 1] +
441                 src[C][x] +
442                 src[C][x + 1] +
443                 src[B][x] +
444                 src[N][x]) * (1.0 / 7.0);
445        }
446    }
447
448    int flops()
449    {
450        return 8;
451    }
452};
453
454
455class Scalar3D
456{
457public:
458    static int coefficients()
459    {
460        return 7;
461    }
462
463    inline void step(double *coeff[7], double *src, double *dst, int offsetY, int offsetZ, int startX, int endX)
464    {
465        for (int x = startX; x < endX; ++x) {
466            dst[x] = 
467                coeff[0][x] * src[x - offsetZ] +
468                coeff[1][x] * src[x - offsetY] +
469                coeff[2][x] * src[x - 1] +
470                coeff[3][x] * src[x] +
471                coeff[4][x] * src[x + 1] +
472                coeff[5][x] * src[x + offsetY] +
473                coeff[6][x] * src[x + offsetZ];
474        }
475    }
476
477    int flops()
478    {
479        return 13;
480    }
481};
482
483class Vectorized3D
484{
485public:
486    static int coefficients()
487    {
488        return 7;
489    }
490
491    inline void step(double **coeff, double *src, double *dst, int offsetY, int offsetZ, int startX, int endX)
492    {
493        int x = startX;
494        Scalar3D scalarUpdater;
495
496        if ((x & 1) == 1) {
497            scalarUpdater.step(coeff, src, dst, offsetY, offsetZ, x, x + 1);
498            x += 1;
499        }
500
501        __m128d same0 = _mm_load_pd(src + x + 0);
502        __m128d neig0 = _mm_loadu_pd(src + x + 1);
503       
504        int paddedEndX = endX - 7;
505        for (; x < paddedEndX; x += 8) {
506            __m128d same1 = _mm_load_pd(src + x + 2);
507            __m128d same2 = _mm_load_pd(src + x + 4);
508            __m128d same3 = _mm_load_pd(src + x + 6);
509            __m128d same4 = _mm_load_pd(src + x + 8);
510
511            __m128d neig1 = _mm_shuffle_pd(same0, same1, (1 << 0) | (0 << 2));
512            __m128d neig2 = _mm_shuffle_pd(same1, same2, (1 << 0) | (0 << 2));
513            __m128d neig3 = _mm_shuffle_pd(same2, same3, (1 << 0) | (0 << 2));
514            __m128d neig4 = _mm_shuffle_pd(same3, same4, (1 << 0) | (0 << 2));
515
516            same0 = _mm_mul_pd(same0, _mm_load_pd(&coeff[3][x + 0]));
517            same1 = _mm_mul_pd(same1, _mm_load_pd(&coeff[3][x + 2]));
518            same2 = _mm_mul_pd(same2, _mm_load_pd(&coeff[3][x + 4]));
519            same3 = _mm_mul_pd(same3, _mm_load_pd(&coeff[3][x + 6]));
520
521            __m128d temp1 = _mm_mul_pd(neig0, _mm_load_pd(&coeff[2][x + 0]));
522            __m128d temp2 = _mm_mul_pd(neig1, _mm_load_pd(&coeff[2][x + 2]));
523            __m128d temp3 = _mm_mul_pd(neig2, _mm_load_pd(&coeff[2][x + 4]));
524            __m128d temp4 = _mm_mul_pd(neig3, _mm_load_pd(&coeff[2][x + 6]));
525
526            same0 = _mm_add_pd(same0, temp1);
527            same1 = _mm_add_pd(same1, temp2);
528            same2 = _mm_add_pd(same2, temp3);
529            same3 = _mm_add_pd(same3, temp4);
530
531            temp1 = _mm_mul_pd(neig1, _mm_load_pd(&coeff[4][x + 0]));
532            temp2 = _mm_mul_pd(neig2, _mm_load_pd(&coeff[4][x + 2]));
533            temp3 = _mm_mul_pd(neig3, _mm_load_pd(&coeff[4][x + 4]));
534            temp4 = _mm_mul_pd(neig4, _mm_load_pd(&coeff[4][x + 6]));
535
536            same0 = _mm_add_pd(same0, temp1);
537            same1 = _mm_add_pd(same1, temp2);
538            same2 = _mm_add_pd(same2, temp3);
539            same3 = _mm_add_pd(same3, temp4);
540
541            neig0 = _mm_load_pd(src + x - offsetZ + 0);
542            neig1 = _mm_load_pd(src + x - offsetZ + 2);
543            neig2 = _mm_load_pd(src + x - offsetZ + 4);
544            neig3 = _mm_load_pd(src + x - offsetZ + 6);
545
546            temp1 = _mm_mul_pd(neig0, _mm_load_pd(&coeff[0][x + 0]));
547            temp2 = _mm_mul_pd(neig1, _mm_load_pd(&coeff[0][x + 2]));
548            temp3 = _mm_mul_pd(neig2, _mm_load_pd(&coeff[0][x + 4]));
549            temp4 = _mm_mul_pd(neig3, _mm_load_pd(&coeff[0][x + 6]));
550
551            same0 = _mm_add_pd(same0, temp1);
552            same1 = _mm_add_pd(same1, temp2);
553            same2 = _mm_add_pd(same2, temp3);
554            same3 = _mm_add_pd(same3, temp4);
555
556            neig0 = _mm_load_pd(src + x - offsetY + 0);
557            neig1 = _mm_load_pd(src + x - offsetY + 2);
558            neig2 = _mm_load_pd(src + x - offsetY + 4);
559            neig3 = _mm_load_pd(src + x - offsetY + 6);
560
561            temp1 = _mm_mul_pd(neig0, _mm_load_pd(&coeff[1][x + 0]));
562            temp2 = _mm_mul_pd(neig1, _mm_load_pd(&coeff[1][x + 2]));
563            temp3 = _mm_mul_pd(neig2, _mm_load_pd(&coeff[1][x + 4]));
564            temp4 = _mm_mul_pd(neig3, _mm_load_pd(&coeff[1][x + 6]));
565
566            same0 = _mm_add_pd(same0, temp1);
567            same1 = _mm_add_pd(same1, temp2);
568            same2 = _mm_add_pd(same2, temp3);
569            same3 = _mm_add_pd(same3, temp4);
570
571            neig0 = _mm_load_pd(src + x + offsetY + 0);
572            neig1 = _mm_load_pd(src + x + offsetY + 2);
573            neig2 = _mm_load_pd(src + x + offsetY + 4);
574            neig3 = _mm_load_pd(src + x + offsetY + 6);
575
576            temp1 = _mm_mul_pd(neig0, _mm_load_pd(&coeff[5][x + 0]));
577            temp2 = _mm_mul_pd(neig1, _mm_load_pd(&coeff[5][x + 2]));
578            temp3 = _mm_mul_pd(neig2, _mm_load_pd(&coeff[5][x + 4]));
579            temp4 = _mm_mul_pd(neig3, _mm_load_pd(&coeff[5][x + 6]));
580
581            same0 = _mm_add_pd(same0, temp1);
582            same1 = _mm_add_pd(same1, temp2);
583            same2 = _mm_add_pd(same2, temp3);
584            same3 = _mm_add_pd(same3, temp4);
585
586            neig0 = _mm_load_pd(src + x + offsetZ + 0);
587            neig1 = _mm_load_pd(src + x + offsetZ + 2);
588            neig2 = _mm_load_pd(src + x + offsetZ + 4);
589            neig3 = _mm_load_pd(src + x + offsetZ + 6);
590
591            temp1 = _mm_mul_pd(neig0, _mm_load_pd(&coeff[6][x + 0]));
592            temp2 = _mm_mul_pd(neig1, _mm_load_pd(&coeff[6][x + 2]));
593            temp3 = _mm_mul_pd(neig2, _mm_load_pd(&coeff[6][x + 4]));
594            temp4 = _mm_mul_pd(neig3, _mm_load_pd(&coeff[6][x + 6]));
595
596            same0 = _mm_add_pd(same0, temp1);
597            same1 = _mm_add_pd(same1, temp2);
598            same2 = _mm_add_pd(same2, temp3);
599            same3 = _mm_add_pd(same3, temp4);
600
601            _mm_store_pd(dst + 0, same0);
602            _mm_store_pd(dst + 2, same1);
603            _mm_store_pd(dst + 4, same2);
604            _mm_store_pd(dst + 6, same3);
605
606            same0 = same4;
607            neig0 = neig4;
608
609            // dst[x] =
610            //     coeff[0][x] * src[x - offsetZ] +
611            //     coeff[1][x] * src[x - offsetY] +
612            //     coeff[2][x] * src[x - 1] +
613            //     coeff[3][x] * src[x] +
614            //     coeff[4][x] * src[x + 1] +
615            //     coeff[5][x] * src[x + offsetY] +
616            //     coeff[6][x] * src[x + offsetZ];
617        }
618
619        scalarUpdater.step(coeff, src, dst, offsetY, offsetZ, x, endX);
620    }
621
622    int flops()
623    {
624        return 13;
625    }
626};
627
628class ExtendedScalar3D
629{
630public:
631    static int coefficients()
632    {
633        return 13;
634    }
635
636    inline void step(double *coeff[13], double *src, double *dst, int offsetY, int offsetZ, int startX, int endX)
637    {
638        for (int x = startX; x < endX; ++x) {
639            dst[x] = 
640                coeff[ 0][x] * src[x - offsetZ - offsetY] +
641                coeff[ 1][x] * src[x - offsetZ] +
642                coeff[ 2][x] * src[x - offsetZ + offsetY] +
643                coeff[ 3][x] * src[x - offsetY] +
644                coeff[ 4][x] * src[x - 1] +
645                coeff[ 5][x] * src[x] +
646                coeff[ 6][x] * src[x + 1] +
647                coeff[ 7][x] * src[x + offsetY] +
648                coeff[ 8][x] * src[x + offsetZ - offsetY] +
649                coeff[ 9][x] * src[x + offsetZ] +
650                coeff[10][x] * src[x + offsetZ + offsetY] +
651                coeff[11][x - offsetZ - offsetY] +
652                coeff[11][x - offsetZ] +
653                coeff[11][x - offsetZ + offsetY] +
654                coeff[11][x - offsetY] +
655                coeff[11][x] +
656                coeff[11][x + offsetY] +
657                coeff[11][x + offsetZ - offsetY] +
658                coeff[11][x + offsetZ] +
659                coeff[11][x + offsetZ + offsetY] +
660                coeff[12][x - offsetZ - offsetY] +
661                coeff[12][x - offsetZ] +
662                coeff[12][x - offsetZ + offsetY] +
663                coeff[12][x - offsetY] +
664                coeff[12][x] +
665                coeff[12][x + offsetY] +
666                coeff[12][x + offsetZ - offsetY] +
667                coeff[12][x + offsetZ] +
668                coeff[12][x + offsetZ + offsetY];
669        }
670    }
671
672    int flops()
673    {
674        return 40;
675    }
676};
677
678class ExtendedVectorized3D
679{
680public:
681    static int coefficients()
682    {
683        return 13;
684    }
685
686    inline void step(double *coeff[13], double *src, double *dst, int offsetY, int offsetZ, int startX, int endX)
687    {
688        int x = startX;
689        ExtendedScalar3D scalarUpdater;
690
691        if ((x & 1) == 1) {
692            scalarUpdater.step(coeff, src, dst, offsetY, offsetZ, x, x + 1);
693            x += 1;
694        }
695
696        __m128d same0 = _mm_load_pd(src + x + 0);
697        __m128d neig0 = _mm_loadu_pd(src + x + 1);
698       
699        int paddedEndX = endX - 7;
700        for (; x < paddedEndX; x += 8) {
701            __m128d same1 = _mm_load_pd(src + x + 2);
702            __m128d same2 = _mm_load_pd(src + x + 4);
703            __m128d same3 = _mm_load_pd(src + x + 6);
704            __m128d same4 = _mm_load_pd(src + x + 8);
705
706            __m128d neig1 = _mm_shuffle_pd(same0, same1, (1 << 0) | (0 << 2));
707            __m128d neig2 = _mm_shuffle_pd(same1, same2, (1 << 0) | (0 << 2));
708            __m128d neig3 = _mm_shuffle_pd(same2, same3, (1 << 0) | (0 << 2));
709            __m128d neig4 = _mm_shuffle_pd(same3, same4, (1 << 0) | (0 << 2));
710
711            same0 = _mm_mul_pd(same0, _mm_load_pd(&coeff[3][x + 0]));
712            same1 = _mm_mul_pd(same1, _mm_load_pd(&coeff[3][x + 2]));
713            same2 = _mm_mul_pd(same2, _mm_load_pd(&coeff[3][x + 4]));
714            same3 = _mm_mul_pd(same3, _mm_load_pd(&coeff[3][x + 6]));
715
716            __m128d temp1 = _mm_mul_pd(neig0, _mm_load_pd(&coeff[2][x + 0]));
717            __m128d temp2 = _mm_mul_pd(neig1, _mm_load_pd(&coeff[2][x + 2]));
718            __m128d temp3 = _mm_mul_pd(neig2, _mm_load_pd(&coeff[2][x + 4]));
719            __m128d temp4 = _mm_mul_pd(neig3, _mm_load_pd(&coeff[2][x + 6]));
720
721            same0 = _mm_add_pd(same0, temp1);
722            same1 = _mm_add_pd(same1, temp2);
723            same2 = _mm_add_pd(same2, temp3);
724            same3 = _mm_add_pd(same3, temp4);
725
726            temp1 = _mm_mul_pd(neig1, _mm_load_pd(&coeff[4][x + 0]));
727            temp2 = _mm_mul_pd(neig2, _mm_load_pd(&coeff[4][x + 2]));
728            temp3 = _mm_mul_pd(neig3, _mm_load_pd(&coeff[4][x + 4]));
729            temp4 = _mm_mul_pd(neig4, _mm_load_pd(&coeff[4][x + 6]));
730
731            same0 = _mm_add_pd(same0, temp1);
732            same1 = _mm_add_pd(same1, temp2);
733            same2 = _mm_add_pd(same2, temp3);
734            same3 = _mm_add_pd(same3, temp4);
735
736            neig0 = _mm_load_pd(src + x - offsetZ + 0);
737            neig1 = _mm_load_pd(src + x - offsetZ + 2);
738            neig2 = _mm_load_pd(src + x - offsetZ + 4);
739            neig3 = _mm_load_pd(src + x - offsetZ + 6);
740
741            temp1 = _mm_mul_pd(neig0, _mm_load_pd(&coeff[0][x + 0]));
742            temp2 = _mm_mul_pd(neig1, _mm_load_pd(&coeff[0][x + 2]));
743            temp3 = _mm_mul_pd(neig2, _mm_load_pd(&coeff[0][x + 4]));
744            temp4 = _mm_mul_pd(neig3, _mm_load_pd(&coeff[0][x + 6]));
745
746            same0 = _mm_add_pd(same0, temp1);
747            same1 = _mm_add_pd(same1, temp2);
748            same2 = _mm_add_pd(same2, temp3);
749            same3 = _mm_add_pd(same3, temp4);
750
751            neig0 = _mm_load_pd(src + x - offsetY + 0);
752            neig1 = _mm_load_pd(src + x - offsetY + 2);
753            neig2 = _mm_load_pd(src + x - offsetY + 4);
754            neig3 = _mm_load_pd(src + x - offsetY + 6);
755
756            temp1 = _mm_mul_pd(neig0, _mm_load_pd(&coeff[1][x + 0]));
757            temp2 = _mm_mul_pd(neig1, _mm_load_pd(&coeff[1][x + 2]));
758            temp3 = _mm_mul_pd(neig2, _mm_load_pd(&coeff[1][x + 4]));
759            temp4 = _mm_mul_pd(neig3, _mm_load_pd(&coeff[1][x + 6]));
760
761            same0 = _mm_add_pd(same0, temp1);
762            same1 = _mm_add_pd(same1, temp2);
763            same2 = _mm_add_pd(same2, temp3);
764            same3 = _mm_add_pd(same3, temp4);
765
766            neig0 = _mm_load_pd(src + x + offsetY + 0);
767            neig1 = _mm_load_pd(src + x + offsetY + 2);
768            neig2 = _mm_load_pd(src + x + offsetY + 4);
769            neig3 = _mm_load_pd(src + x + offsetY + 6);
770
771            temp1 = _mm_mul_pd(neig0, _mm_load_pd(&coeff[5][x + 0]));
772            temp2 = _mm_mul_pd(neig1, _mm_load_pd(&coeff[5][x + 2]));
773            temp3 = _mm_mul_pd(neig2, _mm_load_pd(&coeff[5][x + 4]));
774            temp4 = _mm_mul_pd(neig3, _mm_load_pd(&coeff[5][x + 6]));
775
776            same0 = _mm_add_pd(same0, temp1);
777            same1 = _mm_add_pd(same1, temp2);
778            same2 = _mm_add_pd(same2, temp3);
779            same3 = _mm_add_pd(same3, temp4);
780
781            //xxxxxxxxxxxxx
782            neig0 = _mm_load_pd(src + x + offsetZ + 0);
783            neig1 = _mm_load_pd(src + x + offsetZ + 2);
784            neig2 = _mm_load_pd(src + x + offsetZ + 4);
785            neig3 = _mm_load_pd(src + x + offsetZ + 6);
786
787            temp1 = _mm_mul_pd(neig0, _mm_load_pd(&coeff[6][x + 0]));
788            temp2 = _mm_mul_pd(neig1, _mm_load_pd(&coeff[6][x + 2]));
789            temp3 = _mm_mul_pd(neig2, _mm_load_pd(&coeff[6][x + 4]));
790            temp4 = _mm_mul_pd(neig3, _mm_load_pd(&coeff[6][x + 6]));
791
792            same0 = _mm_add_pd(same0, temp1);
793            same1 = _mm_add_pd(same1, temp2);
794            same2 = _mm_add_pd(same2, temp3);
795            same3 = _mm_add_pd(same3, temp4);
796
797            //xxxxxxxxxxxxx
798            neig0 = _mm_load_pd(src + x - offsetZ - offsetY + 0);
799            neig1 = _mm_load_pd(src + x - offsetZ - offsetY + 2);
800            neig2 = _mm_load_pd(src + x - offsetZ - offsetY + 4);
801            neig3 = _mm_load_pd(src + x - offsetZ - offsetY + 6);
802
803            temp1 = _mm_mul_pd(neig0, _mm_load_pd(&coeff[7][x + 0]));
804            temp2 = _mm_mul_pd(neig1, _mm_load_pd(&coeff[7][x + 2]));
805            temp3 = _mm_mul_pd(neig2, _mm_load_pd(&coeff[7][x + 4]));
806            temp4 = _mm_mul_pd(neig3, _mm_load_pd(&coeff[7][x + 6]));
807
808            same0 = _mm_add_pd(same0, temp1);
809            same1 = _mm_add_pd(same1, temp2);
810            same2 = _mm_add_pd(same2, temp3);
811            same3 = _mm_add_pd(same3, temp4);
812
813            //xxxxxxxxxxxxx
814            neig0 = _mm_load_pd(src + x - offsetZ + offsetY + 0);
815            neig1 = _mm_load_pd(src + x - offsetZ + offsetY + 2);
816            neig2 = _mm_load_pd(src + x - offsetZ + offsetY + 4);
817            neig3 = _mm_load_pd(src + x - offsetZ + offsetY + 6);
818
819            temp1 = _mm_mul_pd(neig0, _mm_load_pd(&coeff[8][x + 0]));
820            temp2 = _mm_mul_pd(neig1, _mm_load_pd(&coeff[8][x + 2]));
821            temp3 = _mm_mul_pd(neig2, _mm_load_pd(&coeff[8][x + 4]));
822            temp4 = _mm_mul_pd(neig3, _mm_load_pd(&coeff[8][x + 6]));
823
824            same0 = _mm_add_pd(same0, temp1);
825            same1 = _mm_add_pd(same1, temp2);
826            same2 = _mm_add_pd(same2, temp3);
827            same3 = _mm_add_pd(same3, temp4);
828
829            //xxxxxxxxxxxxx
830            neig0 = _mm_load_pd(src + x + offsetZ - offsetY + 0);
831            neig1 = _mm_load_pd(src + x + offsetZ - offsetY + 2);
832            neig2 = _mm_load_pd(src + x + offsetZ - offsetY + 4);
833            neig3 = _mm_load_pd(src + x + offsetZ - offsetY + 6);
834
835            temp1 = _mm_mul_pd(neig0, _mm_load_pd(&coeff[9][x + 0]));
836            temp2 = _mm_mul_pd(neig1, _mm_load_pd(&coeff[9][x + 2]));
837            temp3 = _mm_mul_pd(neig2, _mm_load_pd(&coeff[9][x + 4]));
838            temp4 = _mm_mul_pd(neig3, _mm_load_pd(&coeff[9][x + 6]));
839
840            same0 = _mm_add_pd(same0, temp1);
841            same1 = _mm_add_pd(same1, temp2);
842            same2 = _mm_add_pd(same2, temp3);
843            same3 = _mm_add_pd(same3, temp4);
844
845            //xxxxxxxxxxxxx
846            neig0 = _mm_load_pd(src + x + offsetZ + offsetY + 0);
847            neig1 = _mm_load_pd(src + x + offsetZ + offsetY + 2);
848            neig2 = _mm_load_pd(src + x + offsetZ + offsetY + 4);
849            neig3 = _mm_load_pd(src + x + offsetZ + offsetY + 6);
850
851            temp1 = _mm_mul_pd(neig0, _mm_load_pd(&coeff[10][x + 0]));
852            temp2 = _mm_mul_pd(neig1, _mm_load_pd(&coeff[10][x + 2]));
853            temp3 = _mm_mul_pd(neig2, _mm_load_pd(&coeff[10][x + 4]));
854            temp4 = _mm_mul_pd(neig3, _mm_load_pd(&coeff[10][x + 6]));
855
856            same0 = _mm_add_pd(same0, temp1);
857            same1 = _mm_add_pd(same1, temp2);
858            same2 = _mm_add_pd(same2, temp3);
859            same3 = _mm_add_pd(same3, temp4);
860
861            //xxxxxxxxxxxxx
862            neig0 = _mm_load_pd(coeff[11] + x - offsetZ - offsetY + 0);
863            neig1 = _mm_load_pd(coeff[11] + x - offsetZ - offsetY + 2);
864            neig2 = _mm_load_pd(coeff[11] + x - offsetZ - offsetY + 4);
865            neig3 = _mm_load_pd(coeff[11] + x - offsetZ - offsetY + 6);
866
867            same0 = _mm_add_pd(same0, neig0);
868            same1 = _mm_add_pd(same1, neig1);
869            same2 = _mm_add_pd(same2, neig2);
870            same3 = _mm_add_pd(same3, neig3);
871
872            //xxxxxxxxxxxxx
873            neig0 = _mm_load_pd(coeff[11] + x - offsetZ + 0);
874            neig1 = _mm_load_pd(coeff[11] + x - offsetZ + 2);
875            neig2 = _mm_load_pd(coeff[11] + x - offsetZ + 4);
876            neig3 = _mm_load_pd(coeff[11] + x - offsetZ + 6);
877
878            same0 = _mm_add_pd(same0, neig0);
879            same1 = _mm_add_pd(same1, neig1);
880            same2 = _mm_add_pd(same2, neig2);
881            same3 = _mm_add_pd(same3, neig3);
882
883            //xxxxxxxxxxxxx
884            neig0 = _mm_load_pd(coeff[11] + x - offsetZ + offsetY + 0);
885            neig1 = _mm_load_pd(coeff[11] + x - offsetZ + offsetY + 2);
886            neig2 = _mm_load_pd(coeff[11] + x - offsetZ + offsetY + 4);
887            neig3 = _mm_load_pd(coeff[11] + x - offsetZ + offsetY + 6);
888
889            same0 = _mm_add_pd(same0, neig0);
890            same1 = _mm_add_pd(same1, neig1);
891            same2 = _mm_add_pd(same2, neig2);
892            same3 = _mm_add_pd(same3, neig3);
893
894            //xxxxxxxxxxxxx
895            neig0 = _mm_load_pd(coeff[11] + x - offsetY + 0);
896            neig1 = _mm_load_pd(coeff[11] + x - offsetY + 2);
897            neig2 = _mm_load_pd(coeff[11] + x - offsetY + 4);
898            neig3 = _mm_load_pd(coeff[11] + x - offsetY + 6);
899
900            same0 = _mm_add_pd(same0, neig0);
901            same1 = _mm_add_pd(same1, neig1);
902            same2 = _mm_add_pd(same2, neig2);
903            same3 = _mm_add_pd(same3, neig3);
904
905            //xxxxxxxxxxxxx
906            neig0 = _mm_load_pd(coeff[11] + offsetY + x + 0);
907            neig1 = _mm_load_pd(coeff[11] + offsetY + x + 2);
908            neig2 = _mm_load_pd(coeff[11] + offsetY + x + 4);
909            neig3 = _mm_load_pd(coeff[11] + offsetY + x + 6);
910
911            same0 = _mm_add_pd(same0, neig0);
912            same1 = _mm_add_pd(same1, neig1);
913            same2 = _mm_add_pd(same2, neig2);
914            same3 = _mm_add_pd(same3, neig3);
915
916            //xxxxxxxxxxxxx
917            neig0 = _mm_load_pd(coeff[11] + x + 0);
918            neig1 = _mm_load_pd(coeff[11] + x + 2);
919            neig2 = _mm_load_pd(coeff[11] + x + 4);
920            neig3 = _mm_load_pd(coeff[11] + x + 6);
921
922            same0 = _mm_add_pd(same0, neig0);
923            same1 = _mm_add_pd(same1, neig1);
924            same2 = _mm_add_pd(same2, neig2);
925            same3 = _mm_add_pd(same3, neig3);
926
927            //xxxxxxxxxxxxx
928            neig0 = _mm_load_pd(coeff[11] + x + offsetZ - offsetY + 0);
929            neig1 = _mm_load_pd(coeff[11] + x + offsetZ - offsetY + 2);
930            neig2 = _mm_load_pd(coeff[11] + x + offsetZ - offsetY + 4);
931            neig3 = _mm_load_pd(coeff[11] + x + offsetZ - offsetY + 6);
932
933            same0 = _mm_add_pd(same0, neig0);
934            same1 = _mm_add_pd(same1, neig1);
935            same2 = _mm_add_pd(same2, neig2);
936            same3 = _mm_add_pd(same3, neig3);
937
938            //xxxxxxxxxxxxx
939            neig0 = _mm_load_pd(coeff[11] + x + offsetZ + 0);
940            neig1 = _mm_load_pd(coeff[11] + x + offsetZ + 2);
941            neig2 = _mm_load_pd(coeff[11] + x + offsetZ + 4);
942            neig3 = _mm_load_pd(coeff[11] + x + offsetZ + 6);
943
944            same0 = _mm_add_pd(same0, neig0);
945            same1 = _mm_add_pd(same1, neig1);
946            same2 = _mm_add_pd(same2, neig2);
947            same3 = _mm_add_pd(same3, neig3);
948
949            //xxxxxxxxxxxxx
950            neig0 = _mm_load_pd(coeff[11] + x + offsetZ - offsetY + 0);
951            neig1 = _mm_load_pd(coeff[11] + x + offsetZ - offsetY + 2);
952            neig2 = _mm_load_pd(coeff[11] + x + offsetZ - offsetY + 4);
953            neig3 = _mm_load_pd(coeff[11] + x + offsetZ - offsetY + 6);
954
955            same0 = _mm_add_pd(same0, neig0);
956            same1 = _mm_add_pd(same1, neig1);
957            same2 = _mm_add_pd(same2, neig2);
958            same3 = _mm_add_pd(same3, neig3);
959
960            //xxxxxxxxxxxxx
961            neig0 = _mm_load_pd(coeff[11] + x + offsetZ + offsetY + 0);
962            neig1 = _mm_load_pd(coeff[11] + x + offsetZ + offsetY + 2);
963            neig2 = _mm_load_pd(coeff[11] + x + offsetZ + offsetY + 4);
964            neig3 = _mm_load_pd(coeff[11] + x + offsetZ + offsetY + 6);
965
966            same0 = _mm_add_pd(same0, neig0);
967            same1 = _mm_add_pd(same1, neig1);
968            same2 = _mm_add_pd(same2, neig2);
969            same3 = _mm_add_pd(same3, neig3);
970
971            //xxxxxxxxxxxxx
972            neig0 = _mm_load_pd(coeff[12] + x - offsetZ - offsetY + 0);
973            neig1 = _mm_load_pd(coeff[12] + x - offsetZ - offsetY + 2);
974            neig2 = _mm_load_pd(coeff[12] + x - offsetZ - offsetY + 4);
975            neig3 = _mm_load_pd(coeff[12] + x - offsetZ - offsetY + 6);
976
977            same0 = _mm_add_pd(same0, neig0);
978            same1 = _mm_add_pd(same1, neig1);
979            same2 = _mm_add_pd(same2, neig2);
980            same3 = _mm_add_pd(same3, neig3);
981
982            //xxxxxxxxxxxxx
983            neig0 = _mm_load_pd(coeff[12] + x - offsetZ + 0);
984            neig1 = _mm_load_pd(coeff[12] + x - offsetZ + 2);
985            neig2 = _mm_load_pd(coeff[12] + x - offsetZ + 4);
986            neig3 = _mm_load_pd(coeff[12] + x - offsetZ + 6);
987
988            same0 = _mm_add_pd(same0, neig0);
989            same1 = _mm_add_pd(same1, neig1);
990            same2 = _mm_add_pd(same2, neig2);
991            same3 = _mm_add_pd(same3, neig3);
992
993            //xxxxxxxxxxxxx
994            neig0 = _mm_load_pd(coeff[12] + x - offsetZ + offsetY + 0);
995            neig1 = _mm_load_pd(coeff[12] + x - offsetZ + offsetY + 2);
996            neig2 = _mm_load_pd(coeff[12] + x - offsetZ + offsetY + 4);
997            neig3 = _mm_load_pd(coeff[12] + x - offsetZ + offsetY + 6);
998
999            same0 = _mm_add_pd(same0, neig0);
1000            same1 = _mm_add_pd(same1, neig1);
1001            same2 = _mm_add_pd(same2, neig2);
1002            same3 = _mm_add_pd(same3, neig3);
1003
1004            //xxxxxxxxxxxxx
1005            neig0 = _mm_load_pd(coeff[12] + x - offsetY + 0);
1006            neig1 = _mm_load_pd(coeff[12] + x - offsetY + 2);
1007            neig2 = _mm_load_pd(coeff[12] + x - offsetY + 4);
1008            neig3 = _mm_load_pd(coeff[12] + x - offsetY + 6);
1009
1010            same0 = _mm_add_pd(same0, neig0);
1011            same1 = _mm_add_pd(same1, neig1);
1012            same2 = _mm_add_pd(same2, neig2);
1013            same3 = _mm_add_pd(same3, neig3);
1014
1015            //xxxxxxxxxxxxx
1016            neig0 = _mm_load_pd(coeff[12] + offsetY + x + 0);
1017            neig1 = _mm_load_pd(coeff[12] + offsetY + x + 2);
1018            neig2 = _mm_load_pd(coeff[12] + offsetY + x + 4);
1019            neig3 = _mm_load_pd(coeff[12] + offsetY + x + 6);
1020
1021            same0 = _mm_add_pd(same0, neig0);
1022            same1 = _mm_add_pd(same1, neig1);
1023            same2 = _mm_add_pd(same2, neig2);
1024            same3 = _mm_add_pd(same3, neig3);
1025
1026            //xxxxxxxxxxxxx
1027            neig0 = _mm_load_pd(coeff[12] + x + 0);
1028            neig1 = _mm_load_pd(coeff[12] + x + 2);
1029            neig2 = _mm_load_pd(coeff[12] + x + 4);
1030            neig3 = _mm_load_pd(coeff[12] + x + 6);
1031
1032            same0 = _mm_add_pd(same0, neig0);
1033            same1 = _mm_add_pd(same1, neig1);
1034            same2 = _mm_add_pd(same2, neig2);
1035            same3 = _mm_add_pd(same3, neig3);
1036
1037            //xxxxxxxxxxxxx
1038            neig0 = _mm_load_pd(coeff[12] + x + offsetZ - offsetY + 0);
1039            neig1 = _mm_load_pd(coeff[12] + x + offsetZ - offsetY + 2);
1040            neig2 = _mm_load_pd(coeff[12] + x + offsetZ - offsetY + 4);
1041            neig3 = _mm_load_pd(coeff[12] + x + offsetZ - offsetY + 6);
1042
1043            same0 = _mm_add_pd(same0, neig0);
1044            same1 = _mm_add_pd(same1, neig1);
1045            same2 = _mm_add_pd(same2, neig2);
1046            same3 = _mm_add_pd(same3, neig3);
1047
1048            //xxxxxxxxxxxxx
1049            neig0 = _mm_load_pd(coeff[12] + x + offsetZ + 0);
1050            neig1 = _mm_load_pd(coeff[12] + x + offsetZ + 2);
1051            neig2 = _mm_load_pd(coeff[12] + x + offsetZ + 4);
1052            neig3 = _mm_load_pd(coeff[12] + x + offsetZ + 6);
1053
1054            same0 = _mm_add_pd(same0, neig0);
1055            same1 = _mm_add_pd(same1, neig1);
1056            same2 = _mm_add_pd(same2, neig2);
1057            same3 = _mm_add_pd(same3, neig3);
1058
1059            //xxxxxxxxxxxxx
1060            neig0 = _mm_load_pd(coeff[12] + x + offsetZ + offsetY + 0);
1061            neig1 = _mm_load_pd(coeff[12] + x + offsetZ + offsetY + 2);
1062            neig2 = _mm_load_pd(coeff[12] + x + offsetZ + offsetY + 4);
1063            neig3 = _mm_load_pd(coeff[12] + x + offsetZ + offsetY + 6);
1064
1065            same0 = _mm_add_pd(same0, neig0);
1066            same1 = _mm_add_pd(same1, neig1);
1067            same2 = _mm_add_pd(same2, neig2);
1068            same3 = _mm_add_pd(same3, neig3);
1069
1070            //yyyyyyyyyyyyy
1071            _mm_store_pd(dst + x + 0, same0);
1072            _mm_store_pd(dst + x + 2, same1);
1073            _mm_store_pd(dst + x + 4, same2);
1074            _mm_store_pd(dst + x + 6, same3);
1075
1076            same0 = same4;
1077            neig0 = neig4;
1078
1079            // dst[x] =
1080            //     coeff[0][x] * src[x - offsetZ] +
1081            //     coeff[1][x] * src[x - offsetY] +
1082            //     coeff[2][x] * src[x - 1] +
1083            //     coeff[3][x] * src[x] +
1084            //     coeff[4][x] * src[x + 1] +
1085            //     coeff[5][x] * src[x + offsetY] +
1086            //     coeff[6][x] * src[x + offsetZ];
1087        }
1088
1089        scalarUpdater.step(coeff, src, dst, offsetY, offsetZ, x, endX);
1090    }
1091
1092    int flops()
1093    {
1094        return 40;
1095    }
1096};
1097
1098template<int DIM_X, int DIM_Y, int DIM_Z>
1099class ExtendedVectorized3DFixed
1100{
1101public:
1102    static int coefficients()
1103    {
1104        return 13;
1105    }
1106
1107    inline void step(double *coeff[13], double *src, double *dst, int unusedOffsetY, int unusedOffsetZ, int startX, int endX)
1108    {
1109        const int SLICE_SIZE = DIM_X * DIM_Y;
1110        const int TOTAL_SIZE = DIM_X * DIM_Y * DIM_Z;
1111        const int offsetY = DIM_X;
1112        const int offsetZ = SLICE_SIZE;
1113
1114        int x = startX;
1115        ExtendedScalar3D scalarUpdater;
1116
1117        if ((x & 1) == 1) {
1118            scalarUpdater.step(coeff, src, dst, DIM_Y, offsetZ, x, x + 1);
1119            x += 1;
1120        }
1121
1122        __m128d same0 = _mm_load_pd(src + x + 0);
1123        __m128d neig0 = _mm_loadu_pd(src + x + 1);
1124       
1125        int paddedEndX = endX - 7;
1126        for (; x < paddedEndX; x += 8) {
1127            __m128d same1 = _mm_load_pd(src + x + 2);
1128            __m128d same2 = _mm_load_pd(src + x + 4);
1129            __m128d same3 = _mm_load_pd(src + x + 6);
1130            __m128d same4 = _mm_load_pd(src + x + 8);
1131
1132            __m128d neig1 = _mm_shuffle_pd(same0, same1, (1 << 0) | (0 << 2));
1133            __m128d neig2 = _mm_shuffle_pd(same1, same2, (1 << 0) | (0 << 2));
1134            __m128d neig3 = _mm_shuffle_pd(same2, same3, (1 << 0) | (0 << 2));
1135            __m128d neig4 = _mm_shuffle_pd(same3, same4, (1 << 0) | (0 << 2));
1136
1137            same0 = _mm_mul_pd(same0, _mm_load_pd(&coeff[0][TOTAL_SIZE * 3 + x + 0]));
1138            same1 = _mm_mul_pd(same1, _mm_load_pd(&coeff[0][TOTAL_SIZE * 3 + x + 2]));
1139            same2 = _mm_mul_pd(same2, _mm_load_pd(&coeff[0][TOTAL_SIZE * 3 + x + 4]));
1140            same3 = _mm_mul_pd(same3, _mm_load_pd(&coeff[0][TOTAL_SIZE * 3 + x + 6]));
1141
1142            __m128d temp1 = _mm_mul_pd(neig0, _mm_load_pd(&coeff[0][TOTAL_SIZE * 2 + x + 0]));
1143            __m128d temp2 = _mm_mul_pd(neig1, _mm_load_pd(&coeff[0][TOTAL_SIZE * 2 + x + 2]));
1144            __m128d temp3 = _mm_mul_pd(neig2, _mm_load_pd(&coeff[0][TOTAL_SIZE * 2 + x + 4]));
1145            __m128d temp4 = _mm_mul_pd(neig3, _mm_load_pd(&coeff[0][TOTAL_SIZE * 2 + x + 6]));
1146
1147            same0 = _mm_add_pd(same0, temp1);
1148            same1 = _mm_add_pd(same1, temp2);
1149            same2 = _mm_add_pd(same2, temp3);
1150            same3 = _mm_add_pd(same3, temp4);
1151
1152            temp1 = _mm_mul_pd(neig1, _mm_load_pd(&coeff[0][TOTAL_SIZE * 4 + x + 0]));
1153            temp2 = _mm_mul_pd(neig2, _mm_load_pd(&coeff[0][TOTAL_SIZE * 4 + x + 2]));
1154            temp3 = _mm_mul_pd(neig3, _mm_load_pd(&coeff[0][TOTAL_SIZE * 4 + x + 4]));
1155            temp4 = _mm_mul_pd(neig4, _mm_load_pd(&coeff[0][TOTAL_SIZE * 4 + x + 6]));
1156
1157            same0 = _mm_add_pd(same0, temp1);
1158            same1 = _mm_add_pd(same1, temp2);
1159            same2 = _mm_add_pd(same2, temp3);
1160            same3 = _mm_add_pd(same3, temp4);
1161
1162            neig0 = _mm_load_pd(src + x - offsetZ + 0);
1163            neig1 = _mm_load_pd(src + x - offsetZ + 2);
1164            neig2 = _mm_load_pd(src + x - offsetZ + 4);
1165            neig3 = _mm_load_pd(src + x - offsetZ + 6);
1166
1167            temp1 = _mm_mul_pd(neig0, _mm_load_pd(&coeff[0][TOTAL_SIZE * 0 + x + 0]));
1168            temp2 = _mm_mul_pd(neig1, _mm_load_pd(&coeff[0][TOTAL_SIZE * 0 + x + 2]));
1169            temp3 = _mm_mul_pd(neig2, _mm_load_pd(&coeff[0][TOTAL_SIZE * 0 + x + 4]));
1170            temp4 = _mm_mul_pd(neig3, _mm_load_pd(&coeff[0][TOTAL_SIZE * 0 + x + 6]));
1171
1172            same0 = _mm_add_pd(same0, temp1);
1173            same1 = _mm_add_pd(same1, temp2);
1174            same2 = _mm_add_pd(same2, temp3);
1175            same3 = _mm_add_pd(same3, temp4);
1176
1177            neig0 = _mm_load_pd(src + x - offsetY + 0);
1178            neig1 = _mm_load_pd(src + x - offsetY + 2);
1179            neig2 = _mm_load_pd(src + x - offsetY + 4);
1180            neig3 = _mm_load_pd(src + x - offsetY + 6);
1181
1182            temp1 = _mm_mul_pd(neig0, _mm_load_pd(&coeff[0][TOTAL_SIZE * 1 + x + 0]));
1183            temp2 = _mm_mul_pd(neig1, _mm_load_pd(&coeff[0][TOTAL_SIZE * 1 + x + 2]));
1184            temp3 = _mm_mul_pd(neig2, _mm_load_pd(&coeff[0][TOTAL_SIZE * 1 + x + 4]));
1185            temp4 = _mm_mul_pd(neig3, _mm_load_pd(&coeff[0][TOTAL_SIZE * 1 + x + 6]));
1186
1187            same0 = _mm_add_pd(same0, temp1);
1188            same1 = _mm_add_pd(same1, temp2);
1189            same2 = _mm_add_pd(same2, temp3);
1190            same3 = _mm_add_pd(same3, temp4);
1191
1192            neig0 = _mm_load_pd(src + x + offsetY + 0);
1193            neig1 = _mm_load_pd(src + x + offsetY + 2);
1194            neig2 = _mm_load_pd(src + x + offsetY + 4);
1195            neig3 = _mm_load_pd(src + x + offsetY + 6);
1196
1197            temp1 = _mm_mul_pd(neig0, _mm_load_pd(&coeff[0][TOTAL_SIZE * 5 + x + 0]));
1198            temp2 = _mm_mul_pd(neig1, _mm_load_pd(&coeff[0][TOTAL_SIZE * 5 + x + 2]));
1199            temp3 = _mm_mul_pd(neig2, _mm_load_pd(&coeff[0][TOTAL_SIZE * 5 + x + 4]));
1200            temp4 = _mm_mul_pd(neig3, _mm_load_pd(&coeff[0][TOTAL_SIZE * 5 + x + 6]));
1201
1202            same0 = _mm_add_pd(same0, temp1);
1203            same1 = _mm_add_pd(same1, temp2);
1204            same2 = _mm_add_pd(same2, temp3);
1205            same3 = _mm_add_pd(same3, temp4);
1206
1207            //xxxxxxxxxxxxx
1208            neig0 = _mm_load_pd(src + x + offsetZ + 0);
1209            neig1 = _mm_load_pd(src + x + offsetZ + 2);
1210            neig2 = _mm_load_pd(src + x + offsetZ + 4);
1211            neig3 = _mm_load_pd(src + x + offsetZ + 6);
1212
1213            temp1 = _mm_mul_pd(neig0, _mm_load_pd(&coeff[0][TOTAL_SIZE * 6 + x + 0]));
1214            temp2 = _mm_mul_pd(neig1, _mm_load_pd(&coeff[0][TOTAL_SIZE * 6 + x + 2]));
1215            temp3 = _mm_mul_pd(neig2, _mm_load_pd(&coeff[0][TOTAL_SIZE * 6 + x + 4]));
1216            temp4 = _mm_mul_pd(neig3, _mm_load_pd(&coeff[0][TOTAL_SIZE * 6 + x + 6]));
1217
1218            same0 = _mm_add_pd(same0, temp1);
1219            same1 = _mm_add_pd(same1, temp2);
1220            same2 = _mm_add_pd(same2, temp3);
1221            same3 = _mm_add_pd(same3, temp4);
1222
1223            //xxxxxxxxxxxxx
1224            neig0 = _mm_load_pd(src + x - offsetZ - offsetY + 0);
1225            neig1 = _mm_load_pd(src + x - offsetZ - offsetY + 2);
1226            neig2 = _mm_load_pd(src + x - offsetZ - offsetY + 4);
1227            neig3 = _mm_load_pd(src + x - offsetZ - offsetY + 6);
1228
1229            temp1 = _mm_mul_pd(neig0, _mm_load_pd(&coeff[0][TOTAL_SIZE * 7 + x + 0]));
1230            temp2 = _mm_mul_pd(neig1, _mm_load_pd(&coeff[0][TOTAL_SIZE * 7 + x + 2]));
1231            temp3 = _mm_mul_pd(neig2, _mm_load_pd(&coeff[0][TOTAL_SIZE * 7 + x + 4]));
1232            temp4 = _mm_mul_pd(neig3, _mm_load_pd(&coeff[0][TOTAL_SIZE * 7 + x + 6]));
1233
1234            same0 = _mm_add_pd(same0, temp1);
1235            same1 = _mm_add_pd(same1, temp2);
1236            same2 = _mm_add_pd(same2, temp3);
1237            same3 = _mm_add_pd(same3, temp4);
1238
1239            //xxxxxxxxxxxxx
1240            neig0 = _mm_load_pd(src + x - offsetZ + offsetY + 0);
1241            neig1 = _mm_load_pd(src + x - offsetZ + offsetY + 2);
1242            neig2 = _mm_load_pd(src + x - offsetZ + offsetY + 4);
1243            neig3 = _mm_load_pd(src + x - offsetZ + offsetY + 6);
1244
1245            temp1 = _mm_mul_pd(neig0, _mm_load_pd(&coeff[0][TOTAL_SIZE * 8 + x + 0]));
1246            temp2 = _mm_mul_pd(neig1, _mm_load_pd(&coeff[0][TOTAL_SIZE * 8 + x + 2]));
1247            temp3 = _mm_mul_pd(neig2, _mm_load_pd(&coeff[0][TOTAL_SIZE * 8 + x + 4]));
1248            temp4 = _mm_mul_pd(neig3, _mm_load_pd(&coeff[0][TOTAL_SIZE * 8 + x + 6]));
1249
1250            same0 = _mm_add_pd(same0, temp1);
1251            same1 = _mm_add_pd(same1, temp2);
1252            same2 = _mm_add_pd(same2, temp3);
1253            same3 = _mm_add_pd(same3, temp4);
1254
1255            //xxxxxxxxxxxxx
1256            neig0 = _mm_load_pd(src + x + offsetZ - offsetY + 0);
1257            neig1 = _mm_load_pd(src + x + offsetZ - offsetY + 2);
1258            neig2 = _mm_load_pd(src + x + offsetZ - offsetY + 4);
1259            neig3 = _mm_load_pd(src + x + offsetZ - offsetY + 6);
1260
1261            temp1 = _mm_mul_pd(neig0, _mm_load_pd(&coeff[0][TOTAL_SIZE * 9 + x + 0]));
1262            temp2 = _mm_mul_pd(neig1, _mm_load_pd(&coeff[0][TOTAL_SIZE * 9 + x + 2]));
1263            temp3 = _mm_mul_pd(neig2, _mm_load_pd(&coeff[0][TOTAL_SIZE * 9 + x + 4]));
1264            temp4 = _mm_mul_pd(neig3, _mm_load_pd(&coeff[0][TOTAL_SIZE * 9 + x + 6]));
1265
1266            same0 = _mm_add_pd(same0, temp1);
1267            same1 = _mm_add_pd(same1, temp2);
1268            same2 = _mm_add_pd(same2, temp3);
1269            same3 = _mm_add_pd(same3, temp4);
1270
1271            //xxxxxxxxxxxxx
1272            neig0 = _mm_load_pd(src + x + offsetZ + offsetY + 0);
1273            neig1 = _mm_load_pd(src + x + offsetZ + offsetY + 2);
1274            neig2 = _mm_load_pd(src + x + offsetZ + offsetY + 4);
1275            neig3 = _mm_load_pd(src + x + offsetZ + offsetY + 6);
1276
1277            temp1 = _mm_mul_pd(neig0, _mm_load_pd(&coeff[0][TOTAL_SIZE * 10 + x + 0]));
1278            temp2 = _mm_mul_pd(neig1, _mm_load_pd(&coeff[0][TOTAL_SIZE * 10 + x + 2]));
1279            temp3 = _mm_mul_pd(neig2, _mm_load_pd(&coeff[0][TOTAL_SIZE * 10 + x + 4]));
1280            temp4 = _mm_mul_pd(neig3, _mm_load_pd(&coeff[0][TOTAL_SIZE * 10 + x + 6]));
1281
1282            same0 = _mm_add_pd(same0, temp1);
1283            same1 = _mm_add_pd(same1, temp2);
1284            same2 = _mm_add_pd(same2, temp3);
1285            same3 = _mm_add_pd(same3, temp4);
1286
1287            //xxxxxxxxxxxxx
1288            neig0 = _mm_load_pd(coeff[0] + 11 * TOTAL_SIZE + x - offsetZ - offsetY + 0);
1289            neig1 = _mm_load_pd(coeff[0] + 11 * TOTAL_SIZE + x - offsetZ - offsetY + 2);
1290            neig2 = _mm_load_pd(coeff[0] + 11 * TOTAL_SIZE + x - offsetZ - offsetY + 4);
1291            neig3 = _mm_load_pd(coeff[0] + 11 * TOTAL_SIZE + x - offsetZ - offsetY + 6);
1292
1293            same0 = _mm_add_pd(same0, neig0);
1294            same1 = _mm_add_pd(same1, neig1);
1295            same2 = _mm_add_pd(same2, neig2);
1296            same3 = _mm_add_pd(same3, neig3);
1297
1298            //xxxxxxxxxxxxx
1299            neig0 = _mm_load_pd(coeff[0] + 11 * TOTAL_SIZE + x - offsetZ + 0);
1300            neig1 = _mm_load_pd(coeff[0] + 11 * TOTAL_SIZE + x - offsetZ + 2);
1301            neig2 = _mm_load_pd(coeff[0] + 11 * TOTAL_SIZE + x - offsetZ + 4);
1302            neig3 = _mm_load_pd(coeff[0] + 11 * TOTAL_SIZE + x - offsetZ + 6);
1303
1304            same0 = _mm_add_pd(same0, neig0);
1305            same1 = _mm_add_pd(same1, neig1);
1306            same2 = _mm_add_pd(same2, neig2);
1307            same3 = _mm_add_pd(same3, neig3);
1308
1309            //xxxxxxxxxxxxx
1310            neig0 = _mm_load_pd(coeff[0] + 11 * TOTAL_SIZE + x - offsetZ + offsetY + 0);
1311            neig1 = _mm_load_pd(coeff[0] + 11 * TOTAL_SIZE + x - offsetZ + offsetY + 2);
1312            neig2 = _mm_load_pd(coeff[0] + 11 * TOTAL_SIZE + x - offsetZ + offsetY + 4);
1313            neig3 = _mm_load_pd(coeff[0] + 11 * TOTAL_SIZE + x - offsetZ + offsetY + 6);
1314
1315            same0 = _mm_add_pd(same0, neig0);
1316            same1 = _mm_add_pd(same1, neig1);
1317            same2 = _mm_add_pd(same2, neig2);
1318            same3 = _mm_add_pd(same3, neig3);
1319
1320            //xxxxxxxxxxxxx
1321            neig0 = _mm_load_pd(coeff[0] + 11 * TOTAL_SIZE + x - offsetY + 0);
1322            neig1 = _mm_load_pd(coeff[0] + 11 * TOTAL_SIZE + x - offsetY + 2);
1323            neig2 = _mm_load_pd(coeff[0] + 11 * TOTAL_SIZE + x - offsetY + 4);
1324            neig3 = _mm_load_pd(coeff[0] + 11 * TOTAL_SIZE + x - offsetY + 6);
1325
1326            same0 = _mm_add_pd(same0, neig0);
1327            same1 = _mm_add_pd(same1, neig1);
1328            same2 = _mm_add_pd(same2, neig2);
1329            same3 = _mm_add_pd(same3, neig3);
1330
1331            //xxxxxxxxxxxxx
1332            neig0 = _mm_load_pd(coeff[0] + 11 * TOTAL_SIZE + offsetY + x + 0);
1333            neig1 = _mm_load_pd(coeff[0] + 11 * TOTAL_SIZE + offsetY + x + 2);
1334            neig2 = _mm_load_pd(coeff[0] + 11 * TOTAL_SIZE + offsetY + x + 4);
1335            neig3 = _mm_load_pd(coeff[0] + 11 * TOTAL_SIZE + offsetY + x + 6);
1336
1337            same0 = _mm_add_pd(same0, neig0);
1338            same1 = _mm_add_pd(same1, neig1);
1339            same2 = _mm_add_pd(same2, neig2);
1340            same3 = _mm_add_pd(same3, neig3);
1341
1342            //xxxxxxxxxxxxx
1343            neig0 = _mm_load_pd(coeff[0] + 11 * TOTAL_SIZE + x + 0);
1344            neig1 = _mm_load_pd(coeff[0] + 11 * TOTAL_SIZE + x + 2);
1345            neig2 = _mm_load_pd(coeff[0] + 11 * TOTAL_SIZE + x + 4);
1346            neig3 = _mm_load_pd(coeff[0] + 11 * TOTAL_SIZE + x + 6);
1347
1348            same0 = _mm_add_pd(same0, neig0);
1349            same1 = _mm_add_pd(same1, neig1);
1350            same2 = _mm_add_pd(same2, neig2);
1351            same3 = _mm_add_pd(same3, neig3);
1352
1353            //xxxxxxxxxxxxx
1354            neig0 = _mm_load_pd(coeff[0] + 11 * TOTAL_SIZE + x + offsetZ - offsetY + 0);
1355            neig1 = _mm_load_pd(coeff[0] + 11 * TOTAL_SIZE + x + offsetZ - offsetY + 2);
1356            neig2 = _mm_load_pd(coeff[0] + 11 * TOTAL_SIZE + x + offsetZ - offsetY + 4);
1357            neig3 = _mm_load_pd(coeff[0] + 11 * TOTAL_SIZE + x + offsetZ - offsetY + 6);
1358
1359            same0 = _mm_add_pd(same0, neig0);
1360            same1 = _mm_add_pd(same1, neig1);
1361            same2 = _mm_add_pd(same2, neig2);
1362            same3 = _mm_add_pd(same3, neig3);
1363
1364            //xxxxxxxxxxxxx
1365            neig0 = _mm_load_pd(coeff[0] + 11 * TOTAL_SIZE + x + offsetZ + 0);
1366            neig1 = _mm_load_pd(coeff[0] + 11 * TOTAL_SIZE + x + offsetZ + 2);
1367            neig2 = _mm_load_pd(coeff[0] + 11 * TOTAL_SIZE + x + offsetZ + 4);
1368            neig3 = _mm_load_pd(coeff[0] + 11 * TOTAL_SIZE + x + offsetZ + 6);
1369
1370            same0 = _mm_add_pd(same0, neig0);
1371            same1 = _mm_add_pd(same1, neig1);
1372            same2 = _mm_add_pd(same2, neig2);
1373            same3 = _mm_add_pd(same3, neig3);
1374
1375            //xxxxxxxxxxxxx
1376            neig0 = _mm_load_pd(coeff[0] + 11 * TOTAL_SIZE + x + offsetZ + offsetY + 0);
1377            neig1 = _mm_load_pd(coeff[0] + 11 * TOTAL_SIZE + x + offsetZ + offsetY + 2);
1378            neig2 = _mm_load_pd(coeff[0] + 11 * TOTAL_SIZE + x + offsetZ + offsetY + 4);
1379            neig3 = _mm_load_pd(coeff[0] + 11 * TOTAL_SIZE + x + offsetZ + offsetY + 6);
1380
1381            same0 = _mm_add_pd(same0, neig0);
1382            same1 = _mm_add_pd(same1, neig1);
1383            same2 = _mm_add_pd(same2, neig2);
1384            same3 = _mm_add_pd(same3, neig3);
1385
1386            //xxxxxxxxxxxxx
1387            neig0 = _mm_load_pd(coeff[0] + 12 * TOTAL_SIZE + x - offsetZ - offsetY + 0);
1388            neig1 = _mm_load_pd(coeff[0] + 12 * TOTAL_SIZE + x - offsetZ - offsetY + 2);
1389            neig2 = _mm_load_pd(coeff[0] + 12 * TOTAL_SIZE + x - offsetZ - offsetY + 4);
1390            neig3 = _mm_load_pd(coeff[0] + 12 * TOTAL_SIZE + x - offsetZ - offsetY + 6);
1391
1392            same0 = _mm_add_pd(same0, neig0);
1393            same1 = _mm_add_pd(same1, neig1);
1394            same2 = _mm_add_pd(same2, neig2);
1395            same3 = _mm_add_pd(same3, neig3);
1396
1397            //xxxxxxxxxxxxx
1398            neig0 = _mm_load_pd(coeff[0] + 12 * TOTAL_SIZE + x - offsetZ + 0);
1399            neig1 = _mm_load_pd(coeff[0] + 12 * TOTAL_SIZE + x - offsetZ + 2);
1400            neig2 = _mm_load_pd(coeff[0] + 12 * TOTAL_SIZE + x - offsetZ + 4);
1401            neig3 = _mm_load_pd(coeff[0] + 12 * TOTAL_SIZE + x - offsetZ + 6);
1402
1403            same0 = _mm_add_pd(same0, neig0);
1404            same1 = _mm_add_pd(same1, neig1);
1405            same2 = _mm_add_pd(same2, neig2);
1406            same3 = _mm_add_pd(same3, neig3);
1407
1408            //xxxxxxxxxxxxx
1409            neig0 = _mm_load_pd(coeff[0] + 12 * TOTAL_SIZE + x - offsetZ + offsetY + 0);
1410            neig1 = _mm_load_pd(coeff[0] + 12 * TOTAL_SIZE + x - offsetZ + offsetY + 2);
1411            neig2 = _mm_load_pd(coeff[0] + 12 * TOTAL_SIZE + x - offsetZ + offsetY + 4);
1412            neig3 = _mm_load_pd(coeff[0] + 12 * TOTAL_SIZE + x - offsetZ + offsetY + 6);
1413
1414            same0 = _mm_add_pd(same0, neig0);
1415            same1 = _mm_add_pd(same1, neig1);
1416            same2 = _mm_add_pd(same2, neig2);
1417            same3 = _mm_add_pd(same3, neig3);
1418
1419            //xxxxxxxxxxxxx
1420            neig0 = _mm_load_pd(coeff[0] + 12 * TOTAL_SIZE + x - offsetY + 0);
1421            neig1 = _mm_load_pd(coeff[0] + 12 * TOTAL_SIZE + x - offsetY + 2);
1422            neig2 = _mm_load_pd(coeff[0] + 12 * TOTAL_SIZE + x - offsetY + 4);
1423            neig3 = _mm_load_pd(coeff[0] + 12 * TOTAL_SIZE + x - offsetY + 6);
1424
1425            same0 = _mm_add_pd(same0, neig0);
1426            same1 = _mm_add_pd(same1, neig1);
1427            same2 = _mm_add_pd(same2, neig2);
1428            same3 = _mm_add_pd(same3, neig3);
1429
1430            //xxxxxxxxxxxxx
1431            neig0 = _mm_load_pd(coeff[0] + 12 * TOTAL_SIZE + offsetY + x + 0);
1432            neig1 = _mm_load_pd(coeff[0] + 12 * TOTAL_SIZE + offsetY + x + 2);
1433            neig2 = _mm_load_pd(coeff[0] + 12 * TOTAL_SIZE + offsetY + x + 4);
1434            neig3 = _mm_load_pd(coeff[0] + 12 * TOTAL_SIZE + offsetY + x + 6);
1435
1436            same0 = _mm_add_pd(same0, neig0);
1437            same1 = _mm_add_pd(same1, neig1);
1438            same2 = _mm_add_pd(same2, neig2);
1439            same3 = _mm_add_pd(same3, neig3);
1440
1441            //xxxxxxxxxxxxx
1442            neig0 = _mm_load_pd(coeff[0] + 12 * TOTAL_SIZE + x + 0);
1443            neig1 = _mm_load_pd(coeff[0] + 12 * TOTAL_SIZE + x + 2);
1444            neig2 = _mm_load_pd(coeff[0] + 12 * TOTAL_SIZE + x + 4);
1445            neig3 = _mm_load_pd(coeff[0] + 12 * TOTAL_SIZE + x + 6);
1446
1447            same0 = _mm_add_pd(same0, neig0);
1448            same1 = _mm_add_pd(same1, neig1);
1449            same2 = _mm_add_pd(same2, neig2);
1450            same3 = _mm_add_pd(same3, neig3);
1451
1452            //xxxxxxxxxxxxx
1453            neig0 = _mm_load_pd(coeff[0] + 12 * TOTAL_SIZE + x + offsetZ - offsetY + 0);
1454            neig1 = _mm_load_pd(coeff[0] + 12 * TOTAL_SIZE + x + offsetZ - offsetY + 2);
1455            neig2 = _mm_load_pd(coeff[0] + 12 * TOTAL_SIZE + x + offsetZ - offsetY + 4);
1456            neig3 = _mm_load_pd(coeff[0] + 12 * TOTAL_SIZE + x + offsetZ - offsetY + 6);
1457
1458            same0 = _mm_add_pd(same0, neig0);
1459            same1 = _mm_add_pd(same1, neig1);
1460            same2 = _mm_add_pd(same2, neig2);
1461            same3 = _mm_add_pd(same3, neig3);
1462
1463            //xxxxxxxxxxxxx
1464            neig0 = _mm_load_pd(coeff[0] + 12 * TOTAL_SIZE + x + offsetZ + 0);
1465            neig1 = _mm_load_pd(coeff[0] + 12 * TOTAL_SIZE + x + offsetZ + 2);
1466            neig2 = _mm_load_pd(coeff[0] + 12 * TOTAL_SIZE + x + offsetZ + 4);
1467            neig3 = _mm_load_pd(coeff[0] + 12 * TOTAL_SIZE + x + offsetZ + 6);
1468
1469            same0 = _mm_add_pd(same0, neig0);
1470            same1 = _mm_add_pd(same1, neig1);
1471            same2 = _mm_add_pd(same2, neig2);
1472            same3 = _mm_add_pd(same3, neig3);
1473
1474            //xxxxxxxxxxxxx
1475            neig0 = _mm_load_pd(coeff[0] + 12 * TOTAL_SIZE + x + offsetZ + offsetY + 0);
1476            neig1 = _mm_load_pd(coeff[0] + 12 * TOTAL_SIZE + x + offsetZ + offsetY + 2);
1477            neig2 = _mm_load_pd(coeff[0] + 12 * TOTAL_SIZE + x + offsetZ + offsetY + 4);
1478            neig3 = _mm_load_pd(coeff[0] + 12 * TOTAL_SIZE + x + offsetZ + offsetY + 6);
1479
1480            same0 = _mm_add_pd(same0, neig0);
1481            same1 = _mm_add_pd(same1, neig1);
1482            same2 = _mm_add_pd(same2, neig2);
1483            same3 = _mm_add_pd(same3, neig3);
1484
1485            //yyyyyyyyyyyyy
1486            _mm_store_pd(dst + x + 0, same0);
1487            _mm_store_pd(dst + x + 2, same1);
1488            _mm_store_pd(dst + x + 4, same2);
1489            _mm_store_pd(dst + x + 6, same3);
1490
1491            same0 = same4;
1492            neig0 = neig4;
1493
1494            // dst[x] =
1495            //     coeff[0][x] * src[x - offsetZ] +
1496            //     coeff[1][x] * src[x - offsetY] +
1497            //     coeff[2][x] * src[x - 1] +
1498            //     coeff[3][x] * src[x] +
1499            //     coeff[4][x] * src[x + 1] +
1500            //     coeff[5][x] * src[x + offsetY] +
1501            //     coeff[6][x] * src[x + offsetZ];
1502        }
1503
1504        scalarUpdater.step(coeff, src, dst, offsetY, offsetZ, x, endX);
1505    }
1506
1507    int flops()
1508    {
1509        return 40;
1510    }
1511};
1512
1513template<int DIM_X, int DIM_Y, int DIM_Z>
1514class ExtendedVectorized3DNextGen
1515{
1516public:
1517    static int coefficients()
1518    {
1519        return 13;
1520    }
1521
1522    template<class NEIGHBORHOOD>
1523    inline void step(const NEIGHBORHOOD& hood, double *dst, int unusedOffsetY, int unusedOffsetZ, int startX, int endX)
1524    {
1525        const int SLICE_SIZE = DIM_X * DIM_Y;
1526        const int TOTAL_SIZE = DIM_X * DIM_Y * DIM_Z;
1527        const int offsetY = DIM_X;
1528        const int offsetZ = SLICE_SIZE;
1529
1530        int x = startX;
1531
1532        if ((x & 1) == 1) {
1533            stepScalar(hood, dst, DIM_Y, offsetZ, x, x + 1);
1534            x += 1;
1535        }
1536
1537        __m128d same0 = _mm_load_pd(&(hood.template src<0, 0, 0>(x)));
1538        __m128d neig0 = _mm_loadu_pd(&hood.template src<1, 0, 0>(x));
1539
1540        int paddedEndX = endX - 7;
1541        for (; x < paddedEndX; x += 8) {
1542            __m128d same1 = _mm_load_pd(&hood.template src<2, 0, 0>(x));
1543            __m128d same2 = _mm_load_pd(&hood.template src<4, 0, 0>(x));
1544            __m128d same3 = _mm_load_pd(&hood.template src<6, 0, 0>(x));
1545            __m128d same4 = _mm_load_pd(&hood.template src<8, 0, 0>(x));
1546
1547            __m128d neig1 = _mm_shuffle_pd(same0, same1, (1 << 0) | (0 << 2));
1548            __m128d neig2 = _mm_shuffle_pd(same1, same2, (1 << 0) | (0 << 2));
1549            __m128d neig3 = _mm_shuffle_pd(same2, same3, (1 << 0) | (0 << 2));
1550            __m128d neig4 = _mm_shuffle_pd(same3, same4, (1 << 0) | (0 << 2));
1551
1552            same0 = _mm_mul_pd(same0, _mm_load_pd(&hood.template coeff<0, 0, 0, 3>(x)));
1553            same1 = _mm_mul_pd(same1, _mm_load_pd(&hood.template coeff<2, 0, 0, 3>(x)));
1554            same2 = _mm_mul_pd(same2, _mm_load_pd(&hood.template coeff<4, 0, 0, 3>(x)));
1555            same3 = _mm_mul_pd(same3, _mm_load_pd(&hood.template coeff<6, 0, 0, 3>(x)));
1556
1557            __m128d temp1 = _mm_mul_pd(neig0, _mm_load_pd(&hood.template coeff<0, 0, 0, 2>(x)));
1558            __m128d temp2 = _mm_mul_pd(neig1, _mm_load_pd(&hood.template coeff<2, 0, 0, 2>(x)));
1559            __m128d temp3 = _mm_mul_pd(neig2, _mm_load_pd(&hood.template coeff<4, 0, 0, 2>(x)));
1560            __m128d temp4 = _mm_mul_pd(neig3, _mm_load_pd(&hood.template coeff<6, 0, 0, 2>(x)));
1561
1562            same0 = _mm_add_pd(same0, temp1);
1563            same1 = _mm_add_pd(same1, temp2);
1564            same2 = _mm_add_pd(same2, temp3);
1565            same3 = _mm_add_pd(same3, temp4);
1566
1567            temp1 = _mm_mul_pd(neig1, _mm_load_pd(&hood.template coeff<0, 0, 0, 4>(x)));
1568            temp2 = _mm_mul_pd(neig2, _mm_load_pd(&hood.template coeff<2, 0, 0, 4>(x)));
1569            temp3 = _mm_mul_pd(neig3, _mm_load_pd(&hood.template coeff<4, 0, 0, 4>(x)));
1570            temp4 = _mm_mul_pd(neig4, _mm_load_pd(&hood.template coeff<6, 0, 0, 4>(x)));
1571
1572            same0 = _mm_add_pd(same0, temp1);
1573            same1 = _mm_add_pd(same1, temp2);
1574            same2 = _mm_add_pd(same2, temp3);
1575            same3 = _mm_add_pd(same3, temp4);
1576
1577            neig0 = _mm_load_pd(&hood.template src<0, 0, -1>(x));
1578            neig1 = _mm_load_pd(&hood.template src<2, 0, -1>(x));
1579            neig2 = _mm_load_pd(&hood.template src<4, 0, -1>(x));
1580            neig3 = _mm_load_pd(&hood.template src<6, 0, -1>(x));
1581
1582            temp1 = _mm_mul_pd(neig0, _mm_load_pd(&hood.template coeff<0, 0, 0, 0>(x)));
1583            temp2 = _mm_mul_pd(neig1, _mm_load_pd(&hood.template coeff<2, 0, 0, 0>(x)));
1584            temp3 = _mm_mul_pd(neig2, _mm_load_pd(&hood.template coeff<4, 0, 0, 0>(x)));
1585            temp4 = _mm_mul_pd(neig3, _mm_load_pd(&hood.template coeff<6, 0, 0, 0>(x)));
1586
1587            same0 = _mm_add_pd(same0, temp1);
1588            same1 = _mm_add_pd(same1, temp2);
1589            same2 = _mm_add_pd(same2, temp3);
1590            same3 = _mm_add_pd(same3, temp4);
1591
1592            neig0 = _mm_load_pd(&hood.template src<0, -1, 0>(x));
1593            neig1 = _mm_load_pd(&hood.template src<2, -1, 0>(x));
1594            neig2 = _mm_load_pd(&hood.template src<4, -1, 0>(x));
1595            neig3 = _mm_load_pd(&hood.template src<6, -1, 0>(x));
1596
1597            temp1 = _mm_mul_pd(neig0, _mm_load_pd(&hood.template coeff<0, 0, 0, 1>(x)));
1598            temp2 = _mm_mul_pd(neig1, _mm_load_pd(&hood.template coeff<2, 0, 0, 1>(x)));
1599            temp3 = _mm_mul_pd(neig2, _mm_load_pd(&hood.template coeff<4, 0, 0, 1>(x)));
1600            temp4 = _mm_mul_pd(neig3, _mm_load_pd(&hood.template coeff<6, 0, 0, 1>(x)));
1601
1602            same0 = _mm_add_pd(same0, temp1);
1603            same1 = _mm_add_pd(same1, temp2);
1604            same2 = _mm_add_pd(same2, temp3);
1605            same3 = _mm_add_pd(same3, temp4);
1606
1607            neig0 = _mm_load_pd(&hood.template src<0, 1, 0>(x));
1608            neig1 = _mm_load_pd(&hood.template src<2, 1, 0>(x));
1609            neig2 = _mm_load_pd(&hood.template src<4, 1, 0>(x));
1610            neig3 = _mm_load_pd(&hood.template src<6, 1, 0>(x));
1611
1612            temp1 = _mm_mul_pd(neig0, _mm_load_pd(&hood.template coeff<0, 0, 0, 5>(x)));
1613            temp2 = _mm_mul_pd(neig1, _mm_load_pd(&hood.template coeff<2, 0, 0, 5>(x)));
1614            temp3 = _mm_mul_pd(neig2, _mm_load_pd(&hood.template coeff<4, 0, 0, 5>(x)));
1615            temp4 = _mm_mul_pd(neig3, _mm_load_pd(&hood.template coeff<6, 0, 0, 5>(x)));
1616
1617            same0 = _mm_add_pd(same0, temp1);
1618            same1 = _mm_add_pd(same1, temp2);
1619            same2 = _mm_add_pd(same2, temp3);
1620            same3 = _mm_add_pd(same3, temp4);
1621
1622            //xxxxxxxxxxxxx
1623            neig0 = _mm_load_pd(&hood.template src<0, 0, 1>(x));
1624            neig1 = _mm_load_pd(&hood.template src<2, 0, 1>(x));
1625            neig2 = _mm_load_pd(&hood.template src<4, 0, 1>(x));
1626            neig3 = _mm_load_pd(&hood.template src<6, 0, 1>(x));
1627
1628            temp1 = _mm_mul_pd(neig0, _mm_load_pd(&hood.template coeff<0, 0, 0, 6>(x)));
1629            temp2 = _mm_mul_pd(neig1, _mm_load_pd(&hood.template coeff<2, 0, 0, 6>(x)));
1630            temp3 = _mm_mul_pd(neig2, _mm_load_pd(&hood.template coeff<4, 0, 0, 6>(x)));
1631            temp4 = _mm_mul_pd(neig3, _mm_load_pd(&hood.template coeff<6, 0, 0, 6>(x)));
1632
1633            same0 = _mm_add_pd(same0, temp1);
1634            same1 = _mm_add_pd(same1, temp2);
1635            same2 = _mm_add_pd(same2, temp3);
1636            same3 = _mm_add_pd(same3, temp4);
1637
1638            //xxxxxxxxxxxxx
1639            neig0 = _mm_load_pd(&hood.template src<0, -1, -1>(x));
1640            neig1 = _mm_load_pd(&hood.template src<2, -1, -1>(x));
1641            neig2 = _mm_load_pd(&hood.template src<4, -1, -1>(x));
1642            neig3 = _mm_load_pd(&hood.template src<6, -1, -1>(x));
1643
1644            temp1 = _mm_mul_pd(neig0, _mm_load_pd(&hood.template coeff<0, 0, 0, 7>(x)));
1645            temp2 = _mm_mul_pd(neig1, _mm_load_pd(&hood.template coeff<2, 0, 0, 7>(x)));
1646            temp3 = _mm_mul_pd(neig2, _mm_load_pd(&hood.template coeff<4, 0, 0, 7>(x)));
1647            temp4 = _mm_mul_pd(neig3, _mm_load_pd(&hood.template coeff<6, 0, 0, 7>(x)));
1648
1649            same0 = _mm_add_pd(same0, temp1);
1650            same1 = _mm_add_pd(same1, temp2);
1651            same2 = _mm_add_pd(same2, temp3);
1652            same3 = _mm_add_pd(same3, temp4);
1653
1654            //xxxxxxxxxxxxx
1655            neig0 = _mm_load_pd(&hood.template src<0, 1, -1>(x));
1656            neig1 = _mm_load_pd(&hood.template src<2, 1, -1>(x));
1657            neig2 = _mm_load_pd(&hood.template src<4, 1, -1>(x));
1658            neig3 = _mm_load_pd(&hood.template src<6, 1, -1>(x));
1659
1660            temp1 = _mm_mul_pd(neig0, _mm_load_pd(&hood.template coeff<0, 0, 0, 8>(x)));
1661            temp2 = _mm_mul_pd(neig1, _mm_load_pd(&hood.template coeff<2, 0, 0, 8>(x)));
1662            temp3 = _mm_mul_pd(neig2, _mm_load_pd(&hood.template coeff<4, 0, 0, 8>(x)));
1663            temp4 = _mm_mul_pd(neig3, _mm_load_pd(&hood.template coeff<6, 0, 0, 8>(x)));
1664
1665            same0 = _mm_add_pd(same0, temp1);
1666            same1 = _mm_add_pd(same1, temp2);
1667            same2 = _mm_add_pd(same2, temp3);
1668            same3 = _mm_add_pd(same3, temp4);
1669
1670            //xxxxxxxxxxxxx
1671            neig0 = _mm_load_pd(&hood.template src<0, -1, 1>(x));
1672            neig1 = _mm_load_pd(&hood.template src<2, -1, 1>(x));
1673            neig2 = _mm_load_pd(&hood.template src<4, -1, 1>(x));
1674            neig3 = _mm_load_pd(&hood.template src<6, -1, 1>(x));
1675
1676            temp1 = _mm_mul_pd(neig0, _mm_load_pd(&hood.template coeff<0, 0, 0, 9>(x)));
1677            temp2 = _mm_mul_pd(neig1, _mm_load_pd(&hood.template coeff<2, 0, 0, 9>(x)));
1678            temp3 = _mm_mul_pd(neig2, _mm_load_pd(&hood.template coeff<4, 0, 0, 9>(x)));
1679            temp4 = _mm_mul_pd(neig3, _mm_load_pd(&hood.template coeff<6, 0, 0, 9>(x)));
1680
1681            same0 = _mm_add_pd(same0, temp1);
1682            same1 = _mm_add_pd(same1, temp2);
1683            same2 = _mm_add_pd(same2, temp3);
1684            same3 = _mm_add_pd(same3, temp4);
1685
1686            //xxxxxxxxxxxxx
1687            neig0 = _mm_load_pd(&hood.template src<0, 1, 1>(x));
1688            neig1 = _mm_load_pd(&hood.template src<2, 1, 1>(x));
1689            neig2 = _mm_load_pd(&hood.template src<4, 1, 1>(x));
1690            neig3 = _mm_load_pd(&hood.template src<6, 1, 1>(x));
1691
1692            temp1 = _mm_mul_pd(neig0, _mm_load_pd(&hood.template coeff<0, 0, 0, 10>(x)));
1693            temp2 = _mm_mul_pd(neig1, _mm_load_pd(&hood.template coeff<2, 0, 0, 10>(x)));
1694            temp3 = _mm_mul_pd(neig2, _mm_load_pd(&hood.template coeff<4, 0, 0, 10>(x)));
1695            temp4 = _mm_mul_pd(neig3, _mm_load_pd(&hood.template coeff<6, 0, 0, 10>(x)));
1696
1697            same0 = _mm_add_pd(same0, temp1);
1698            same1 = _mm_add_pd(same1, temp2);
1699            same2 = _mm_add_pd(same2, temp3);
1700            same3 = _mm_add_pd(same3, temp4);
1701
1702            //xxxxxxxxxxxxx
1703            neig0 = _mm_load_pd(&hood.template coeff<0, -1, -1, 11>(x));
1704            neig1 = _mm_load_pd(&hood.template coeff<2, -1, -1, 11>(x));
1705            neig2 = _mm_load_pd(&hood.template coeff<4, -1, -1, 11>(x));
1706            neig3 = _mm_load_pd(&hood.template coeff<6, -1, -1, 11>(x));
1707
1708            same0 = _mm_add_pd(same0, neig0);
1709            same1 = _mm_add_pd(same1, neig1);
1710            same2 = _mm_add_pd(same2, neig2);
1711            same3 = _mm_add_pd(same3, neig3);
1712
1713            //xxxxxxxxxxxxx
1714            neig0 = _mm_load_pd(&hood.template coeff<0, 0, -1, 11>(x));
1715            neig1 = _mm_load_pd(&hood.template coeff<2, 0, -1, 11>(x));
1716            neig2 = _mm_load_pd(&hood.template coeff<4, 0, -1, 11>(x));
1717            neig3 = _mm_load_pd(&hood.template coeff<6, 0, -1, 11>(x));
1718
1719            same0 = _mm_add_pd(same0, neig0);
1720            same1 = _mm_add_pd(same1, neig1);
1721            same2 = _mm_add_pd(same2, neig2);
1722            same3 = _mm_add_pd(same3, neig3);
1723
1724            //xxxxxxxxxxxxx
1725            neig0 = _mm_load_pd(&hood.template coeff<0, 1, -1, 11>(x));
1726            neig1 = _mm_load_pd(&hood.template coeff<2, 1, -1, 11>(x));
1727            neig2 = _mm_load_pd(&hood.template coeff<4, 1, -1, 11>(x));
1728            neig3 = _mm_load_pd(&hood.template coeff<6, 1, -1, 11>(x));
1729
1730            same0 = _mm_add_pd(same0, neig0);
1731            same1 = _mm_add_pd(same1, neig1);
1732            same2 = _mm_add_pd(same2, neig2);
1733            same3 = _mm_add_pd(same3, neig3);
1734
1735            //xxxxxxxxxxxxx
1736            neig0 = _mm_load_pd(&hood.template coeff<0, -1, 0, 11>(x));
1737            neig1 = _mm_load_pd(&hood.template coeff<2, -1, 0, 11>(x));
1738            neig2 = _mm_load_pd(&hood.template coeff<4, -1, 0, 11>(x));
1739            neig3 = _mm_load_pd(&hood.template coeff<6, -1, 0, 11>(x));
1740
1741            same0 = _mm_add_pd(same0, neig0);
1742            same1 = _mm_add_pd(same1, neig1);
1743            same2 = _mm_add_pd(same2, neig2);
1744            same3 = _mm_add_pd(same3, neig3);
1745
1746            //xxxxxxxxxxxxx
1747            neig0 = _mm_load_pd(&hood.template coeff<0, 1, 0, 11>(x));
1748            neig1 = _mm_load_pd(&hood.template coeff<2, 1, 0, 11>(x));
1749            neig2 = _mm_load_pd(&hood.template coeff<4, 1, 0, 11>(x));
1750            neig3 = _mm_load_pd(&hood.template coeff<6, 1, 0, 11>(x));
1751
1752            same0 = _mm_add_pd(same0, neig0);
1753            same1 = _mm_add_pd(same1, neig1);
1754            same2 = _mm_add_pd(same2, neig2);
1755            same3 = _mm_add_pd(same3, neig3);
1756
1757            //xxxxxxxxxxxxx
1758            neig0 = _mm_load_pd(&hood.template coeff<0, 0, 0, 11>(x));
1759            neig1 = _mm_load_pd(&hood.template coeff<2, 0, 0, 11>(x));
1760            neig2 = _mm_load_pd(&hood.template coeff<4, 0, 0, 11>(x));
1761            neig3 = _mm_load_pd(&hood.template coeff<6, 0, 0, 11>(x));
1762
1763            same0 = _mm_add_pd(same0, neig0);
1764            same1 = _mm_add_pd(same1, neig1);
1765            same2 = _mm_add_pd(same2, neig2);
1766            same3 = _mm_add_pd(same3, neig3);
1767
1768            //xxxxxxxxxxxxx
1769            neig0 = _mm_load_pd(&hood.template coeff<0, -1, 1, 11>(x));
1770            neig1 = _mm_load_pd(&hood.template coeff<2, -1, 1, 11>(x));
1771            neig2 = _mm_load_pd(&hood.template coeff<4, -1, 1, 11>(x));
1772            neig3 = _mm_load_pd(&hood.template coeff<6, -1, 1, 11>(x));
1773
1774            same0 = _mm_add_pd(same0, neig0);
1775            same1 = _mm_add_pd(same1, neig1);
1776            same2 = _mm_add_pd(same2, neig2);
1777            same3 = _mm_add_pd(same3, neig3);
1778
1779            //xxxxxxxxxxxxx
1780            neig0 = _mm_load_pd(&hood.template coeff<0, 0, 1, 11>(x));
1781            neig1 = _mm_load_pd(&hood.template coeff<2, 0, 1, 11>(x));
1782            neig2 = _mm_load_pd(&hood.template coeff<4, 0, 1, 11>(x));
1783            neig3 = _mm_load_pd(&hood.template coeff<6, 0, 1, 11>(x));
1784
1785            same0 = _mm_add_pd(same0, neig0);
1786            same1 = _mm_add_pd(same1, neig1);
1787            same2 = _mm_add_pd(same2, neig2);
1788            same3 = _mm_add_pd(same3, neig3);
1789
1790            //xxxxxxxxxxxxx
1791            neig0 = _mm_load_pd(&hood.template coeff<0, 1, 1, 11>(x));
1792            neig1 = _mm_load_pd(&hood.template coeff<2, 1, 1, 11>(x));
1793            neig2 = _mm_load_pd(&hood.template coeff<4, 1, 1, 11>(x));
1794            neig3 = _mm_load_pd(&hood.template coeff<6, 1, 1, 11>(x));
1795
1796            same0 = _mm_add_pd(same0, neig0);
1797            same1 = _mm_add_pd(same1, neig1);
1798            same2 = _mm_add_pd(same2, neig2);
1799            same3 = _mm_add_pd(same3, neig3);
1800
1801            //xxxxxxxxxxxxx
1802            neig0 = _mm_load_pd(&hood.template coeff<0, -1, -1, 12>(x));
1803            neig1 = _mm_load_pd(&hood.template coeff<2, -1, -1, 12>(x));
1804            neig2 = _mm_load_pd(&hood.template coeff<4, -1, -1, 12>(x));
1805            neig3 = _mm_load_pd(&hood.template coeff<6, -1, -1, 12>(x));
1806
1807            same0 = _mm_add_pd(same0, neig0);
1808            same1 = _mm_add_pd(same1, neig1);
1809            same2 = _mm_add_pd(same2, neig2);
1810            same3 = _mm_add_pd(same3, neig3);
1811
1812            //xxxxxxxxxxxxx
1813            neig0 = _mm_load_pd(&hood.template coeff<0, 0, -1, 12>(x));
1814            neig1 = _mm_load_pd(&hood.template coeff<2, 0, -1, 12>(x));
1815            neig2 = _mm_load_pd(&hood.template coeff<4, 0, -1, 12>(x));
1816            neig3 = _mm_load_pd(&hood.template coeff<6, 0, -1, 12>(x));
1817
1818            same0 = _mm_add_pd(same0, neig0);
1819            same1 = _mm_add_pd(same1, neig1);
1820            same2 = _mm_add_pd(same2, neig2);
1821            same3 = _mm_add_pd(same3, neig3);
1822
1823            //xxxxxxxxxxxxx
1824            neig0 = _mm_load_pd(&hood.template coeff<0, 1, -1, 12>(x));
1825            neig1 = _mm_load_pd(&hood.template coeff<2, 1, -1, 12>(x));
1826            neig2 = _mm_load_pd(&hood.template coeff<4, 1, -1, 12>(x));
1827            neig3 = _mm_load_pd(&hood.template coeff<6, 1, -1, 12>(x));
1828
1829            same0 = _mm_add_pd(same0, neig0);
1830            same1 = _mm_add_pd(same1, neig1);
1831            same2 = _mm_add_pd(same2, neig2);
1832            same3 = _mm_add_pd(same3, neig3);
1833
1834            //xxxxxxxxxxxxx
1835            neig0 = _mm_load_pd(&hood.template coeff<0, -1, 0, 12>(x));
1836            neig1 = _mm_load_pd(&hood.template coeff<2, -1, 0, 12>(x));
1837            neig2 = _mm_load_pd(&hood.template coeff<4, -1, 0, 12>(x));
1838            neig3 = _mm_load_pd(&hood.template coeff<6, -1, 0, 12>(x));
1839
1840            same0 = _mm_add_pd(same0, neig0);
1841            same1 = _mm_add_pd(same1, neig1);
1842            same2 = _mm_add_pd(same2, neig2);
1843            same3 = _mm_add_pd(same3, neig3);
1844
1845            //xxxxxxxxxxxxx
1846            neig0 = _mm_load_pd(&hood.template coeff<0, 1, 0, 12>(x));
1847            neig1 = _mm_load_pd(&hood.template coeff<2, 1, 0, 12>(x));
1848            neig2 = _mm_load_pd(&hood.template coeff<4, 1, 0, 12>(x));
1849            neig3 = _mm_load_pd(&hood.template coeff<6, 1, 0, 12>(x));
1850
1851            same0 = _mm_add_pd(same0, neig0);
1852            same1 = _mm_add_pd(same1, neig1);
1853            same2 = _mm_add_pd(same2, neig2);
1854            same3 = _mm_add_pd(same3, neig3);
1855
1856            //xxxxxxxxxxxxx
1857            neig0 = _mm_load_pd(&hood.template coeff<0, 0, 0, 12>(x));
1858            neig1 = _mm_load_pd(&hood.template coeff<2, 0, 0, 12>(x));
1859            neig2 = _mm_load_pd(&hood.template coeff<4, 0, 0, 12>(x));
1860            neig3 = _mm_load_pd(&hood.template coeff<6, 0, 0, 12>(x));
1861
1862            same0 = _mm_add_pd(same0, neig0);
1863            same1 = _mm_add_pd(same1, neig1);
1864            same2 = _mm_add_pd(same2, neig2);
1865            same3 = _mm_add_pd(same3, neig3);
1866
1867            //xxxxxxxxxxxxx
1868            neig0 = _mm_load_pd(&hood.template coeff<0, -1, 1, 12>(x));
1869            neig1 = _mm_load_pd(&hood.template coeff<2, -1, 1, 12>(x));
1870            neig2 = _mm_load_pd(&hood.template coeff<4, -1, 1, 12>(x));
1871            neig3 = _mm_load_pd(&hood.template coeff<6, -1, 1, 12>(x));
1872
1873            same0 = _mm_add_pd(same0, neig0);
1874            same1 = _mm_add_pd(same1, neig1);
1875            same2 = _mm_add_pd(same2, neig2);
1876            same3 = _mm_add_pd(same3, neig3);
1877
1878            //xxxxxxxxxxxxx
1879            neig0 = _mm_load_pd(&hood.template coeff<0, 0, 1, 12>(x));
1880            neig1 = _mm_load_pd(&hood.template coeff<2, 0, 1, 12>(x));
1881            neig2 = _mm_load_pd(&hood.template coeff<4, 0, 1, 12>(x));
1882            neig3 = _mm_load_pd(&hood.template coeff<6, 0, 1, 12>(x));
1883
1884            same0 = _mm_add_pd(same0, neig0);
1885            same1 = _mm_add_pd(same1, neig1);
1886            same2 = _mm_add_pd(same2, neig2);
1887            same3 = _mm_add_pd(same3, neig3);
1888
1889            //xxxxxxxxxxxxx
1890            neig0 = _mm_load_pd(&hood.template coeff<0, 1, 1, 12>(x));
1891            neig1 = _mm_load_pd(&hood.template coeff<2, 1, 1, 12>(x));
1892            neig2 = _mm_load_pd(&hood.template coeff<4, 1, 1, 12>(x));
1893            neig3 = _mm_load_pd(&hood.template coeff<6, 1, 1, 12>(x));
1894
1895            same0 = _mm_add_pd(same0, neig0);
1896            same1 = _mm_add_pd(same1, neig1);
1897            same2 = _mm_add_pd(same2, neig2);
1898            same3 = _mm_add_pd(same3, neig3);
1899
1900            //yyyyyyyyyyyyy
1901            _mm_store_pd(dst + x + 0, same0);
1902            _mm_store_pd(dst + x + 2, same1);
1903            _mm_store_pd(dst + x + 4, same2);
1904            _mm_store_pd(dst + x + 6, same3);
1905
1906            same0 = same4;
1907            neig0 = neig4;
1908
1909            // dst[x] =
1910            //     coeff[0][x] * src[x - offsetZ] +
1911            //     coeff[1][x] * src[x - offsetY] +
1912            //     coeff[2][x] * src[x - 1] +
1913            //     coeff[3][x] * src[x] +
1914            //     coeff[4][x] * src[x + 1] +
1915            //     coeff[5][x] * src[x + offsetY] +
1916            //     coeff[6][x] * src[x + offsetZ];
1917        }
1918
1919        stepScalar(hood, dst, offsetY, offsetZ, x, endX);
1920    }
1921
1922    template<class NEIGHBORHOOD>
1923    inline void stepScalar(const NEIGHBORHOOD& hood, double *dst, int offsetY, int offsetZ, int startX, int endX)
1924    {
1925        for (int x = startX; x < endX; ++x) {
1926            dst[x] = 
1927                hood.template coeff<0, 0, 0, 0>(x) * hood.template src<0,    -1, -1>(x) +
1928                hood.template coeff<0, 0, 0, 1>(x) * hood.template src<0,     0, -1>(x) +
1929                hood.template coeff<0, 0, 0, 2>(x) * hood.template src<0,     1, -1>(x) +
1930                hood.template coeff<0, 0, 0, 3>(x) * hood.template src<0,    -1,  0>(x) + 
1931                hood.template coeff<0, 0, 0, 4>(x) * hood.template src<0 - 1, 0,  0>(x) +
1932                hood.template coeff<0, 0, 0, 5>(x) * hood.template src<0,     0,  0>(x) +
1933                hood.template coeff<0, 0, 0, 6>(x) * hood.template src<0 + 1, 0,  0>(x) +
1934                hood.template coeff<0, 0, 0, 7>(x) * hood.template src<0,     1,  0>(x) +
1935                hood.template coeff<0, 0, 0, 8>(x) * hood.template src<0,    -1,  1>(x) +
1936                hood.template coeff<0, 0, 0, 9>(x) * hood.template src<0,     0,  1>(x) +
1937                hood.template coeff<0, 0, 0,10>(x) * hood.template src<0,     1,  1>(x) +
1938
1939                hood.template coeff<0, -1, -1, 11>(x) +
1940                hood.template coeff<0,  0, -1, 11>(x) +
1941                hood.template coeff<0,  1, -1, 11>(x) +
1942                hood.template coeff<0, -1,  0, 11>(x) +
1943                hood.template coeff<-1, 0,  0, 11>(x) +
1944                hood.template coeff<0,  0,  0, 11>(x) +
1945                hood.template coeff<1,  0,  0, 11>(x) +
1946                hood.template coeff<0,  1,  0, 11>(x) +
1947                hood.template coeff<0, -1,  1, 11>(x) +
1948                hood.template coeff<0,  0,  1, 11>(x) +
1949                hood.template coeff<0,  1,  1, 11>(x) +
1950
1951                hood.template coeff<0, -1, -1, 12>(x) +
1952                hood.template coeff<0,  0, -1, 12>(x) +
1953                hood.template coeff<0,  1, -1, 12>(x) +
1954                hood.template coeff<0, -1,  0, 12>(x) +
1955                hood.template coeff<-1, 0,  0, 12>(x) +
1956                hood.template coeff<0,  0,  0, 12>(x) +
1957                hood.template coeff<1,  0,  0, 12>(x) +
1958                hood.template coeff<0,  1,  0, 12>(x) +
1959                hood.template coeff<0, -1,  1, 12>(x) +
1960                hood.template coeff<0,  0,  1, 12>(x) +
1961                hood.template coeff<0,  1,  1, 12>(x);
1962        }
1963    }
1964
1965    int flops()
1966    {
1967        return 40;
1968    }
1969};
1970
1971class ExtendedVectorized3DDefiant
1972{
1973public:
1974    static int coefficients()
1975    {
1976        return 13;
1977    }
1978
1979    template<class NEIGHBORHOOD>
1980    inline void step(const NEIGHBORHOOD& hood, double *dst, int startX, int endX)
1981    {
1982        int x = startX;
1983
1984        if ((x & 1) == 1) {
1985            stepScalar(hood, dst, x, x + 1);
1986            x += 1;
1987        }
1988
1989        __m128d same0 = _mm_load_pd(&(hood[FC<0, 0, 0>()].srcB(x)));
1990        __m128d neig0 = _mm_loadu_pd(&(hood[FC<1, 0, 0>()].srcB(x)));
1991
1992        int paddedEndX = endX - 7;
1993        for (; x < paddedEndX; x += 8) {
1994            __m128d same1 = _mm_load_pd(&hood[FC<2, 0, 0>()].srcB(x));
1995            __m128d same2 = _mm_load_pd(&hood[FC<4, 0, 0>()].srcB(x));
1996            __m128d same3 = _mm_load_pd(&hood[FC<6, 0, 0>()].srcB(x));
1997            __m128d same4 = _mm_load_pd(&hood[FC<8, 0, 0>()].srcB(x));
1998
1999            __m128d neig1 = _mm_shuffle_pd(same0, same1, (1 << 0) | (0 << 2));
2000            __m128d neig2 = _mm_shuffle_pd(same1, same2, (1 << 0) | (0 << 2));
2001            __m128d neig3 = _mm_shuffle_pd(same2, same3, (1 << 0) | (0 << 2));
2002            __m128d neig4 = _mm_shuffle_pd(same3, same4, (1 << 0) | (0 << 2));
2003
2004            same0 = _mm_mul_pd(same0, _mm_load_pd(&hood[FC<0, 0, 0>()].coeffB(FC< 3>(), x)));
2005            same1 = _mm_mul_pd(same1, _mm_load_pd(&hood[FC<2, 0, 0>()].coeffB(FC< 3>(), x)));
2006            same2 = _mm_mul_pd(same2, _mm_load_pd(&hood[FC<4, 0, 0>()].coeffB(FC< 3>(), x)));
2007            same3 = _mm_mul_pd(same3, _mm_load_pd(&hood[FC<6, 0, 0>()].coeffB(FC< 3>(), x)));
2008
2009            __m128d temp1 = _mm_mul_pd(neig0, _mm_load_pd(&hood[FC<0, 0, 0>()].coeffB(FC< 2>(), x)));
2010            __m128d temp2 = _mm_mul_pd(neig1, _mm_load_pd(&hood[FC<2, 0, 0>()].coeffB(FC< 2>(), x)));
2011            __m128d temp3 = _mm_mul_pd(neig2, _mm_load_pd(&hood[FC<4, 0, 0>()].coeffB(FC< 2>(), x)));
2012            __m128d temp4 = _mm_mul_pd(neig3, _mm_load_pd(&hood[FC<6, 0, 0>()].coeffB(FC< 2>(), x)));
2013
2014            same0 = _mm_add_pd(same0, temp1);
2015            same1 = _mm_add_pd(same1, temp2);
2016            same2 = _mm_add_pd(same2, temp3);
2017            same3 = _mm_add_pd(same3, temp4);
2018
2019            temp1 = _mm_mul_pd(neig1, _mm_load_pd(&hood[FC<0, 0, 0>()].coeffB(FC< 4>(), x)));
2020            temp2 = _mm_mul_pd(neig2, _mm_load_pd(&hood[FC<2, 0, 0>()].coeffB(FC< 4>(), x)));
2021            temp3 = _mm_mul_pd(neig3, _mm_load_pd(&hood[FC<4, 0, 0>()].coeffB(FC< 4>(), x)));
2022            temp4 = _mm_mul_pd(neig4, _mm_load_pd(&hood[FC<6, 0, 0>()].coeffB(FC< 4>(), x)));
2023
2024            same0 = _mm_add_pd(same0, temp1);
2025            same1 = _mm_add_pd(same1, temp2);
2026            same2 = _mm_add_pd(same2, temp3);
2027            same3 = _mm_add_pd(same3, temp4);
2028
2029            neig0 = _mm_load_pd(&hood[FC<0, 0, -1>()].srcB(x));
2030            neig1 = _mm_load_pd(&hood[FC<2, 0, -1>()].srcB(x));
2031            neig2 = _mm_load_pd(&hood[FC<4, 0, -1>()].srcB(x));
2032            neig3 = _mm_load_pd(&hood[FC<6, 0, -1>()].srcB(x));
2033
2034            temp1 = _mm_mul_pd(neig0, _mm_load_pd(&hood[FC<0, 0, 0>()].coeffB(FC< 0>(), x)));
2035            temp2 = _mm_mul_pd(neig1, _mm_load_pd(&hood[FC<2, 0, 0>()].coeffB(FC< 0>(), x)));
2036            temp3 = _mm_mul_pd(neig2, _mm_load_pd(&hood[FC<4, 0, 0>()].coeffB(FC< 0>(), x)));
2037            temp4 = _mm_mul_pd(neig3, _mm_load_pd(&hood[FC<6, 0, 0>()].coeffB(FC< 0>(), x)));
2038
2039            same0 = _mm_add_pd(same0, temp1);
2040            same1 = _mm_add_pd(same1, temp2);
2041            same2 = _mm_add_pd(same2, temp3);
2042            same3 = _mm_add_pd(same3, temp4);
2043
2044            neig0 = _mm_load_pd(&hood[FC<0, -1, 0>()].srcB(x));
2045            neig1 = _mm_load_pd(&hood[FC<2, -1, 0>()].srcB(x));
2046            neig2 = _mm_load_pd(&hood[FC<4, -1, 0>()].srcB(x));
2047            neig3 = _mm_load_pd(&hood[FC<6, -1, 0>()].srcB(x));
2048
2049            temp1 = _mm_mul_pd(neig0, _mm_load_pd(&hood[FC<0, 0, 0>()].coeffB(FC< 1>(), x)));
2050            temp2 = _mm_mul_pd(neig1, _mm_load_pd(&hood[FC<2, 0, 0>()].coeffB(FC< 1>(), x)));
2051            temp3 = _mm_mul_pd(neig2, _mm_load_pd(&hood[FC<4, 0, 0>()].coeffB(FC< 1>(), x)));
2052            temp4 = _mm_mul_pd(neig3, _mm_load_pd(&hood[FC<6, 0, 0>()].coeffB(FC< 1>(), x)));
2053
2054            same0 = _mm_add_pd(same0, temp1);
2055            same1 = _mm_add_pd(same1, temp2);
2056            same2 = _mm_add_pd(same2, temp3);
2057            same3 = _mm_add_pd(same3, temp4);
2058
2059            neig0 = _mm_load_pd(&hood[FC<0, 1, 0>()].srcB(x));
2060            neig1 = _mm_load_pd(&hood[FC<2, 1, 0>()].srcB(x));
2061            neig2 = _mm_load_pd(&hood[FC<4, 1, 0>()].srcB(x));
2062            neig3 = _mm_load_pd(&hood[FC<6, 1, 0>()].srcB(x));
2063
2064            temp1 = _mm_mul_pd(neig0, _mm_load_pd(&hood[FC<0, 0, 0>()].coeffB(FC< 5>(), x)));
2065            temp2 = _mm_mul_pd(neig1, _mm_load_pd(&hood[FC<2, 0, 0>()].coeffB(FC< 5>(), x)));
2066            temp3 = _mm_mul_pd(neig2, _mm_load_pd(&hood[FC<4, 0, 0>()].coeffB(FC< 5>(), x)));
2067            temp4 = _mm_mul_pd(neig3, _mm_load_pd(&hood[FC<6, 0, 0>()].coeffB(FC< 5>(), x)));
2068
2069            same0 = _mm_add_pd(same0, temp1);
2070            same1 = _mm_add_pd(same1, temp2);
2071            same2 = _mm_add_pd(same2, temp3);
2072            same3 = _mm_add_pd(same3, temp4);
2073
2074            //xxxxxxxxxxxxx
2075            neig0 = _mm_load_pd(&hood[FC<0, 0, 1>()].srcB(x));
2076            neig1 = _mm_load_pd(&hood[FC<2, 0, 1>()].srcB(x));
2077            neig2 = _mm_load_pd(&hood[FC<4, 0, 1>()].srcB(x));
2078            neig3 = _mm_load_pd(&hood[FC<6, 0, 1>()].srcB(x));
2079
2080            temp1 = _mm_mul_pd(neig0, _mm_load_pd(&hood[FC<0, 0, 0>()].coeffB(FC< 6>(), x)));
2081            temp2 = _mm_mul_pd(neig1, _mm_load_pd(&hood[FC<2, 0, 0>()].coeffB(FC< 6>(), x)));
2082            temp3 = _mm_mul_pd(neig2, _mm_load_pd(&hood[FC<4, 0, 0>()].coeffB(FC< 6>(), x)));
2083            temp4 = _mm_mul_pd(neig3, _mm_load_pd(&hood[FC<6, 0, 0>()].coeffB(FC< 6>(), x)));
2084
2085            same0 = _mm_add_pd(same0, temp1);
2086            same1 = _mm_add_pd(same1, temp2);
2087            same2 = _mm_add_pd(same2, temp3);
2088            same3 = _mm_add_pd(same3, temp4);
2089
2090            //xxxxxxxxxxxxx
2091            neig0 = _mm_load_pd(&hood[FC<0, -1, -1>()].srcB(x));
2092            neig1 = _mm_load_pd(&hood[FC<2, -1, -1>()].srcB(x));
2093            neig2 = _mm_load_pd(&hood[FC<4, -1, -1>()].srcB(x));
2094            neig3 = _mm_load_pd(&hood[FC<6, -1, -1>()].srcB(x));
2095
2096            temp1 = _mm_mul_pd(neig0, _mm_load_pd(&hood[FC<0, 0, 0>()].coeffB(FC< 7>(), x)));
2097            temp2 = _mm_mul_pd(neig1, _mm_load_pd(&hood[FC<2, 0, 0>()].coeffB(FC< 7>(), x)));
2098            temp3 = _mm_mul_pd(neig2, _mm_load_pd(&hood[FC<4, 0, 0>()].coeffB(FC< 7>(), x)));
2099            temp4 = _mm_mul_pd(neig3, _mm_load_pd(&hood[FC<6, 0, 0>()].coeffB(FC< 7>(), x)));
2100
2101            same0 = _mm_add_pd(same0, temp1);
2102            same1 = _mm_add_pd(same1, temp2);
2103            same2 = _mm_add_pd(same2, temp3);
2104            same3 = _mm_add_pd(same3, temp4);
2105
2106            //xxxxxxxxxxxxx
2107            neig0 = _mm_load_pd(&hood[FC<0, 1, -1>()].srcB(x));
2108            neig1 = _mm_load_pd(&hood[FC<2, 1, -1>()].srcB(x));
2109            neig2 = _mm_load_pd(&hood[FC<4, 1, -1>()].srcB(x));
2110            neig3 = _mm_load_pd(&hood[FC<6, 1, -1>()].srcB(x));
2111
2112            temp1 = _mm_mul_pd(neig0, _mm_load_pd(&hood[FC<0, 0, 0>()].coeffB(FC< 8>(), x)));
2113            temp2 = _mm_mul_pd(neig1, _mm_load_pd(&hood[FC<2, 0, 0>()].coeffB(FC< 8>(), x)));
2114            temp3 = _mm_mul_pd(neig2, _mm_load_pd(&hood[FC<4, 0, 0>()].coeffB(FC< 8>(), x)));
2115            temp4 = _mm_mul_pd(neig3, _mm_load_pd(&hood[FC<6, 0, 0>()].coeffB(FC< 8>(), x)));
2116
2117            same0 = _mm_add_pd(same0, temp1);
2118            same1 = _mm_add_pd(same1, temp2);
2119            same2 = _mm_add_pd(same2, temp3);
2120            same3 = _mm_add_pd(same3, temp4);
2121
2122            //xxxxxxxxxxxxx
2123            neig0 = _mm_load_pd(&hood[FC<0, -1, 1>()].srcB(x));
2124            neig1 = _mm_load_pd(&hood[FC<2, -1, 1>()].srcB(x));
2125            neig2 = _mm_load_pd(&hood[FC<4, -1, 1>()].srcB(x));
2126            neig3 = _mm_load_pd(&hood[FC<6, -1, 1>()].srcB(x));
2127
2128            temp1 = _mm_mul_pd(neig0, _mm_load_pd(&hood[FC<0, 0, 0>()].coeffB(FC< 9>(), x)));
2129            temp2 = _mm_mul_pd(neig1, _mm_load_pd(&hood[FC<2, 0, 0>()].coeffB(FC< 9>(), x)));
2130            temp3 = _mm_mul_pd(neig2, _mm_load_pd(&hood[FC<4, 0, 0>()].coeffB(FC< 9>(), x)));
2131            temp4 = _mm_mul_pd(neig3, _mm_load_pd(&hood[FC<6, 0, 0>()].coeffB(FC< 9>(), x)));
2132
2133            same0 = _mm_add_pd(same0, temp1);
2134            same1 = _mm_add_pd(same1, temp2);
2135            same2 = _mm_add_pd(same2, temp3);
2136            same3 = _mm_add_pd(same3, temp4);
2137
2138            //xxxxxxxxxxxxx
2139            neig0 = _mm_load_pd(&hood[FC<0, 1, 1>()].srcB(x));
2140            neig1 = _mm_load_pd(&hood[FC<2, 1, 1>()].srcB(x));
2141            neig2 = _mm_load_pd(&hood[FC<4, 1, 1>()].srcB(x));
2142            neig3 = _mm_load_pd(&hood[FC<6, 1, 1>()].srcB(x));
2143
2144            temp1 = _mm_mul_pd(neig0, _mm_load_pd(&hood[FC<0, 0, 0>()].coeffB(FC<10>(), x)));
2145            temp2 = _mm_mul_pd(neig1, _mm_load_pd(&hood[FC<2, 0, 0>()].coeffB(FC<10>(), x)));
2146            temp3 = _mm_mul_pd(neig2, _mm_load_pd(&hood[FC<4, 0, 0>()].coeffB(FC<10>(), x)));
2147            temp4 = _mm_mul_pd(neig3, _mm_load_pd(&hood[FC<6, 0, 0>()].coeffB(FC<10>(), x)));
2148
2149            same0 = _mm_add_pd(same0, temp1);
2150            same1 = _mm_add_pd(same1, temp2);
2151            same2 = _mm_add_pd(same2, temp3);
2152            same3 = _mm_add_pd(same3, temp4);
2153
2154            //xxxxxxxxxxxxx
2155            neig0 = _mm_load_pd(&hood[FC<0, -1, -1>()].coeffB(FC<11>(), x));
2156            neig1 = _mm_load_pd(&hood[FC<2, -1, -1>()].coeffB(FC<11>(), x));
2157            neig2 = _mm_load_pd(&hood[FC<4, -1, -1>()].coeffB(FC<11>(), x));
2158            neig3 = _mm_load_pd(&hood[FC<6, -1, -1>()].coeffB(FC<11>(), x));
2159
2160            same0 = _mm_add_pd(same0, neig0);
2161            same1 = _mm_add_pd(same1, neig1);
2162            same2 = _mm_add_pd(same2, neig2);
2163            same3 = _mm_add_pd(same3, neig3);
2164
2165            //xxxxxxxxxxxxx
2166            neig0 = _mm_load_pd(&hood[FC<0, 0, -1>()].coeffB(FC<11>(), x));
2167            neig1 = _mm_load_pd(&hood[FC<2, 0, -1>()].coeffB(FC<11>(), x));
2168            neig2 = _mm_load_pd(&hood[FC<4, 0, -1>()].coeffB(FC<11>(), x));
2169            neig3 = _mm_load_pd(&hood[FC<6, 0, -1>()].coeffB(FC<11>(), x));
2170
2171            same0 = _mm_add_pd(same0, neig0);
2172            same1 = _mm_add_pd(same1, neig1);
2173            same2 = _mm_add_pd(same2, neig2);
2174            same3 = _mm_add_pd(same3, neig3);
2175
2176            //xxxxxxxxxxxxx
2177            neig0 = _mm_load_pd(&hood[FC<0, 1, -1>()].coeffB(FC<11>(), x));
2178            neig1 = _mm_load_pd(&hood[FC<2, 1, -1>()].coeffB(FC<11>(), x));
2179            neig2 = _mm_load_pd(&hood[FC<4, 1, -1>()].coeffB(FC<11>(), x));
2180            neig3 = _mm_load_pd(&hood[FC<6, 1, -1>()].coeffB(FC<11>(), x));
2181
2182            same0 = _mm_add_pd(same0, neig0);
2183            same1 = _mm_add_pd(same1, neig1);
2184            same2 = _mm_add_pd(same2, neig2);
2185            same3 = _mm_add_pd(same3, neig3);
2186
2187            //xxxxxxxxxxxxx
2188            neig0 = _mm_load_pd(&hood[FC<0, -1, 0>()].coeffB(FC<11>(), x));
2189            neig1 = _mm_load_pd(&hood[FC<2, -1, 0>()].coeffB(FC<11>(), x));
2190            neig2 = _mm_load_pd(&hood[FC<4, -1, 0>()].coeffB(FC<11>(), x));
2191            neig3 = _mm_load_pd(&hood[FC<6, -1, 0>()].coeffB(FC<11>(), x));
2192
2193            same0 = _mm_add_pd(same0, neig0);
2194            same1 = _mm_add_pd(same1, neig1);
2195            same2 = _mm_add_pd(same2, neig2);
2196            same3 = _mm_add_pd(same3, neig3);
2197
2198            //xxxxxxxxxxxxx
2199            neig0 = _mm_load_pd(&hood[FC<0, 1, 0>()].coeffB(FC<11>(), x));
2200            neig1 = _mm_load_pd(&hood[FC<2, 1, 0>()].coeffB(FC<11>(), x));
2201            neig2 = _mm_load_pd(&hood[FC<4, 1, 0>()].coeffB(FC<11>(), x));
2202            neig3 = _mm_load_pd(&hood[FC<6, 1, 0>()].coeffB(FC<11>(), x));
2203
2204            same0 = _mm_add_pd(same0, neig0);
2205            same1 = _mm_add_pd(same1, neig1);
2206            same2 = _mm_add_pd(same2, neig2);
2207            same3 = _mm_add_pd(same3, neig3);
2208
2209            //xxxxxxxxxxxxx
2210            neig0 = _mm_load_pd(&hood[FC<0, 0, 0>()].coeffB(FC<11>(), x));
2211            neig1 = _mm_load_pd(&hood[FC<2, 0, 0>()].coeffB(FC<11>(), x));
2212            neig2 = _mm_load_pd(&hood[FC<4, 0, 0>()].coeffB(FC<11>(), x));
2213            neig3 = _mm_load_pd(&hood[FC<6, 0, 0>()].coeffB(FC<11>(), x));
2214
2215            same0 = _mm_add_pd(same0, neig0);
2216            same1 = _mm_add_pd(same1, neig1);
2217            same2 = _mm_add_pd(same2, neig2);
2218            same3 = _mm_add_pd(same3, neig3);
2219
2220            //xxxxxxxxxxxxx
2221            neig0 = _mm_load_pd(&hood[FC<0, -1, 1>()].coeffB(FC<11>(), x));
2222            neig1 = _mm_load_pd(&hood[FC<2, -1, 1>()].coeffB(FC<11>(), x));
2223            neig2 = _mm_load_pd(&hood[FC<4, -1, 1>()].coeffB(FC<11>(), x));
2224            neig3 = _mm_load_pd(&hood[FC<6, -1, 1>()].coeffB(FC<11>(), x));
2225
2226            same0 = _mm_add_pd(same0, neig0);
2227            same1 = _mm_add_pd(same1, neig1);
2228            same2 = _mm_add_pd(same2, neig2);
2229            same3 = _mm_add_pd(same3, neig3);
2230
2231            //xxxxxxxxxxxxx
2232            neig0 = _mm_load_pd(&hood[FC<0, 0, 1>()].coeffB(FC<11>(), x));
2233            neig1 = _mm_load_pd(&hood[FC<2, 0, 1>()].coeffB(FC<11>(), x));
2234            neig2 = _mm_load_pd(&hood[FC<4, 0, 1>()].coeffB(FC<11>(), x));
2235            neig3 = _mm_load_pd(&hood[FC<6, 0, 1>()].coeffB(FC<11>(), x));
2236
2237            same0 = _mm_add_pd(same0, neig0);
2238            same1 = _mm_add_pd(same1, neig1);
2239            same2 = _mm_add_pd(same2, neig2);
2240            same3 = _mm_add_pd(same3, neig3);
2241
2242            //xxxxxxxxxxxxx
2243            neig0 = _mm_load_pd(&hood[FC<0, 1, 1>()].coeffB(FC<11>(), x));
2244            neig1 = _mm_load_pd(&hood[FC<2, 1, 1>()].coeffB(FC<11>(), x));
2245            neig2 = _mm_load_pd(&hood[FC<4, 1, 1>()].coeffB(FC<11>(), x));
2246            neig3 = _mm_load_pd(&hood[FC<6, 1, 1>()].coeffB(FC<11>(), x));
2247
2248            same0 = _mm_add_pd(same0, neig0);
2249            same1 = _mm_add_pd(same1, neig1);
2250            same2 = _mm_add_pd(same2, neig2);
2251            same3 = _mm_add_pd(same3, neig3);
2252
2253            //xxxxxxxxxxxxx
2254            neig0 = _mm_load_pd(&hood[FC<0, -1, -1>()].coeffB(FC<12>(), x));
2255            neig1 = _mm_load_pd(&hood[FC<2, -1, -1>()].coeffB(FC<12>(), x));
2256            neig2 = _mm_load_pd(&hood[FC<4, -1, -1>()].coeffB(FC<12>(), x));
2257            neig3 = _mm_load_pd(&hood[FC<6, -1, -1>()].coeffB(FC<12>(), x));
2258
2259            same0 = _mm_add_pd(same0, neig0);
2260            same1 = _mm_add_pd(same1, neig1);
2261            same2 = _mm_add_pd(same2, neig2);
2262            same3 = _mm_add_pd(same3, neig3);
2263
2264            //xxxxxxxxxxxxx
2265            neig0 = _mm_load_pd(&hood[FC<0, 0, -1>()].coeffB(FC<12>(), x));
2266            neig1 = _mm_load_pd(&hood[FC<2, 0, -1>()].coeffB(FC<12>(), x));
2267            neig2 = _mm_load_pd(&hood[FC<4, 0, -1>()].coeffB(FC<12>(), x));
2268            neig3 = _mm_load_pd(&hood[FC<6, 0, -1>()].coeffB(FC<12>(), x));
2269
2270            same0 = _mm_add_pd(same0, neig0);
2271            same1 = _mm_add_pd(same1, neig1);
2272            same2 = _mm_add_pd(same2, neig2);
2273            same3 = _mm_add_pd(same3, neig3);
2274
2275            //xxxxxxxxxxxxx
2276            neig0 = _mm_load_pd(&hood[FC<0, 1, -1>()].coeffB(FC<12>(), x));
2277            neig1 = _mm_load_pd(&hood[FC<2, 1, -1>()].coeffB(FC<12>(), x));
2278            neig2 = _mm_load_pd(&hood[FC<4, 1, -1>()].coeffB(FC<12>(), x));
2279            neig3 = _mm_load_pd(&hood[FC<6, 1, -1>()].coeffB(FC<12>(), x));
2280
2281            same0 = _mm_add_pd(same0, neig0);
2282            same1 = _mm_add_pd(same1, neig1);
2283            same2 = _mm_add_pd(same2, neig2);
2284            same3 = _mm_add_pd(same3, neig3);
2285
2286            //xxxxxxxxxxxxx
2287            neig0 = _mm_load_pd(&hood[FC<0, -1, 0>()].coeffB(FC<12>(), x));
2288            neig1 = _mm_load_pd(&hood[FC<2, -1, 0>()].coeffB(FC<12>(), x));
2289            neig2 = _mm_load_pd(&hood[FC<4, -1, 0>()].coeffB(FC<12>(), x));
2290            neig3 = _mm_load_pd(&hood[FC<6, -1, 0>()].coeffB(FC<12>(), x));
2291
2292            same0 = _mm_add_pd(same0, neig0);
2293            same1 = _mm_add_pd(same1, neig1);
2294            same2 = _mm_add_pd(same2, neig2);
2295            same3 = _mm_add_pd(same3, neig3);
2296
2297            //xxxxxxxxxxxxx
2298            neig0 = _mm_load_pd(&hood[FC<0, 1, 0>()].coeffB(FC<12>(), x));
2299            neig1 = _mm_load_pd(&hood[FC<2, 1, 0>()].coeffB(FC<12>(), x));
2300            neig2 = _mm_load_pd(&hood[FC<4, 1, 0>()].coeffB(FC<12>(), x));
2301            neig3 = _mm_load_pd(&hood[FC<6, 1, 0>()].coeffB(FC<12>(), x));
2302
2303            same0 = _mm_add_pd(same0, neig0);
2304            same1 = _mm_add_pd(same1, neig1);
2305            same2 = _mm_add_pd(same2, neig2);
2306            same3 = _mm_add_pd(same3, neig3);
2307
2308            //xxxxxxxxxxxxx
2309            neig0 = _mm_load_pd(&hood[FC<0, 0, 0>()].coeffB(FC<12>(), x));
2310            neig1 = _mm_load_pd(&hood[FC<2, 0, 0>()].coeffB(FC<12>(), x));
2311            neig2 = _mm_load_pd(&hood[FC<4, 0, 0>()].coeffB(FC<12>(), x));
2312            neig3 = _mm_load_pd(&hood[FC<6, 0, 0>()].coeffB(FC<12>(), x));
2313
2314            same0 = _mm_add_pd(same0, neig0);
2315            same1 = _mm_add_pd(same1, neig1);
2316            same2 = _mm_add_pd(same2, neig2);
2317            same3 = _mm_add_pd(same3, neig3);
2318
2319            //xxxxxxxxxxxxx
2320            neig0 = _mm_load_pd(&hood[FC<0, -1, 1>()].coeffB(FC<12>(), x));
2321            neig1 = _mm_load_pd(&hood[FC<2, -1, 1>()].coeffB(FC<12>(), x));
2322            neig2 = _mm_load_pd(&hood[FC<4, -1, 1>()].coeffB(FC<12>(), x));
2323            neig3 = _mm_load_pd(&hood[FC<6, -1, 1>()].coeffB(FC<12>(), x));
2324
2325            same0 = _mm_add_pd(same0, neig0);
2326            same1 = _mm_add_pd(same1, neig1);
2327            same2 = _mm_add_pd(same2, neig2);
2328            same3 = _mm_add_pd(same3, neig3);
2329
2330            //xxxxxxxxxxxxx
2331            neig0 = _mm_load_pd(&hood[FC<0, 0, 1>()].coeffB(FC<12>(), x));
2332            neig1 = _mm_load_pd(&hood[FC<2, 0, 1>()].coeffB(FC<12>(), x));
2333            neig2 = _mm_load_pd(&hood[FC<4, 0, 1>()].coeffB(FC<12>(), x));
2334            neig3 = _mm_load_pd(&hood[FC<6, 0, 1>()].coeffB(FC<12>(), x));
2335
2336            same0 = _mm_add_pd(same0, neig0);
2337            same1 = _mm_add_pd(same1, neig1);
2338            same2 = _mm_add_pd(same2, neig2);
2339            same3 = _mm_add_pd(same3, neig3);
2340
2341            //xxxxxxxxxxxxx
2342            neig0 = _mm_load_pd(&hood[FC<0, 1, 1>()].coeffB(FC<12>(), x));
2343            neig1 = _mm_load_pd(&hood[FC<2, 1, 1>()].coeffB(FC<12>(), x));
2344            neig2 = _mm_load_pd(&hood[FC<4, 1, 1>()].coeffB(FC<12>(), x));
2345            neig3 = _mm_load_pd(&hood[FC<6, 1, 1>()].coeffB(FC<12>(), x));
2346
2347            same0 = _mm_add_pd(same0, neig0);
2348            same1 = _mm_add_pd(same1, neig1);
2349            same2 = _mm_add_pd(same2, neig2);
2350            same3 = _mm_add_pd(same3, neig3);
2351
2352            //yyyyyyyyyyyyy
2353            _mm_store_pd(dst + x + 0, same0);
2354            _mm_store_pd(dst + x + 2, same1);
2355            _mm_store_pd(dst + x + 4, same2);
2356            _mm_store_pd(dst + x + 6, same3);
2357
2358            same0 = same4;
2359            neig0 = neig4;
2360
2361            // dst[x] =
2362            //     coeff[0][x] * src[x - offsetZ] +
2363            //     coeff[1][x] * src[x - offsetY] +
2364            //     coeff[2][x] * src[x - 1] +
2365            //     coeff[3][x] * src[x] +
2366            //     coeff[4][x] * src[x + 1] +
2367            //     coeff[5][x] * src[x + offsetY] +
2368            //     coeff[6][x] * src[x + offsetZ];
2369        }
2370
2371        stepScalar(hood, dst, x, endX);
2372    }
2373
2374    template<class NEIGHBORHOOD>
2375    inline void stepScalar(const NEIGHBORHOOD& hood, double *dst, int startX, int endX)
2376    {
2377        for (int x = startX; x < endX; ++x) {
2378            dst[x] = 
2379                hood[FC<0, 0, 0>()].coeffB(FC< 0>(), x) * hood[FC< 0, -1, -1>()].srcB(x) +
2380                hood[FC<0, 0, 0>()].coeffB(FC< 1>(), x) * hood[FC< 0,  0, -1>()].srcB(x) +
2381                hood[FC<0, 0, 0>()].coeffB(FC< 2>(), x) * hood[FC< 0,  1, -1>()].srcB(x) +
2382                hood[FC<0, 0, 0>()].coeffB(FC< 3>(), x) * hood[FC< 0, -1,  0>()].srcB(x) +
2383                hood[FC<0, 0, 0>()].coeffB(FC< 4>(), x) * hood[FC<-1,  0,  0>()].srcB(x) +
2384                hood[FC<0, 0, 0>()].coeffB(FC< 5>(), x) * hood[FC< 0,  0,  0>()].srcB(x) +
2385                hood[FC<0, 0, 0>()].coeffB(FC< 6>(), x) * hood[FC< 1,  0,  0>()].srcB(x) +
2386                hood[FC<0, 0, 0>()].coeffB(FC< 7>(), x) * hood[FC< 0,  1,  0>()].srcB(x) +
2387                hood[FC<0, 0, 0>()].coeffB(FC< 8>(), x) * hood[FC< 0, -1,  1>()].srcB(x) +
2388                hood[FC<0, 0, 0>()].coeffB(FC< 9>(), x) * hood[FC< 0,  0,  1>()].srcB(x) +
2389                hood[FC<0, 0, 0>()].coeffB(FC<10>(), x) * hood[FC< 0,  1,  1>()].srcB(x) +
2390
2391                hood[FC<0, -1, -1>()].coeffB(FC<11>(), x) +
2392                hood[FC<0,  0, -1>()].coeffB(FC<11>(), x) +
2393                hood[FC<0,  1, -1>()].coeffB(FC<11>(), x) +
2394                hood[FC<0, -1,  0>()].coeffB(FC<11>(), x) +
2395                hood[FC<-1, 0,  0>()].coeffB(FC<11>(), x) +
2396                hood[FC<0,  0,  0>()].coeffB(FC<11>(), x) +
2397                hood[FC<1,  0,  0>()].coeffB(FC<11>(), x) +
2398                hood[FC<0,  1,  0>()].coeffB(FC<11>(), x) +
2399                hood[FC<0, -1,  1>()].coeffB(FC<11>(), x) +
2400                hood[FC<0,  0,  1>()].coeffB(FC<11>(), x) +
2401                hood[FC<0,  1,  1>()].coeffB(FC<11>(), x) +
2402
2403                hood[FC<0, -1, -1>()].coeffB(FC<12>(), x) +
2404                hood[FC<0,  0, -1>()].coeffB(FC<12>(), x) +
2405                hood[FC<0,  1, -1>()].coeffB(FC<12>(), x) +
2406                hood[FC<0, -1,  0>()].coeffB(FC<12>(), x) +
2407                hood[FC<-1, 0,  0>()].coeffB(FC<12>(), x) +
2408                hood[FC<0,  0,  0>()].coeffB(FC<12>(), x) +
2409                hood[FC<1,  0,  0>()].coeffB(FC<12>(), x) +
2410                hood[FC<0,  1,  0>()].coeffB(FC<12>(), x) +
2411                hood[FC<0, -1,  1>()].coeffB(FC<12>(), x) +
2412                hood[FC<0,  0,  1>()].coeffB(FC<12>(), x) +
2413                hood[FC<0,  1,  1>()].coeffB(FC<12>(), x);
2414        }
2415    }
2416
2417    int flops()
2418    {
2419        return 40;
2420    }
2421};
2422
2423int main(int argc, char *argv[])
2424{
2425    // Benchmark<Scalar3D, 3>().exercise();
2426    // Benchmark<Vectorized3D, 3>().exercise();
2427    // Benchmark<VectorizedSSEMelbourneShuffle2D>().exercise();
2428    // Benchmark<ExtendedVectorized3D, 3>().exercise();
2429    // Benchmark<ExtendedVectorized3DFixed, 3>().exercise();
2430    // Benchmark<ExtendedVectorized3DNextGen<4, 4, MY_SIZE>, 3>().exercise();
2431    Benchmark<ExtendedVectorized3DDefiant, 3>().exercise();
2432    // Benchmark<Jacobi3D, 3>().exercise();
2433
2434
2435    // std::vector<cl::Platform> platforms;
2436    // cl::Platform::get(&platforms);
2437
2438    // for (int i = 0; i < platforms.size(); ++i) {
2439    //     std::string str;
2440    //     platforms[i].getInfo(CL_PLATFORM_NAME, &str);
2441    //     std::cout << "Platform[" << i << "] = " << str << std::endl;
2442    // }
2443
2444    return 0;
2445}
Note: See TracBrowser for help on using the browser.