1 // --------------------------------------------------------------------------
4 // Markus Wittmann, 2016-2017
5 // RRZE, University of Erlangen-Nuremberg, Germany
6 // markus.wittmann -at- fau.de or hpc -at- rrze.fau.de
9 // LSS, University of Erlangen-Nuremberg, Germany
11 // This file is part of the Lattice Boltzmann Benchmark Kernels (LbmBenchKernels).
13 // LbmBenchKernels is free software: you can redistribute it and/or modify
14 // it under the terms of the GNU General Public License as published by
15 // the Free Software Foundation, either version 3 of the License, or
16 // (at your option) any later version.
18 // LbmBenchKernels is distributed in the hope that it will be useful,
19 // but WITHOUT ANY WARRANTY; without even the implied warranty of
20 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 // GNU General Public License for more details.
23 // You should have received a copy of the GNU General Public License
24 // along with LbmBenchKernels. If not, see <http://www.gnu.org/licenses/>.
26 // --------------------------------------------------------------------------
27 #include "BenchKernelD3Q19AaVecCommon.h"
42 static void KernelEven(LatticeDesc * ld, KernelData * kd, CaseData * cd);
43 static void KernelOddVecSl(LatticeDesc * ld, KernelData * kd, CaseData * cd);
46 void DumpPdfs(LatticeDesc * ld, KernelData * kd, int zStart, int zStop, int iter, const char * prefix, int dir)
48 int * gDims = kd->GlobalDims;
56 int localZStart = zStart;
57 int localZStop = zStop;
59 if (localZStart == -1) localZStart = 0;
60 if (localZStop == -1) localZStop = gDims[2] - 1;
62 printf("D iter: %d dir: %d %s\n", iter, dir, D3Q19_NAMES[dir]);
64 // for (int dir = 0; dir < 19; ++dir) {
65 for (int z = localZStop; z >= localZStart; --z) {
66 printf("D [%2d][%2d][%s] plane % 2d\n", iter, dir, prefix, z);
68 for(int y = 0; y < nY; ++y) {
69 // for(int y = 2; y < nY - 2; ++y) {
70 printf("D [%2d][%2d][%s] %2d ", iter, dir, prefix, y);
72 for(int x = 0; x < nX; ++x) {
74 if (1) { // ld->Lattice[L_INDEX_4(ld->Dims, x, y, z)] != LAT_CELL_OBSTACLE) {
76 #define I(x, y, z, dir) P_INDEX_5(gDims, (x), (y), (z), (dir))
77 pdfs[dir] = kd->PdfsActive[I(x, y, z, dir)];
84 printf("%.16e ", pdfs[dir]);
85 // printf("%08.0f ", pdfs[dir]);
95 void FNAME(D3Q19AaVecSlKernel)(LatticeDesc * ld, KernelData * kd, CaseData * cd)
101 Assert(cd->Omega > 0.0);
102 Assert(cd->Omega < 2.0);
104 KernelDataAa * kda = KDA(kd);
106 PdfT * src = kd->PdfsActive;
108 int maxIterations = cd->MaxIterations;
112 kd->PdfsActive = src;
113 VtkWrite(ld, kd, cd, -1);
118 kd->PdfsActive = src;
119 KernelStatistics(kd, ld, cd, 0);
122 Assert((maxIterations % 2) == 0);
127 #pragma omp parallel default(none) shared(kda, kd, ld, cd, src, maxIterations)
130 for (int iter = 0; iter < maxIterations; iter += 2) {
132 // --------------------------------------------------------------------
134 // --------------------------------------------------------------------
136 X_LIKWID_START("aa-vec-even");
138 KernelEven(ld, kd, cd);
143 X_LIKWID_STOP("aa-vec-even");
145 // Fixup bounce back PDFs.
149 #ifdef INTEL_OPT_DIRECTIVES
152 for (int i = 0; i < kd->nBounceBackPdfs; ++i) {
153 src[kd->BounceBackPdfsSrc[i]] = src[kd->BounceBackPdfsDst[i]];
160 // save current iteration
161 kda->Iteration = iter;
164 kd->PdfsActive = src;
165 KernelAddBodyForce(kd, ld, cd);
169 if (cd->VtkOutput && (iter % cd->VtkModulus) == 0) {
170 kd->PdfsActive = src;
171 VtkWrite(ld, kd, cd, iter);
176 kd->PdfsActive = src;
177 KernelStatistics(kd, ld, cd, iter);
185 // --------------------------------------------------------------------
187 // --------------------------------------------------------------------
189 X_LIKWID_START("aa-vec-odd");
192 KernelOddVecSl(ld, kd, cd);
197 // Stop counters before bounce back. Else computing loop balance will
200 X_LIKWID_STOP("aa-vec-odd");
202 // Fixup bounce back PDFs.
206 #ifdef INTEL_OPT_DIRECTIVES
209 for (int i = 0; i < kd->nBounceBackPdfs; ++i) {
210 src[kd->BounceBackPdfsDst[i]] = src[kd->BounceBackPdfsSrc[i]];
217 // save current iteration
218 kda->Iteration = iter + 1;
221 kd->PdfsActive = src;
222 KernelAddBodyForce(kd, ld, cd);
226 if (cd->VtkOutput && ((iter + 1) % cd->VtkModulus) == 0) {
227 kd->PdfsActive = src;
228 VtkWrite(ld, kd, cd, iter + 1);
233 kd->PdfsActive = src;
234 KernelStatistics(kd, ld, cd, iter + 1);
240 } // for (int iter = 0; ...
248 kd->PdfsActive = src;
249 VtkWrite(ld, kd, cd, maxIterations);
257 static void KernelEven(LatticeDesc * ld, KernelData * kd, CaseData * cd) // {{{
263 Assert(cd->Omega > F(0.0));
264 Assert(cd->Omega < F(2.0));
266 KernelDataAa * kda = KDA(kd);
268 int nX = ld->Dims[0];
269 int nY = ld->Dims[1];
270 int nZ = ld->Dims[2];
272 int * gDims = kd->GlobalDims;
274 int oX = kd->Offsets[0];
275 int oY = kd->Offsets[1];
276 int oZ = kd->Offsets[2];
279 blk[0] = kda->Blk[0];
280 blk[1] = kda->Blk[1];
281 blk[2] = kda->Blk[2];
283 PdfT omega = cd->Omega;
284 PdfT omegaEven = omega;
286 PdfT magicParam = F(1.0) / F(12.0);
287 PdfT omegaOdd = F(1.0) / (F(0.5) + magicParam / (F(1.0) / omega - F(0.5)));
289 const PdfT w_0 = F(1.0) / F( 3.0);
290 const PdfT w_1 = F(1.0) / F(18.0);
291 const PdfT w_2 = F(1.0) / F(36.0);
293 const PdfT w_1_x3 = w_1 * F(3.0); const PdfT w_1_nine_half = w_1 * F(9.0) / F(2.0);
294 const PdfT w_2_x3 = w_2 * F(3.0); const PdfT w_2_nine_half = w_2 * F(9.0) / F(2.0);
297 VPDFT VONE_HALF = VSET(F(0.5));
298 VPDFT VTHREE_HALF = VSET(F(3.0) / F(2.0));
300 VPDFT vw_1_indep, vw_2_indep;
301 VPDFT vw_0 = VSET(w_0);
302 VPDFT vw_1 = VSET(w_1);
303 VPDFT vw_2 = VSET(w_2);
305 VPDFT vw_1_x3 = VSET(w_1_x3);
306 VPDFT vw_2_x3 = VSET(w_2_x3);
307 VPDFT vw_1_nine_half = VSET(w_1_nine_half);
308 VPDFT vw_2_nine_half = VSET(w_2_nine_half);
310 VPDFT vui, vux, vuy, vuz, vdens;
312 VPDFT vevenPart, voddPart, vdir_indep_trm;
314 VPDFT vomegaEven = VSET(omegaEven);
315 VPDFT vomegaOdd = VSET(omegaOdd);
317 VPDFT vpdf_a, vpdf_b;
319 // Declare pdf_N, pdf_E, pdf_S, pdf_W, ...
320 #define X(name, idx, idxinv, x, y, z) VPDFT JOIN(vpdf_,name); PdfT * JOIN(ppdf_,name);
324 PdfT * src = kd->Pdfs[0];
330 nThreads = omp_get_max_threads();
331 threadId = omp_get_thread_num();
334 const int nodesPlane = gDims[1] * gDims[2];
335 const int nodesCol = gDims[2];
337 #define I(x, y, z, dir) P_INDEX_5(gDims, (x), (y), (z), (dir))
339 // TODO: make inline function out of macros.
341 #define IMPLODE(_x, _y, _z) (nodesPlane * (_x) + nodesCol * (_y) + (_z))
342 #define EXPLODE(index, _x, _y, _z) _x = index / (nodesPlane); _y = (index - nodesPlane * (_x)) / nodesCol; _z = index - nodesPlane * (_x) - nodesCol * (_y);
348 int indexStart = IMPLODE(startX, startY, startZ);
349 int indexEnd = IMPLODE(startX + nX - 1, startY + nY - 1, startZ + nZ - 1);
351 // How many cells as multiples of VSIZE do we have (rounded up)?
352 int nVCells = (indexEnd - indexStart + 1 + VSIZE - 1) / VSIZE;
354 int threadStart = nVCells / nThreads * threadId;
355 int threadEnd = nVCells / nThreads * (threadId + 1);
357 if (nVCells % nThreads > threadId) {
358 threadStart += threadId;
359 threadEnd += threadId + 1;
362 threadStart += nVCells % nThreads;
363 threadEnd += nVCells % nThreads;
366 threadStart *= VSIZE;
369 // As threadStart/End is now in the granularity of cells we add the start offset.
370 threadStart += indexStart;
371 threadEnd += indexStart;
373 EXPLODE(threadStart, startX, startY, startZ);
378 #define X(name, idx, idxinv, _x, _y, _z) JOIN(ppdf_,name) = &src[I(startX, startY, startZ, idx)];
382 // printf("e thread %d idx start: %d end: %d thread start: %d end: %d\n",
383 // threadId, indexStart, indexEnd, threadStart, threadEnd);
386 for (int i = threadStart; i < threadEnd; i += VSIZE) { // LOOP aa-vec-sl-even
388 // Load PDFs of local cell: pdf_N = src[I(x, y, z, D3Q19_N)]; ...
389 // #define X(name, idx, idxinv, _x, _y, _z) JOIN(vpdf_,name) = VLDU(&src[I(x, y, z, idx)]);
390 #define X(name, idx, idxinv, _x, _y, _z) JOIN(vpdf_,name) = VLDU(JOIN(ppdf_,name));
395 vux = VSUB(VSUB(VSUB(VSUB(VSUB(VADD(VADD(vpdf_E,VADD(vpdf_NE,vpdf_SE)),VADD(vpdf_TE,vpdf_BE)),vpdf_W),vpdf_NW),vpdf_SW),vpdf_TW),vpdf_BW);
396 vuy = VSUB(VSUB(VSUB(VSUB(VSUB(VADD(VADD(vpdf_N,VADD(vpdf_NE,vpdf_NW)),VADD(vpdf_TN,vpdf_BN)),vpdf_S),vpdf_SE),vpdf_SW),vpdf_TS),vpdf_BS);
397 vuz = VSUB(VSUB(VSUB(VSUB(VSUB(VADD(VADD(vpdf_T,VADD(vpdf_TE,vpdf_TW)),VADD(vpdf_TN,vpdf_TS)),vpdf_B),vpdf_BE),vpdf_BW),vpdf_BN),vpdf_BS);
399 vdens = VADD(VADD(VADD(VADD(VADD(VADD(VADD(VADD(VADD(vpdf_C,VADD(vpdf_N,vpdf_E)),VADD(vpdf_S,vpdf_W)),VADD(vpdf_NE,vpdf_SE)),
400 VADD(vpdf_SW,vpdf_NW)),VADD(vpdf_T,vpdf_TN)),VADD(vpdf_TE,vpdf_TS)),VADD(vpdf_TW,vpdf_B)),
401 VADD(vpdf_BN,vpdf_BE)),VADD(vpdf_BS,vpdf_BW));
403 vdir_indep_trm = VSUB(vdens,VMUL(VADD(VADD(VMUL(vux,vux),VMUL(vuy,vuy)),VMUL(vuz,vuz)),VTHREE_HALF));
405 VSTU(ppdf_C, VSUB(vpdf_C,VMUL(vomegaEven,VSUB(vpdf_C,VMUL(vw_0,vdir_indep_trm)))));
407 vw_1_indep = VMUL(vw_1,vdir_indep_trm);
408 vw_2_indep = VMUL(vw_2,vdir_indep_trm);
410 #if defined(LOOP_1) || defined(LOOP_2)
411 #error Loop macros are not allowed to be defined here.
414 #define LOOP_1(_dir1, _dir2, _vel) \
416 vpdf_a = JOIN(vpdf_,_dir1); \
417 vpdf_b = JOIN(vpdf_,_dir2); \
419 vevenPart = VMUL(vomegaEven, VSUB(VSUB(VMUL(VONE_HALF, VADD(vpdf_a, vpdf_b)), VMUL(vui, VMUL(vui, vw_1_nine_half))), vw_1_indep)); \
420 voddPart = VMUL(vomegaOdd, VSUB( VMUL(VONE_HALF, VSUB(vpdf_a, vpdf_b)), VMUL(vui, vw_1_x3))); \
422 VSTU(JOIN(ppdf_,_dir2), VSUB(VSUB(vpdf_a, vevenPart), voddPart)); \
423 VSTU(JOIN(ppdf_,_dir1), VADD(VSUB(vpdf_b, vevenPart), voddPart));
425 #define LOOP_2(_dir1, _dir2, _expr) \
427 vpdf_a = JOIN(vpdf_,_dir1); \
428 vpdf_b = JOIN(vpdf_,_dir2); \
430 vevenPart = VMUL(vomegaEven, VSUB(VSUB(VMUL(VONE_HALF, VADD(vpdf_a, vpdf_b)), VMUL(vui, VMUL(vui, vw_2_nine_half))), vw_2_indep)); \
431 voddPart = VMUL(vomegaOdd, VSUB( VMUL(VONE_HALF, VSUB(vpdf_a, vpdf_b)), VMUL(vui, vw_2_x3))); \
433 VSTU(JOIN(ppdf_,_dir2), VSUB(VSUB(vpdf_a, vevenPart), voddPart)); \
434 VSTU(JOIN(ppdf_,_dir1), VADD(VSUB(vpdf_b, vevenPart), voddPart));
440 LOOP_2(NW, SE, VSUB(vuy, vux));
441 LOOP_2(NE, SW, VADD(vuy, vux));
442 LOOP_2(TW, BE, VSUB(vuz, vux));
443 LOOP_2(TE, BW, VADD(vuz, vux));
444 LOOP_2(TS, BN, VSUB(vuz, vuy));
445 LOOP_2(TN, BS, VADD(vuz, vuy));
450 #define X(name, idx, idxinv, _x, _y, _z) JOIN(ppdf_,name) += VSIZE;
461 static void KernelOddVecSl(LatticeDesc * ld, KernelData * kd, CaseData * cd) // {{{
467 Assert(cd->Omega > 0.0);
468 Assert(cd->Omega < F(2.0));
470 KernelDataAa * kda = KDA(kd);
472 int nX = ld->Dims[0];
473 int nY = ld->Dims[1];
474 int nZ = ld->Dims[2];
476 int * gDims = kd->GlobalDims;
478 int oX = kd->Offsets[0];
479 int oY = kd->Offsets[1];
480 int oZ = kd->Offsets[2];
483 blk[0] = kda->Blk[0];
484 blk[1] = kda->Blk[1];
485 blk[2] = kda->Blk[2];
487 PdfT omega = cd->Omega;
488 PdfT omegaEven = omega;
490 PdfT magicParam = F(1.0) / F(12.0);
491 PdfT omegaOdd = F(1.0) / (F(0.5) + magicParam / (F(1.0) / omega - F(0.5)));
493 const PdfT w_0 = F(1.0) / F( 3.0);
494 const PdfT w_1 = F(1.0) / F(18.0);
495 const PdfT w_2 = F(1.0) / F(36.0);
497 const PdfT w_1_x3 = w_1 * F(3.0); const PdfT w_1_nine_half = w_1 * F(9.0) / F(2.0);
498 const PdfT w_2_x3 = w_2 * F(3.0); const PdfT w_2_nine_half = w_2 * F(9.0) / F(2.0);
500 VPDFT VONE_HALF = VSET(F(0.5));
501 VPDFT VTHREE_HALF = VSET(F(3.0) / F(2.0));
503 VPDFT vw_1_indep, vw_2_indep;
504 VPDFT vw_0 = VSET(w_0);
505 VPDFT vw_1 = VSET(w_1);
506 VPDFT vw_2 = VSET(w_2);
508 VPDFT vw_1_x3 = VSET(w_1_x3);
509 VPDFT vw_2_x3 = VSET(w_2_x3);
510 VPDFT vw_1_nine_half = VSET(w_1_nine_half);
511 VPDFT vw_2_nine_half = VSET(w_2_nine_half);
513 VPDFT vui, vux, vuy, vuz, vdens;
515 VPDFT vevenPart, voddPart, vdir_indep_trm;
517 VPDFT vomegaEven = VSET(omegaEven);
518 VPDFT vomegaOdd = VSET(omegaOdd);
520 VPDFT vpdf_a, vpdf_b;
522 // Declare pdf_N, pdf_E, pdf_S, pdf_W, ...
523 #define X(name, idx, idxinv, x, y, z) VPDFT JOIN(vpdf_,name); PdfT * JOIN(ppdf_,idx);
527 PdfT * src = kd->Pdfs[0];
533 nThreads = omp_get_max_threads();
534 threadId = omp_get_thread_num();
537 const int nodesPlane = gDims[1] * gDims[2];
538 const int nodesCol = gDims[2];
540 #define I(x, y, z, dir) P_INDEX_5(gDims, (x), (y), (z), (dir))
542 // TODO: make inline function out of macros.
544 #define IMPLODE(_x, _y, _z) (nodesPlane * (_x) + nodesCol * (_y) + (_z))
545 #define EXPLODE(index, _x, _y, _z) _x = index / (nodesPlane); _y = (index - nodesPlane * (_x)) / nodesCol; _z = index - nodesPlane * (_x) - nodesCol * (_y);
551 int indexStart = IMPLODE(startX, startY, startZ);
552 int indexEnd = IMPLODE(startX + nX - 1, startY + nY - 1, startZ + nZ - 1);
554 // How many multiples of VSIZE cells (rounded up) do we have?
555 int nVCells = (indexEnd - indexStart + 1 + VSIZE - 1) / VSIZE;
557 int threadStart = nVCells / nThreads * threadId;
558 int threadEnd = nVCells / nThreads * (threadId + 1);
560 if (nVCells % nThreads > threadId) {
561 threadStart += threadId;
562 threadEnd += threadId + 1;
565 threadStart += nVCells % nThreads;
566 threadEnd += nVCells % nThreads;
569 threadStart *= VSIZE;
572 // As threadStart/End is now in the granularity of cells we add the start offset.
573 threadStart += indexStart;
574 threadEnd += indexStart;
576 EXPLODE(threadStart, startX, startY, startZ);
581 // printf("o thread %d idx start: %d end: %d thread start: %d end: %d\n",
582 // threadId, indexStart, indexEnd, threadStart, threadEnd);
584 #define X(name, idx, idxinv, _x, _y, _z) JOIN(ppdf_,idx) = &src[I(startX + _x, startY + _y, startZ + _z, idx)];
590 #define X(name, idx, idxinv, x, y, z) PdfT * JOIN(ppdf_start_,idx), * JOIN(ppdf_end_,idx);
594 #define X(name, idx, idxinv, _x, _y, _z) JOIN(ppdf_start_,idx) = &src[I(startX + _x, startY + _y, startZ + _z, idx)];
598 #define X(name, idx, idxinv, _x, _y, _z) JOIN(ppdf_end_,idx) = &src[I(startX + nX - 1 + _x, startY + nY - 1 + _y, startZ + nZ - 1 + _z, idx)];
603 #define X(name, idx, idxinv, _x, _y, _z) printf("%2s ppdf_%d = %p (%d %d %d) (%d %d %d)\n", STRINGIFY(name), idx, JOIN(ppdf_,idx), \
604 startX , startY , startZ , startX + _x, startY + _y, startZ + _z);
609 #endif // DEBUG_EXTENDED
612 for (int i = threadStart; i < threadEnd; i += VSIZE) { // LOOP aa-vec-sl-odd
615 #define X(name, idx, idxinv, _x, _y, _z) Assert((unsigned long)(JOIN(ppdf_,idx)) >= (unsigned long)(JOIN(ppdf_start_,idx))); Assert((unsigned long)(JOIN(ppdf_,idx)) <= (unsigned long)(JOIN(ppdf_end_,idx)));
620 #define X(name, idx, idxinv, _x, _y, _z) JOIN(vpdf_,name) = VLDU(JOIN(ppdf_,idxinv));
624 vux = VSUB(VSUB(VSUB(VSUB(VSUB(VADD(VADD(vpdf_E,VADD(vpdf_NE,vpdf_SE)),VADD(vpdf_TE,vpdf_BE)),vpdf_W),vpdf_NW),vpdf_SW),vpdf_TW),vpdf_BW);
625 vuy = VSUB(VSUB(VSUB(VSUB(VSUB(VADD(VADD(vpdf_N,VADD(vpdf_NE,vpdf_NW)),VADD(vpdf_TN,vpdf_BN)),vpdf_S),vpdf_SE),vpdf_SW),vpdf_TS),vpdf_BS);
626 vuz = VSUB(VSUB(VSUB(VSUB(VSUB(VADD(VADD(vpdf_T,VADD(vpdf_TE,vpdf_TW)),VADD(vpdf_TN,vpdf_TS)),vpdf_B),vpdf_BE),vpdf_BW),vpdf_BN),vpdf_BS);
628 vdens = VADD(VADD(VADD(VADD(VADD(VADD(VADD(VADD(VADD(vpdf_C,VADD(vpdf_N,vpdf_E)),VADD(vpdf_S,vpdf_W)),VADD(vpdf_NE,vpdf_SE)),
629 VADD(vpdf_SW,vpdf_NW)),VADD(vpdf_T,vpdf_TN)),VADD(vpdf_TE,vpdf_TS)),VADD(vpdf_TW,vpdf_B)),VADD(vpdf_BN,vpdf_BE)),VADD(vpdf_BS,vpdf_BW));
631 vdir_indep_trm = VSUB(vdens,VMUL(VADD(VADD(VMUL(vux,vux),VMUL(vuy,vuy)),VMUL(vuz,vuz)),VTHREE_HALF));
633 // ppdf_18 is the pointer to the center pdfs.
634 VSTU(ppdf_18, VSUB(vpdf_C,VMUL(vomegaEven,VSUB(vpdf_C,VMUL(vw_0,vdir_indep_trm)))));
636 vw_1_indep = VMUL(vw_1,vdir_indep_trm);
637 vw_2_indep = VMUL(vw_2,vdir_indep_trm);
639 #if defined(LOOP_1) || defined(LOOP_2)
640 #error Loop macros are not allowed to be defined here.
643 #define LOOP_1(_dir1, _dir2, _idx1, _idx2, _vel) \
645 vpdf_a = JOIN(vpdf_,_dir1); \
646 vpdf_b = JOIN(vpdf_,_dir2); \
648 vevenPart = VMUL(vomegaEven, VSUB(VSUB(VMUL(VONE_HALF, VADD(vpdf_a, vpdf_b)), VMUL(vui, VMUL(vui, vw_1_nine_half))), vw_1_indep)); \
649 voddPart = VMUL(vomegaOdd, VSUB( VMUL(VONE_HALF, VSUB(vpdf_a, vpdf_b)), VMUL(vui, vw_1_x3))); \
651 VSTU(JOIN(ppdf_,_idx1), VSUB(VSUB(vpdf_a, vevenPart), voddPart)); \
652 VSTU(JOIN(ppdf_,_idx2), VADD(VSUB(vpdf_b, vevenPart), voddPart));
654 #define LOOP_2(_dir1, _dir2, _idx1, _idx2, _expr) \
656 vpdf_a = JOIN(vpdf_,_dir1); \
657 vpdf_b = JOIN(vpdf_,_dir2); \
659 vevenPart = VMUL(vomegaEven, VSUB(VSUB(VMUL(VONE_HALF, VADD(vpdf_a, vpdf_b)), VMUL(vui, VMUL(vui, vw_2_nine_half))), vw_2_indep)); \
660 voddPart = VMUL(vomegaOdd, VSUB( VMUL(VONE_HALF, VSUB(vpdf_a, vpdf_b)), VMUL(vui, vw_2_x3))); \
662 VSTU(JOIN(ppdf_,_idx1), VSUB(VSUB(vpdf_a, vevenPart), voddPart)); \
663 VSTU(JOIN(ppdf_,_idx2), VADD(VSUB(vpdf_b, vevenPart), voddPart));
666 LOOP_1(N, S, D3Q19_N, D3Q19_S, vuy);
667 LOOP_1(E, W, D3Q19_E, D3Q19_W, vux);
668 LOOP_1(T, B, D3Q19_T, D3Q19_B, vuz);
670 LOOP_2(NW, SE, D3Q19_NW, D3Q19_SE, VSUB(vuy, vux));
671 LOOP_2(NE, SW, D3Q19_NE, D3Q19_SW, VADD(vuy, vux));
672 LOOP_2(TW, BE, D3Q19_TW, D3Q19_BE, VSUB(vuz, vux));
673 LOOP_2(TE, BW, D3Q19_TE, D3Q19_BW, VADD(vuz, vux));
674 LOOP_2(TS, BN, D3Q19_TS, D3Q19_BN, VSUB(vuz, vuy));
675 LOOP_2(TN, BS, D3Q19_TN, D3Q19_BS, VADD(vuz, vuy));
677 #define X(name, idx, idxinv, _x, _y, _z) JOIN(ppdf_,idx) += VSIZE;