src/BenchKernelD3Q19ListAaPvGatherHybrid.c

   1 // --------------------------------------------------------------------------
   2 //
   3 // Copyright
   4 //   Markus Wittmann, 2016-2017
   5 //   RRZE, University of Erlangen-Nuremberg, Germany
   6 //   markus.wittmann -at- fau.de or hpc -at- rrze.fau.de
   7 //
   8 //   Viktor Haag, 2016
   9 //   LSS, University of Erlangen-Nuremberg, Germany
  10 //
  11 //   Michael Hussnaetter, 2017-2018
  12 //   University of Erlangen-Nuremberg, Germany
  13 //   michael.hussnaetter -at- fau.de
  14 //
  15 //  This file is part of the Lattice Boltzmann Benchmark Kernels (LbmBenchKernels).
  16 //
  17 //  LbmBenchKernels is free software: you can redistribute it and/or modify
  18 //  it under the terms of the GNU General Public License as published by
  19 //  the Free Software Foundation, either version 3 of the License, or
  20 //  (at your option) any later version.
  21 //
  22 //  LbmBenchKernels is distributed in the hope that it will be useful,
  23 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
  24 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  25 //  GNU General Public License for more details.
  26 //
  27 //  You should have received a copy of the GNU General Public License
  28 //  along with LbmBenchKernels.  If not, see <http://www.gnu.org/licenses/>.
  29 //
  30 // --------------------------------------------------------------------------
  31 #include "BenchKernelD3Q19ListAaPvGatherHybridCommon.h"
  32
  33 #include "Base.h"
  34 #include "Memory.h"
  35 #include "Vtk.h"
  36 #include "Vector.h"
  37
  38 #include <inttypes.h>
  39 #include <math.h>
  40
  41 #ifdef _OPENMP
  42 #include <omp.h>
  43 #endif
  44
  45 #ifdef LIKWID_PERFMON
  46 #include <likwid.h>
  47 #else
  48 #define LIKWID_MARKER_INIT
  49 #define LIKWID_MARKER_THREADINIT
  50 #define LIKWID_MARKER_SWITCH
  51 #define LIKWID_MARKER_REGISTER(regionTag)
  52 #define LIKWID_MARKER_START(regionTag)
  53 #define LIKWID_MARKER_STOP(regionTag)
  54 #define LIKWID_MARKER_CLOSE
  55 #define LIKWID_MARKER_GET(regionTag, nevents, events, time, count)
  56 #endif
  57
  58 //enable software prefetchting for vectorized gather/scatter loop in odd kernel
  59 #ifndef SOFTWARE_PREFETCH_LOOKAHEAD_L2
  60 #define SOFTWARE_PREFETCH_LOOKAHEAD_L2 (0) //prefetchting X loop iterations ahead
  61 #endif
  62
  63 #ifndef SOFTWARE_PREFETCH_LOOKAHEAD_L1
  64 #define SOFTWARE_PREFETCH_LOOKAHEAD_L1 (0) //prefetchting X loop iterations ahead
  65 #endif
  66
  67 static void KernelEven(LatticeDesc * ld, KernelData * kernelData, CaseData * cd, int * threadIndices);
  68 static void KernelOdd( LatticeDesc * ld, KernelData * kernelData, CaseData * cd);
  69
  70 void FNAME(D3Q19ListAaPvGatherHybridKernel)(LatticeDesc * ld, KernelData * kernelData, CaseData * cd)
  71 {
  72         Assert(ld != NULL);
  73         Assert(kernelData != NULL);
  74         Assert(cd != NULL);
  75
  76         Assert(cd->Omega > 0.0);
  77         Assert(cd->Omega < 2.0);
  78
  79 #if defined(VTK_OUTPUT) || defined(STATISTICS) || defined(VERIFICATION)
  80         KernelData * kd = (KernelData *)kernelData;
  81 #endif
  82         KernelDataList * kdl = KDL(kernelData);
  83
  84         int maxIterations = cd->MaxIterations;
  85         int nFluid = kdl->nFluid;
  86
  87         printf("\n");
  88 #if (SOFTWARE_PREFETCH_LOOKAHEAD_L2 > 0) || (SOFTWARE_PREFETCH_LOOKAHEAD_L1 > 0)
  89         printf("# Software prefetching enabled:\n");
  90         printf("#   Gather/Scatter + Load/Store prefetch lookahead L2: \t%d\n", SOFTWARE_PREFETCH_LOOKAHEAD_L2);
  91         printf("#   Gather/Scatter + Load/Store prefetch lookahead L1: \t%d\n", SOFTWARE_PREFETCH_LOOKAHEAD_L1);
  92 #else
  93         printf("# Software prefetching disabled.\n");
  94 #endif
  95         printf("\n");
  96
  97         int nThreads = 1;
  98 #ifdef _OPENMP
  99         nThreads = omp_get_max_threads();
 100 #endif
 101
 102         int * threadStartIndices = (int *)malloc(sizeof(int) * (nThreads + 1));
 103
 104         for (int i = 0; i < nThreads; ++i) {
 105                 threadStartIndices[i] = i * (nFluid / nThreads) + MinI(i, nFluid % nThreads);
 106         }
 107
 108         threadStartIndices[nThreads] = nFluid;
 109
 110 #ifdef VTK_OUTPUT
 111         if (cd->VtkOutput) {
 112                 kd->PdfsActive = kd->Pdfs[0];
 113                 VtkWrite(ld, kd, cd, -1);
 114         }
 115 #endif
 116
 117 #ifdef STATISTICS
 118         kd->PdfsActive = kd->Pdfs[0];
 119         KernelStatistics(kd, ld, cd, 0);
 120 #endif
 121
 122         LIKWID_MARKER_INIT;
 123
 124         // TODO: outer openmp parallel
 125
 126         X_KERNEL_START(kernelData);
 127
 128         LIKWID_MARKER_START("OuterLoop");
 129
 130         for(int iter = 0; iter < maxIterations; iter += 2) {
 131
 132                 // even time step
 133
 134 #ifdef _OPENMP
 135 #pragma omp parallel default(none) shared(ld, kernelData, cd, threadStartIndices)
 136 #endif
 137                 {
 138                         //LIKWID_MARKER_START("KernelEven");
 139                         KernelEven(ld, kernelData, cd, threadStartIndices);
 140                         //LIKWID_MARKER_STOP("KernelEven");
 141                 }
 142
 143
 144 #ifdef VERIFICATION
 145                 kdl->Iteration = iter;
 146                 kd->PdfsActive = kd->Pdfs[0];
 147                 KernelAddBodyForce(kd, ld, cd);
 148 #endif
 149
 150                 // odd time step
 151
 152 #ifdef _OPENMP
 153 #pragma omp parallel default(none) shared(ld, kernelData, cd, threadStartIndices)
 154 #endif
 155                 {
 156                         //LIKWID_MARKER_START("KernelOdd");
 157                         KernelOdd(ld, kernelData, cd);
 158                         //LIKWID_MARKER_STOP("KernelOdd");
 159                 }
 160
 161
 162 #ifdef VERIFICATION
 163                 kdl->Iteration = iter + 1;
 164                 kd->PdfsActive = kd->Pdfs[0];
 165                 KernelAddBodyForce(kd, ld, cd);
 166 #endif
 167
 168 #ifdef VTK_OUTPUT
 169                 if (cd->VtkOutput && (iter % cd->VtkModulus) == 0) {
 170                         kdl->Iteration = iter + 1;
 171                         kd->PdfsActive = kd->Pdfs[0];
 172                         VtkWrite(ld, kd, cd, iter);
 173                 }
 174 #endif
 175
 176 #ifdef STATISTICS
 177                 kdl->Iteration = iter + 1;
 178                 kd->PdfsActive = kd->Pdfs[0];
 179                 KernelStatistics(kd, ld, cd, iter);
 180 #endif
 181
 182         } // for (int iter = 0; ...
 183         LIKWID_MARKER_STOP("OuterLoop");
 184
 185         X_KERNEL_END(kernelData);
 186
 187 #ifdef VTK_OUTPUT
 188         if (cd->VtkOutput) {
 189                 kd->PdfsActive = kd->Pdfs[0];
 190                 VtkWrite(ld, kd, cd, maxIterations);
 191         }
 192 #endif
 193
 194 #ifdef STATISTICS
 195         kd->PdfsActive = kd->Pdfs[0];
 196         KernelStatistics(kd, ld, cd, maxIterations);
 197 #endif
 198
 199         LIKWID_MARKER_CLOSE;
 200         free(threadStartIndices);
 201
 202         return;
 203 }
 204
 205 static void KernelEven(LatticeDesc * ld, KernelData * kernelData, CaseData * cd, int * threadIndices)
 206 {
 207         Assert(ld != NULL);
 208         Assert(kernelData != NULL);
 209         Assert(cd != NULL);
 210
 211         Assert(cd->Omega > 0.0);
 212         Assert(cd->Omega < 2.0);
 213
 214         KernelData * kd = (KernelData *)kernelData;
 215         KernelDataList * kdl = KDL(kernelData);
 216
 217         PdfT omega = cd->Omega;
 218         PdfT omegaEven = omega;
 219
 220         PdfT magicParam = 1.0 / 12.0;
 221         PdfT omegaOdd = 1.0 / (0.5 + magicParam / (1.0 / omega - 0.5));
 222
 223         PdfT evenPart = 0.0;
 224         PdfT oddPart = 0.0;
 225         PdfT dir_indep_trm = 0.0;
 226
 227         const PdfT w_0 = 1.0 /  3.0;
 228         const PdfT w_1 = 1.0 / 18.0;
 229         const PdfT w_2 = 1.0 / 36.0;
 230
 231         const PdfT w_1_x3 = w_1 * 3.0;  const PdfT w_1_nine_half = w_1 * 9.0 / 2.0;     PdfT w_1_indep = 0.0;
 232         const PdfT w_2_x3 = w_2 * 3.0;  const PdfT w_2_nine_half = w_2 * 9.0 / 2.0;     PdfT w_2_indep = 0.0;
 233
 234         PdfT ux, uy, uz, ui;
 235         PdfT dens;
 236
 237         VPDFT VONE_HALF = VSET(0.5);
 238         VPDFT VTHREE_HALF = VSET(3.0 / 2.0);
 239
 240         VPDFT vw_1_indep, vw_2_indep;
 241         VPDFT vw_0 = VSET(w_0);
 242         VPDFT vw_1 = VSET(w_1);
 243         VPDFT vw_2 = VSET(w_2);
 244
 245         VPDFT vw_1_x3 = VSET(w_1_x3);
 246         VPDFT vw_2_x3 = VSET(w_2_x3);
 247         VPDFT vw_1_nine_half = VSET(w_1_nine_half);
 248         VPDFT vw_2_nine_half = VSET(w_2_nine_half);
 249
 250         VPDFT vui, vux, vuy, vuz, vdens;
 251
 252         VPDFT vevenPart, voddPart, vdir_indep_trm;
 253
 254         VPDFT vomegaEven = VSET(omegaEven);
 255         VPDFT vomegaOdd  = VSET(omegaOdd);
 256
 257         // Declare pdf_N, pdf_E, pdf_S, pdf_W, ...
 258         #define X(name, idx, idxinv, x, y, z) \
 259                 PdfT JOIN(pdf_,name); \
 260                 PdfT * JOIN(ppdf_,name); \
 261                 VPDFT JOIN(vpdf_,name);
 262                 D3Q19_LIST
 263         #undef X
 264
 265         PdfT * src = kd->Pdfs[0];
 266
 267         int nCells = kdl->nCells;
 268
 269         int threadId = 0;
 270 #ifdef _OPENMP
 271         threadId =  omp_get_thread_num();
 272 #endif
 273
 274         int nFluidThread = threadIndices[threadId + 1] - threadIndices[threadId];
 275         int nFluidVec    = nFluidThread - (nFluidThread % VSIZE);
 276
 277         int indexStartVec = threadIndices[threadId];
 278         int indexStopVec  = threadIndices[threadId] + nFluidVec;
 279         int indexStop     = threadIndices[threadId] + nFluidThread;
 280
 281         #define I(index, dir)   P_INDEX_3((nCells), (index), (dir))
 282
 283         #define X(name, idx, idxinv, _x, _y, _z)        JOIN(ppdf_,name) = &(src[I(indexStartVec, idx)]);
 284                         D3Q19_LIST
 285         #undef X
 286
 287         for (int index = indexStartVec; index < indexStopVec; index += VSIZE) {
 288
 289                 #if (SOFTWARE_PREFETCH_LOOKAHEAD_L2 > 0)
 290                         #define X(name, idx, idxinv, _x, _y, _z) _mm_prefetch((char const *)(JOIN(ppdf_,name) + SOFTWARE_PREFETCH_LOOKAHEAD_L2 * VSIZE), _MM_HINT_T1);
 291                                 D3Q19_LIST
 292                         #undef X
 293                 #endif
 294
 295                 #if (SOFTWARE_PREFETCH_LOOKAHEAD_L1 > 0)
 296                         #define X(name, idx, idxinv, _x, _y, _z) _mm_prefetch((char const *)(JOIN(ppdf_,name) + SOFTWARE_PREFETCH_LOOKAHEAD_L1 * VSIZE), _MM_HINT_T0);
 297                                 D3Q19_LIST
 298                         #undef X
 299                 #endif
 300
 301                 #define X(name, idx, idxinv, _x, _y, _z)        JOIN(vpdf_,name) = VLDU(JOIN(ppdf_,name));
 302                                 D3Q19_LIST
 303                 #undef X
 304
 305                 //vux = vpdf_E + vpdf_NE + vpdf_SE + vpdf_TE + vpdf_BE -
 306                 //           vpdf_W - vpdf_NW - vpdf_SW - vpdf_TW - vpdf_BW;
 307                 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);
 308                 //vuy = vpdf_N + vpdf_NE + vpdf_NW + vpdf_TN + vpdf_BN -
 309                 //           vpdf_S - vpdf_SE - vpdf_SW - vpdf_TS - vpdf_BS;
 310                 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);
 311                 //vuz = vpdf_T + vpdf_TE + vpdf_TW + vpdf_TN + vpdf_TS -
 312                 //           vpdf_B - vpdf_BE - vpdf_BW - vpdf_BN - vpdf_BS;
 313                 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);
 314
 315                 //vdens = vpdf_C +
 316                 //          vpdf_N  + vpdf_E  + vpdf_S  + vpdf_W  +
 317                 //          vpdf_NE + vpdf_SE + vpdf_SW + vpdf_NW +
 318                 //          vpdf_T  + vpdf_TN + vpdf_TE + vpdf_TS + vpdf_TW +
 319                 //          vpdf_B  + vpdf_BN + vpdf_BE + vpdf_BS + vpdf_BW;
 320                 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)),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));
 321
 322                 //vdir_indep_trm = vdens - (vux * vux + vuy * vuy + vuz * vuz) * VTHREE_HALF;
 323                 vdir_indep_trm = VSUB(vdens,VMUL(VADD(VADD(VMUL(vux,vux),VMUL(vuy,vuy)),VMUL(vuz,vuz)),VTHREE_HALF));
 324
 325                 //src[I(index, D3Q19_C)]  =[UA] vpdf_C - vomegaEven * (vpdf_C - vw_0 * vdir_indep_trm);
 326                 VSTU(ppdf_C,VSUB(vpdf_C,VMUL(vomegaEven,VSUB(vpdf_C,VMUL(vw_0,vdir_indep_trm)))));
 327
 328                 //vw_1_indep = vw_1 * vdir_indep_trm;
 329                 vw_1_indep = VMUL(vw_1,vdir_indep_trm);
 330
 331                 #define COLLIDE_AA_V(tmpVui, dir1, dir2) \
 332                         vui = tmpVui; \
 333                         vevenPart = VMUL(vomegaEven,VSUB(VSUB(VMUL(VONE_HALF,VADD(JOIN(vpdf_,dir1),JOIN(vpdf_,dir2))),VMUL(vui,VMUL(vui,vw_1_nine_half))),vw_1_indep));\
 334                         voddPart = VMUL(vomegaOdd,VSUB(VMUL(VONE_HALF,VSUB(JOIN(vpdf_,dir1),JOIN(vpdf_,dir2))),VMUL(vui,vw_1_x3)));\
 335                         VSTU(JOIN(ppdf_,dir2),VSUB(VSUB(JOIN(vpdf_,dir1),vevenPart),voddPart));\
 336                         VSTU(JOIN(ppdf_,dir1),VADD(VSUB(JOIN(vpdf_,dir2),vevenPart),voddPart));
 337
 338                 COLLIDE_AA_V(vuy, N, S)
 339                 COLLIDE_AA_V(vux, E, W)
 340                 COLLIDE_AA_V(vuz, T, B)
 341
 342                 #undef COLLIDE_AA_V
 343
 344                 //vw_2_indep = vw_2 * vdir_indep_trm;
 345                 vw_2_indep = VMUL(vw_2,vdir_indep_trm);
 346
 347                 // collide axis unaligned pdfs vectorized
 348                 #define COLLIDE_UA_V(tmpVui, dir1, dir2) \
 349                         vui = tmpVui; \
 350                         vevenPart = VMUL(vomegaEven,VSUB(VSUB(VMUL(VONE_HALF,VADD(JOIN(vpdf_,dir1),JOIN(vpdf_,dir2))),VMUL(vui,VMUL(vui,vw_2_nine_half))),vw_2_indep));\
 351                         voddPart = VMUL(vomegaOdd,VSUB(VMUL(VONE_HALF,VSUB(JOIN(vpdf_,dir1),JOIN(vpdf_,dir2))),VMUL(vui,vw_2_x3)));\
 352                         VSTU(JOIN(ppdf_,dir2),VSUB(VSUB(JOIN(vpdf_,dir1),vevenPart),voddPart)); \
 353                         VSTU(JOIN(ppdf_,dir1),VADD(VSUB(JOIN(vpdf_,dir2),vevenPart),voddPart));
 354
 355                 COLLIDE_UA_V(VSUB(vuy,vux), NW, SE)
 356                 COLLIDE_UA_V(VADD(vux,vuy), NE, SW)
 357                 COLLIDE_UA_V(VSUB(vuz,vux), TW, BE)
 358                 COLLIDE_UA_V(VADD(vux,vuz), TE, BW)
 359                 COLLIDE_UA_V(VSUB(vuz,vuy), TS, BN)
 360                 COLLIDE_UA_V(VADD(vuy,vuz), TN, BS)
 361
 362                 #undef COLLIDE_UA_V
 363
 364                 #define X(name, idx, idxinv, _x, _y, _z)        JOIN(ppdf_,name) +=VSIZE;
 365                         D3Q19_LIST
 366                 #undef X
 367         } // loop over fluid nodes
 368
 369         for (int index = indexStopVec; index < indexStop; ++index) {
 370
 371                 #define X(name, idx, idxinv, _x, _y, _z)        JOIN(pdf_,name) = *(JOIN(ppdf_,name));
 372                         D3Q19_LIST
 373                 #undef X
 374
 375                 ux = pdf_E + pdf_NE + pdf_SE + pdf_TE + pdf_BE -
 376                         pdf_W - pdf_NW - pdf_SW - pdf_TW - pdf_BW;
 377                 uy = pdf_N + pdf_NE + pdf_NW + pdf_TN + pdf_BN -
 378                         pdf_S - pdf_SE - pdf_SW - pdf_TS - pdf_BS;
 379                 uz = pdf_T + pdf_TE + pdf_TW + pdf_TN + pdf_TS -
 380                         pdf_B - pdf_BE - pdf_BW - pdf_BN - pdf_BS;
 381
 382                 dens = pdf_C +
 383                         pdf_N  + pdf_E  + pdf_S  + pdf_W  +
 384                         pdf_NE + pdf_SE + pdf_SW + pdf_NW +
 385                         pdf_T  + pdf_TN + pdf_TE + pdf_TS + pdf_TW +
 386                         pdf_B  + pdf_BN + pdf_BE + pdf_BS + pdf_BW;
 387
 388                 dir_indep_trm = dens - (ux * ux + uy * uy + uz * uz)*3.0/2.0;
 389
 390                 // direction: w_0
 391                 *ppdf_C  = pdf_C - omegaEven*(pdf_C - w_0*dir_indep_trm);
 392
 393                 // direction: w_1
 394                 w_1_indep = w_1*dir_indep_trm;
 395
 396                 #define COLLIDE_AA_S(tmpUi, dir1, dir2) \
 397                         ui = tmpUi; \
 398                         evenPart = omegaEven * (0.5 * (JOIN(pdf_,dir1) + JOIN(pdf_,dir2)) - ui * ui * w_1_nine_half - w_1_indep); \
 399                         oddPart  = omegaOdd  * (0.5 * (JOIN(pdf_,dir1) - JOIN(pdf_,dir2)) - ui * w_1_x3); \
 400                         *(JOIN(ppdf_,dir2))  = JOIN(pdf_,dir1) - evenPart - oddPart; \
 401                         *(JOIN(ppdf_,dir1))  = JOIN(pdf_,dir2) - evenPart + oddPart;
 402
 403                 COLLIDE_AA_S(uy, N, S)
 404                 COLLIDE_AA_S(ux, E, W)
 405                 COLLIDE_AA_S(uz, T, B)
 406
 407                 #undef COLLIDE_AA_S
 408
 409                 // direction: w_2
 410                 w_2_indep = w_2*dir_indep_trm;
 411
 412                 #define COLLIDE_UA_S(tmpUi, dir1, dir2) \
 413                         ui = tmpUi; \
 414                         evenPart = omegaEven * (0.5 * (JOIN(pdf_,dir1) + JOIN(pdf_,dir2)) - ui * ui * w_2_nine_half - w_2_indep); \
 415                         oddPart  = omegaOdd  * (0.5 * (JOIN(pdf_,dir1) - JOIN(pdf_,dir2)) - ui * w_2_x3); \
 416                         *(JOIN(ppdf_,dir2)) = JOIN(pdf_,dir1) - evenPart - oddPart; \
 417                         *(JOIN(ppdf_,dir1)) = JOIN(pdf_,dir2) - evenPart + oddPart;
 418
 419                 COLLIDE_UA_S((-ux + uy), NW, SE)
 420                 COLLIDE_UA_S(( ux + uy), NE, SW)
 421                 COLLIDE_UA_S((-ux + uz), TW, BE)
 422                 COLLIDE_UA_S(( ux + uz), TE, BW)
 423                 COLLIDE_UA_S((-uy + uz), TS, BN)
 424                 COLLIDE_UA_S(( uy + uz), TN, BS)
 425
 426                 #undef COLLIDE_UA_S
 427
 428                 #define X(name, idx, idxinv, _x, _y, _z)        JOIN(ppdf_,name)++;
 429                         D3Q19_LIST
 430                 #undef X
 431         } // loop over fluid nodes
 432
 433         #undef I
 434
 435         return;
 436 }
 437
 438
 439 static void KernelOdd(LatticeDesc * ld, KernelData * kernelData, CaseData * cd)
 440 {
 441
 442         Assert(ld != NULL);
 443         Assert(kernelData != NULL);
 444         Assert(cd != NULL);
 445
 446         Assert(cd->Omega > 0.0);
 447         Assert(cd->Omega < 2.0);
 448
 449         KernelData * kd = (KernelData *)kernelData;
 450         KernelDataList * kdl = KDL(kernelData);
 451         KernelDataListRia * kdlr = KDLR(kernelData);
 452
 453         PdfT omega = cd->Omega;
 454         PdfT omegaEven = omega;
 455
 456         PdfT magicParam = 1.0 / 12.0;
 457         PdfT omegaOdd = 1.0 / (0.5 + magicParam / (1.0 / omega - 0.5));
 458
 459         PdfT evenPart = 0.0;
 460         PdfT oddPart = 0.0;
 461         PdfT dir_indep_trm = 0.0;
 462
 463         const PdfT w_0 = 1.0 /  3.0;
 464         const PdfT w_1 = 1.0 / 18.0;
 465         const PdfT w_2 = 1.0 / 36.0;
 466
 467         const PdfT w_1_x3 = w_1 * 3.0;  const PdfT w_1_nine_half = w_1 * 9.0 / 2.0;     PdfT w_1_indep = 0.0;
 468         const PdfT w_2_x3 = w_2 * 3.0;  const PdfT w_2_nine_half = w_2 * 9.0 / 2.0;     PdfT w_2_indep = 0.0;
 469
 470         PdfT ux, uy, uz, ui;
 471         PdfT dens;
 472
 473         VPDFT VONE_HALF = VSET(0.5);
 474         VPDFT VTHREE_HALF = VSET(3.0 / 2.0);
 475
 476         VPDFT vw_1_indep, vw_2_indep;
 477         VPDFT vw_0 = VSET(w_0);
 478         VPDFT vw_1 = VSET(w_1);
 479         VPDFT vw_2 = VSET(w_2);
 480
 481         VPDFT vw_1_x3 = VSET(w_1_x3);
 482         VPDFT vw_2_x3 = VSET(w_2_x3);
 483         VPDFT vw_1_nine_half = VSET(w_1_nine_half);
 484         VPDFT vw_2_nine_half = VSET(w_2_nine_half);
 485
 486         VPDFT vux, vuy, vuz, vui;
 487         VPDFT vdens;
 488
 489         VPDFT vevenPart, voddPart, vdir_indep_trm;
 490
 491         VPDFT vomegaEven = VSET(omegaEven);
 492         VPDFT vomegaOdd  = VSET(omegaOdd);
 493
 494
 495         // Declare pdf_N, pdf_E, pdf_S, pdf_W, ...
 496         #define X(name, idx, idxinv, x, y, z) \
 497                 PdfT JOIN(pdf_,name); \
 498                 PdfT * JOIN(ppdf_,name) = NULL; \
 499                 VPDFT JOIN(vpdf_,name);
 500                 D3Q19_LIST
 501         #undef X
 502         #define X(name, idx, idxinv, x, y, z)   __m256i JOIN(vgatheridx_,name);
 503                 D3Q19_LIST_WO_C
 504         #undef X
 505
 506         PdfT * src = kd->Pdfs[0];
 507
 508         int nCells = kdl->nCells;
 509
 510         int adjListIndex;
 511         uint32_t * adjList = kdl->AdjList;
 512
 513         int * oddKernelThreadStartIndices = kdlr->oddKernelThreadStartIndices;
 514         int * loopStartIndices = kdlr->loopStartIndices;
 515
 516         int threadId = 0;
 517
 518 #ifdef _OPENMP
 519         threadId = omp_get_thread_num();
 520 #endif
 521
 522         int threadStartIndex = oddKernelThreadStartIndices[threadId + 0];
 523         int threadStopIndex  = oddKernelThreadStartIndices[threadId + 1];
 524
 525         int indexStart    = loopStartIndices[threadStartIndex];
 526         int indexStop     = loopStartIndices[threadStopIndex ];
 527
 528         int consecStartIndex = threadStartIndex + 1;
 529         int consecStopIndex  = threadStopIndex  - 1;
 530
 531         #define I(index, dir)   P_INDEX_3((nCells), (index), (dir))
 532         #define ADJ_LIST(dir) adjList[adjListIndex + (dir * VSIZE)]
 533
 534         // scalar peel loop
 535         #pragma novector
 536         #pragma loop_count max=7
 537         for (int index = indexStart; index < loopStartIndices[consecStartIndex]; ++index) {
 538
 539                 adjListIndex = (index - (index % VSIZE)) * N_D3Q19_IDX + (index % VSIZE);
 540                 #define X(name, idx, idxinv, _x, _y, _z)        JOIN(ppdf_,name) = &(src[ADJ_LIST(idxinv)]);
 541                         D3Q19_LIST_WO_C
 542                 #undef X
 543                 #define X(name, idx, idxinv, _x, _y, _z)        JOIN(pdf_,name) = *(JOIN(ppdf_,name));
 544                         D3Q19_LIST_WO_C
 545                 #undef X
 546
 547                 ppdf_C = &(src[P_INDEX_3(nCells, index, D3Q19_C)]);
 548                 pdf_C = *ppdf_C;
 549
 550                 // macroscopic values
 551                 ux = pdf_E + pdf_NE + pdf_SE + pdf_TE + pdf_BE -
 552                         pdf_W - pdf_NW - pdf_SW - pdf_TW - pdf_BW;
 553                 uy = pdf_N + pdf_NE + pdf_NW + pdf_TN + pdf_BN -
 554                         pdf_S - pdf_SE - pdf_SW - pdf_TS - pdf_BS;
 555                 uz = pdf_T + pdf_TE + pdf_TW + pdf_TN + pdf_TS -
 556                         pdf_B - pdf_BE - pdf_BW - pdf_BN - pdf_BS;
 557
 558                 dens = pdf_C +
 559                         pdf_N  + pdf_E  + pdf_S  + pdf_W  +
 560                         pdf_NE + pdf_SE + pdf_SW + pdf_NW +
 561                         pdf_T  + pdf_TN + pdf_TE + pdf_TS + pdf_TW +
 562                         pdf_B  + pdf_BN + pdf_BE + pdf_BS + pdf_BW;
 563
 564                 dir_indep_trm = dens - (ux * ux + uy * uy + uz * uz)*3.0/2.0;
 565
 566                 // collide direction: w_0
 567                 *ppdf_C  = pdf_C - omegaEven * (pdf_C - w_0 * dir_indep_trm);
 568
 569                 // collide direction: w_1
 570                 w_1_indep = w_1 * dir_indep_trm;
 571
 572                 #define COLLIDE_AA_S(tmpUi, dir1, dir2) \
 573                         ui = tmpUi; \
 574                         evenPart = omegaEven * (0.5 * (JOIN(pdf_,dir1) + JOIN(pdf_,dir2)) - ui * ui * w_1_nine_half - w_1_indep); \
 575                         oddPart  = omegaOdd  * (0.5 * (JOIN(pdf_,dir1) - JOIN(pdf_,dir2)) - ui * w_1_x3); \
 576                         *(JOIN(ppdf_,dir2))  = JOIN(pdf_,dir1) - evenPart - oddPart; \
 577                         *(JOIN(ppdf_,dir1))  = JOIN(pdf_,dir2) - evenPart + oddPart;
 578
 579                 COLLIDE_AA_S(uy, N, S)
 580                 COLLIDE_AA_S(ux, E, W)
 581                 COLLIDE_AA_S(uz, T, B)
 582
 583                 #undef COLLIDE_AA_S
 584
 585                 // collide direction: w_2
 586                 w_2_indep = w_2 * dir_indep_trm;
 587
 588                 #define COLLIDE_UA_S(tmpUi, dir1, dir2) \
 589                         ui = tmpUi; \
 590                         evenPart = omegaEven * (0.5 * (JOIN(pdf_,dir1) + JOIN(pdf_,dir2)) - ui * ui * w_2_nine_half - w_2_indep); \
 591                         oddPart  = omegaOdd  * (0.5 * (JOIN(pdf_,dir1) - JOIN(pdf_,dir2)) - ui * w_2_x3); \
 592                         *(JOIN(ppdf_,dir2)) = JOIN(pdf_,dir1) - evenPart - oddPart; \
 593                         *(JOIN(ppdf_,dir1)) = JOIN(pdf_,dir2) - evenPart + oddPart;
 594
 595                 COLLIDE_UA_S((-ux + uy), NW, SE)
 596                 COLLIDE_UA_S(( ux + uy), NE, SW)
 597                 COLLIDE_UA_S((-ux + uz), TW, BE)
 598                 COLLIDE_UA_S(( ux + uz), TE, BW)
 599                 COLLIDE_UA_S((-uy + uz), TS, BN)
 600                 COLLIDE_UA_S(( uy + uz), TN, BS)
 601
 602                 #undef COLLIDE_UA_S
 603
 604         } // scalar peel loop
 605
 606         // vectorized loop
 607         for (int consecIndex = consecStartIndex; consecIndex < consecStopIndex; /* increment in loop body */) {
 608
 609                 int index = loopStartIndices[consecIndex];
 610                 adjListIndex = (index - (index % VSIZE)) * N_D3Q19_IDX + (index % VSIZE);
 611                 #define X(name, idx, idxinv, _x, _y, _z) JOIN(ppdf_,name) = &(src[ADJ_LIST(idxinv)]);
 612                         D3Q19_LIST_WO_C
 613                 #undef X
 614                 ppdf_C = &(src[P_INDEX_3(nCells, index, D3Q19_C)]);
 615
 616                 // loop for consecutive access pattern -> load & store
 617                 ++consecIndex;
 618
 619                 for (; index < loopStartIndices[consecIndex]; index+=VSIZE){
 620
 621                         #if (SOFTWARE_PREFETCH_LOOKAHEAD_L2 > 0)
 622                                 #define X(name, idx, idxinv, _x, _y, _z) _mm_prefetch((char const *)(JOIN(ppdf_,name) + SOFTWARE_PREFETCH_LOOKAHEAD_L2 * VSIZE), _MM_HINT_T1);
 623                                         D3Q19_LIST
 624                                 #undef X
 625                         #endif
 626                         #if (SOFTWARE_PREFETCH_LOOKAHEAD_L1 > 0)
 627                                 #define X(name, idx, idxinv, _x, _y, _z) _mm_prefetch((char const *)(JOIN(ppdf_,name) + SOFTWARE_PREFETCH_LOOKAHEAD_L1 * VSIZE), _MM_HINT_T0);
 628                                         D3Q19_LIST
 629                                 #undef X
 630                         #endif
 631
 632                         #define X(name, idx, idxinv, _x, _y, _z)        JOIN(vpdf_,name) = VLDU(JOIN(ppdf_,name));
 633                                 D3Q19_LIST
 634                         #undef X
 635
 636                         // macroscopic values
 637                         //vux = vpdf_E + vpdf_NE + vpdf_SE + vpdf_TE + vpdf_BE -
 638                         //      vpdf_W - vpdf_NW - vpdf_SW - vpdf_TW - vpdf_BW;
 639                         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);
 640                         //vuy = vpdf_N + vpdf_NE + vpdf_NW + vpdf_TN + vpdf_BN -
 641                         //      vpdf_S - vpdf_SE - vpdf_SW - vpdf_TS - vpdf_BS;
 642                         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);
 643                         //vuz = vpdf_T + vpdf_TE + vpdf_TW + vpdf_TN + vpdf_TS -
 644                         //      vpdf_B - vpdf_BE - vpdf_BW - vpdf_BN - vpdf_BS;
 645                         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);
 646
 647                         //vdens = vpdf_C +
 648                         //        vpdf_N  + vpdf_E  + vpdf_S  + vpdf_W  +
 649                         //        vpdf_NE + vpdf_SE + vpdf_SW + vpdf_NW +
 650                         //        vpdf_T  + vpdf_TN + vpdf_TE + vpdf_TS + vpdf_TW +
 651                         //        vpdf_B  + vpdf_BN + vpdf_BE + vpdf_BS + vpdf_BW;
 652                         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)),
 653                                                                                 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));
 654
 655                         //vdir_indep_trm = vdens - (vux * vux + vuy * vuy + vuz * vuz) * VTHREE_HALF;
 656                         vdir_indep_trm = VSUB(vdens,VMUL(VADD(VADD(VMUL(vux,vux),VMUL(vuy,vuy)),VMUL(vuz,vuz)),VTHREE_HALF));
 657
 658                         // collide direction: w_0
 659                         //src[I(index, D3Q19_C)]  =[UA] vpdf_C - vomegaEven * (vpdf_C - vw_0 * vdir_indep_trm);
 660                         VSTU(ppdf_C,VSUB(vpdf_C,VMUL(vomegaEven,VSUB(vpdf_C,VMUL(vw_0,vdir_indep_trm)))));
 661
 662                         // collide direction: w_1
 663                         vw_1_indep = VMUL(vw_1,vdir_indep_trm);
 664
 665                         #define COLLIDE_AA_V(tmpVui, dir1, dir2) \
 666                                 vui = tmpVui; \
 667                                 vevenPart = VMUL(vomegaEven,VSUB(VSUB(VMUL(VONE_HALF,VADD(JOIN(vpdf_,dir1),JOIN(vpdf_,dir2))),VMUL(vui,VMUL(vui,vw_1_nine_half))),vw_1_indep));\
 668                                 voddPart = VMUL(vomegaOdd,VSUB(VMUL(VONE_HALF,VSUB(JOIN(vpdf_,dir1),JOIN(vpdf_,dir2))),VMUL(vui,vw_1_x3)));\
 669                                 VSTU(JOIN(ppdf_,dir2),VSUB(VSUB(JOIN(vpdf_,dir1),vevenPart),voddPart));\
 670                                 VSTU(JOIN(ppdf_,dir1),VADD(VSUB(JOIN(vpdf_,dir2),vevenPart),voddPart));
 671
 672                         COLLIDE_AA_V(vuy, N, S)
 673                         COLLIDE_AA_V(vux, E, W)
 674                         COLLIDE_AA_V(vuz, T, B)
 675
 676                         #undef COLLIDE_AA_V
 677
 678                         // collide direction: w_2
 679                         vw_2_indep = VMUL(vw_2,vdir_indep_trm);
 680
 681                         #define COLLIDE_UA_V(tmpVui, dir1, dir2) \
 682                                 vui = tmpVui; \
 683                                 vevenPart = VMUL(vomegaEven,VSUB(VSUB(VMUL(VONE_HALF,VADD(JOIN(vpdf_,dir1),JOIN(vpdf_,dir2))),VMUL(vui,VMUL(vui,vw_2_nine_half))),vw_2_indep));\
 684                                 voddPart = VMUL(vomegaOdd,VSUB(VMUL(VONE_HALF,VSUB(JOIN(vpdf_,dir1),JOIN(vpdf_,dir2))),VMUL(vui,vw_2_x3)));\
 685                                 VSTU(JOIN(ppdf_,dir2),VSUB(VSUB(JOIN(vpdf_,dir1),vevenPart),voddPart)); \
 686                                 VSTU(JOIN(ppdf_,dir1),VADD(VSUB(JOIN(vpdf_,dir2),vevenPart),voddPart));
 687
 688                         COLLIDE_UA_V(VSUB(vuy,vux), NW, SE)
 689                         COLLIDE_UA_V(VADD(vux,vuy), NE, SW)
 690                         COLLIDE_UA_V(VSUB(vuz,vux), TW, BE)
 691                         COLLIDE_UA_V(VADD(vux,vuz), TE, BW)
 692                         COLLIDE_UA_V(VSUB(vuz,vuy), TS, BN)
 693                         COLLIDE_UA_V(VADD(vuy,vuz), TN, BS)
 694
 695                         #undef COLLIDE_UA_V
 696
 697                         //increment pointer
 698                         #define X(name, idx, idxinv, _x, _y, _z) JOIN(ppdf_,name) += VSIZE;
 699                                 D3Q19_LIST
 700                         #undef X
 701                 }
 702
 703                 // loop for irregular access pattern -> gather & scatter
 704                 index = loopStartIndices[consecIndex];
 705                 ++consecIndex;
 706                 for (; index < loopStartIndices[consecIndex]; index+=VSIZE){
 707
 708
 709                         #if (SOFTWARE_PREFETCH_LOOKAHEAD_L2 > 0)
 710                                 int const indexPrefetchL2 = index + VSIZE * SOFTWARE_PREFETCH_LOOKAHEAD_L2;
 711                                 // make sure that adjList access is never out of bounds since it is an actual memory access and no prefetch
 712                                 if (indexPrefetchL2 < indexStop){
 713                                         adjListIndex = (indexPrefetchL2 - (indexPrefetchL2 % VSIZE)) * N_D3Q19_IDX + (indexPrefetchL2 % VSIZE);
 714                                         #define X(name, idx, idxinv, _x, _y, _z) VPG32(VLIU(&ADJ_LIST(idxinv)), (char const *) src, 8, _MM_HINT_T1);
 715                                                 D3Q19_LIST_WO_C
 716                                         #undef X
 717
 718                                         _mm_prefetch((char const *) &(src[P_INDEX_3(nCells, indexPrefetchL2, D3Q19_C)]), _MM_HINT_T1);
 719                                 }
 720                         #endif
 721
 722                         #if (SOFTWARE_PREFETCH_LOOKAHEAD_L1 > 0)
 723                                 int const indexPrefetchL1 = index + VSIZE * SOFTWARE_PREFETCH_LOOKAHEAD_L1;
 724                                 // make sure that adjList access is never out of bounds since it is an actual memory access and no prefetch
 725                                 if (indexPrefetchL1 < indexStop){
 726                                         adjListIndex = (indexPrefetchL1 - (indexPrefetchL1 % VSIZE)) * N_D3Q19_IDX + (indexPrefetchL1 % VSIZE);
 727                                         #define X(name, idx, idxinv, _x, _y, _z) VPG32(VLIU(&ADJ_LIST(idxinv)), (char const *) src, 8, _MM_HINT_T0);
 728                                                 D3Q19_LIST_WO_C
 729                                         #undef X
 730
 731                                         _mm_prefetch((char const *) &(src[P_INDEX_3(nCells, indexPrefetchL1, D3Q19_C)]), _MM_HINT_T0);
 732                                 }
 733                         #endif
 734
 735                         adjListIndex = (index - (index % VSIZE)) * N_D3Q19_IDX + (index % VSIZE);
 736                         #define X(name, idx, idxinv, _x, _y, _z)        JOIN(vgatheridx_,name) = VLIU(&(ADJ_LIST(idxinv)));
 737                                 D3Q19_LIST_WO_C
 738                         #undef X
 739
 740                         #define X(name, idx, idxinv, _x, _y, _z)        JOIN(vpdf_,name) = VG32(JOIN(vgatheridx_,name), src, 8);
 741                                 D3Q19_LIST_WO_C
 742                         #undef X
 743
 744                         ppdf_C = &(src[P_INDEX_3(nCells, index, D3Q19_C)]);
 745                         vpdf_C = VLDU(ppdf_C);
 746
 747                         // macroscopic values
 748                         //vux = vpdf_E + vpdf_NE + vpdf_SE + vpdf_TE + vpdf_BE -
 749                         //      vpdf_W - vpdf_NW - vpdf_SW - vpdf_TW - vpdf_BW;
 750                         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);
 751                         //vuy = vpdf_N + vpdf_NE + vpdf_NW + vpdf_TN + vpdf_BN -
 752                         //      vpdf_S - vpdf_SE - vpdf_SW - vpdf_TS - vpdf_BS;
 753                         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);
 754                         //vuz = vpdf_T + vpdf_TE + vpdf_TW + vpdf_TN + vpdf_TS -
 755                         //      vpdf_B - vpdf_BE - vpdf_BW - vpdf_BN - vpdf_BS;
 756                         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);
 757
 758                         //vdens = vpdf_C +
 759                         //        vpdf_N  + vpdf_E  + vpdf_S  + vpdf_W  +
 760                         //        vpdf_NE + vpdf_SE + vpdf_SW + vpdf_NW +
 761                         //        vpdf_T  + vpdf_TN + vpdf_TE + vpdf_TS + vpdf_TW +
 762                         //        vpdf_B  + vpdf_BN + vpdf_BE + vpdf_BS + vpdf_BW;
 763                         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)),
 764                                                                                 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));
 765
 766                         //vdir_indep_trm = vdens - (vux * vux + vuy * vuy + vuz * vuz) * VTHREE_HALF;
 767                         vdir_indep_trm = VSUB(vdens,VMUL(VADD(VADD(VMUL(vux,vux),VMUL(vuy,vuy)),VMUL(vuz,vuz)),VTHREE_HALF));
 768
 769                         // collide direction: w_0
 770                         //src[I(index, D3Q19_C)]  =[UA] vpdf_C - vomegaEven * (vpdf_C - vw_0 * vdir_indep_trm);
 771                         VSTU(ppdf_C,VSUB(vpdf_C,VMUL(vomegaEven,VSUB(vpdf_C,VMUL(vw_0,vdir_indep_trm)))));
 772
 773
 774                         // collide axis aligend pdfs vectorized
 775                         #define SCAT(offsets, vsrc) VS32(src, offsets, vsrc, 8)
 776
 777                         // collide direction: w_1
 778                         vw_1_indep = VMUL(vw_1,vdir_indep_trm);
 779
 780                         #define COLLIDE_AA_V(tmpVui, dir1, dir2) \
 781                                 vui = tmpVui; \
 782                                 vevenPart = VMUL(vomegaEven,VSUB(VSUB(VMUL(VONE_HALF,VADD(JOIN(vpdf_,dir1),JOIN(vpdf_,dir2))),VMUL(vui,VMUL(vui,vw_1_nine_half))),vw_1_indep));\
 783                                 voddPart = VMUL(vomegaOdd,VSUB(VMUL(VONE_HALF,VSUB(JOIN(vpdf_,dir1),JOIN(vpdf_,dir2))),VMUL(vui,vw_1_x3)));\
 784                                 SCAT(JOIN(vgatheridx_,dir2),VSUB(VSUB(JOIN(vpdf_,dir1),vevenPart),voddPart));\
 785                                 SCAT(JOIN(vgatheridx_,dir1),VADD(VSUB(JOIN(vpdf_,dir2),vevenPart),voddPart));
 786
 787                         COLLIDE_AA_V(vuy, N, S)
 788                         COLLIDE_AA_V(vux, E, W)
 789                         COLLIDE_AA_V(vuz, T, B)
 790
 791                         #undef COLLIDE_AA_V
 792
 793                         // collide direction: w_2
 794                         vw_2_indep = VMUL(vw_2,vdir_indep_trm);
 795
 796                         #define COLLIDE_UA_V(tmpVui, dir1, dir2) \
 797                                 vui = tmpVui; \
 798                                 vevenPart = VMUL(vomegaEven,VSUB(VSUB(VMUL(VONE_HALF,VADD(JOIN(vpdf_,dir1),JOIN(vpdf_,dir2))),VMUL(vui,VMUL(vui,vw_2_nine_half))),vw_2_indep));\
 799                                 voddPart = VMUL(vomegaOdd,VSUB(VMUL(VONE_HALF,VSUB(JOIN(vpdf_,dir1),JOIN(vpdf_,dir2))),VMUL(vui,vw_2_x3)));\
 800                                 SCAT(JOIN(vgatheridx_,dir2),VSUB(VSUB(JOIN(vpdf_,dir1),vevenPart),voddPart)); \
 801                                 SCAT(JOIN(vgatheridx_,dir1),VADD(VSUB(JOIN(vpdf_,dir2),vevenPart),voddPart));
 802
 803                         COLLIDE_UA_V(VSUB(vuy,vux), NW, SE)
 804                         COLLIDE_UA_V(VADD(vux,vuy), NE, SW)
 805                         COLLIDE_UA_V(VSUB(vuz,vux), TW, BE)
 806                         COLLIDE_UA_V(VADD(vux,vuz), TE, BW)
 807                         COLLIDE_UA_V(VSUB(vuz,vuy), TS, BN)
 808                         COLLIDE_UA_V(VADD(vuy,vuz), TN, BS)
 809
 810                         #undef COLLIDE_UA_V
 811                         #undef SCAT
 812                 } // gather & scatter loop
 813
 814         } // Vectorized Loop
 815
 816         // Scalar remainder loop
 817         #pragma novector
 818         #pragma loop_count max=7
 819         for (int index = loopStartIndices[consecStopIndex]; index < indexStop; ++index) {
 820
 821                 adjListIndex = (index - (index % VSIZE)) * N_D3Q19_IDX + (index % VSIZE);
 822                 #define X(name, idx, idxinv, _x, _y, _z)        JOIN(ppdf_,name) = &(src[ADJ_LIST(idxinv)]);
 823                         D3Q19_LIST_WO_C
 824                 #undef X
 825                 #define X(name, idx, idxinv, _x, _y, _z)        JOIN(pdf_,name) = *(JOIN(ppdf_,name));
 826                         D3Q19_LIST_WO_C
 827                 #undef X
 828
 829                 ppdf_C = &(src[P_INDEX_3(nCells, index, D3Q19_C)]);
 830                 pdf_C = *ppdf_C;
 831
 832                 // macroscopic values
 833                 ux = pdf_E + pdf_NE + pdf_SE + pdf_TE + pdf_BE -
 834                         pdf_W - pdf_NW - pdf_SW - pdf_TW - pdf_BW;
 835                 uy = pdf_N + pdf_NE + pdf_NW + pdf_TN + pdf_BN -
 836                         pdf_S - pdf_SE - pdf_SW - pdf_TS - pdf_BS;
 837                 uz = pdf_T + pdf_TE + pdf_TW + pdf_TN + pdf_TS -
 838                         pdf_B - pdf_BE - pdf_BW - pdf_BN - pdf_BS;
 839
 840                 dens = pdf_C +
 841                         pdf_N  + pdf_E  + pdf_S  + pdf_W  +
 842                         pdf_NE + pdf_SE + pdf_SW + pdf_NW +
 843                         pdf_T  + pdf_TN + pdf_TE + pdf_TS + pdf_TW +
 844                         pdf_B  + pdf_BN + pdf_BE + pdf_BS + pdf_BW;
 845
 846                 dir_indep_trm = dens - (ux * ux + uy * uy + uz * uz)*3.0/2.0;
 847
 848                 // collide direction: w_0
 849                 *ppdf_C  = pdf_C - omegaEven * (pdf_C - w_0 * dir_indep_trm);
 850
 851                 // collide direction: w_1
 852                 w_1_indep = w_1 * dir_indep_trm;
 853
 854                 #define COLLIDE_AA_S(tmpUi, dir1, dir2) \
 855                         ui = tmpUi; \
 856                         evenPart = omegaEven * (0.5 * (JOIN(pdf_,dir1) + JOIN(pdf_,dir2)) - ui * ui * w_1_nine_half - w_1_indep); \
 857                         oddPart  = omegaOdd  * (0.5 * (JOIN(pdf_,dir1) - JOIN(pdf_,dir2)) - ui * w_1_x3); \
 858                         *(JOIN(ppdf_,dir2))  = JOIN(pdf_,dir1) - evenPart - oddPart; \
 859                         *(JOIN(ppdf_,dir1))  = JOIN(pdf_,dir2) - evenPart + oddPart;
 860
 861                 COLLIDE_AA_S(uy, N, S)
 862                 COLLIDE_AA_S(ux, E, W)
 863                 COLLIDE_AA_S(uz, T, B)
 864
 865                 #undef COLLIDE_AA_S
 866
 867                 // collide direction: w_2
 868                 w_2_indep = w_2 * dir_indep_trm;
 869
 870                 #define COLLIDE_UA_S(tmpUi, dir1, dir2) \
 871                         ui = tmpUi; \
 872                         evenPart = omegaEven * (0.5 * (JOIN(pdf_,dir1) + JOIN(pdf_,dir2)) - ui * ui * w_2_nine_half - w_2_indep); \
 873                         oddPart  = omegaOdd  * (0.5 * (JOIN(pdf_,dir1) - JOIN(pdf_,dir2)) - ui * w_2_x3); \
 874                         *(JOIN(ppdf_,dir2)) = JOIN(pdf_,dir1) - evenPart - oddPart; \
 875                         *(JOIN(ppdf_,dir1)) = JOIN(pdf_,dir2) - evenPart + oddPart;
 876
 877                 COLLIDE_UA_S((-ux + uy), NW, SE)
 878                 COLLIDE_UA_S(( ux + uy), NE, SW)
 879                 COLLIDE_UA_S((-ux + uz), TW, BE)
 880                 COLLIDE_UA_S(( ux + uz), TE, BW)
 881                 COLLIDE_UA_S((-uy + uz), TS, BN)
 882                 COLLIDE_UA_S(( uy + uz), TN, BS)
 883
 884                 #undef COLLIDE_UA_S
 885
 886         } // remainder loop over fluid nodes
 887
 888 #undef SOFTWARE_PREFETCH_LOAD_STORE
 889 #undef SOFTWARE_PREFETCH_LOAD_STORE_HINT
 890 #undef SOFTWARE_PREFETCH_LOAD_STORE_LOOKAHEAD
 891 #undef SOFTWARE_PREFETCH_GATHER_SCATTER
 892 #undef SOFTWARE_PREFETCH_GATHER_SCATTER_HINT
 893 #undef SOFTWARE_PREFETCH_GATHER_SCATTER_LOOKAHEAD
 894 #undef ADJ_LIST
 895 #undef I
 896 }