--- /dev/null
+// --------------------------------------------------------------------------
+//
+// Copyright
+// Markus Wittmann, 2016-2017
+// RRZE, University of Erlangen-Nuremberg, Germany
+// markus.wittmann -at- fau.de or hpc -at- rrze.fau.de
+//
+// Viktor Haag, 2016
+// LSS, University of Erlangen-Nuremberg, Germany
+//
+// Michael Hussnaetter, 2017-2018
+// University of Erlangen-Nuremberg, Germany
+// michael.hussnaetter -at- fau.de
+//
+// This file is part of the Lattice Boltzmann Benchmark Kernels (LbmBenchKernels).
+//
+// LbmBenchKernels is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// LbmBenchKernels is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with LbmBenchKernels. If not, see <http://www.gnu.org/licenses/>.
+//
+// --------------------------------------------------------------------------
+#include "BenchKernelD3Q19ListAaPvGatherHybridCommon.h"
+
+#include "Memory.h"
+#include "Vtk.h"
+#include "Vector.h"
+
+#include <math.h>
+
+#ifdef _OPENMP
+#include <omp.h>
+#endif
+
+#define PAGE_4K 4096
+
+#if ALLOC_ADJ_LIST_IN_HBM == 1
+#define ADJ_LIST_ALLOCATOR HbwAllocAligned
+#define ADJ_LIST_FREE HbwFree
+#else
+#define ADJ_LIST_ALLOCATOR MemAllocAligned
+#define ADJ_LIST_FREE MemFree
+#endif
+
+#if ALLOC_PDF_IN_HBM == 1
+#define PDF_ALLOCATOR HbwAllocAligned
+#define PDF_FREE HbwFree
+#else
+#define PDF_ALLOCATOR MemAllocAligned
+#define PDF_FREE MemFree
+#endif
+
+// Forward definition.
+void FNAME(D3Q19ListAaPvGatherHybridKernel)(LatticeDesc * ld, struct KernelData_ * kd, CaseData * cd);
+
+
+
+// -----------------------------------------------------------------------
+// Functions which are used as callback by the kernel to read or write
+// PDFs and nodes.
+
+static void FNAME(BCGetPdf)(KernelData * kd, int x, int y, int z, int dir, PdfT * pdf)
+{
+ Assert(kd != NULL);
+ Assert(kd->PdfsActive != NULL);
+ Assert(kd->PdfsActive == kd->Pdfs[0] || kd->PdfsActive == kd->Pdfs[1]);
+ Assert(pdf != NULL);
+
+ Assert(x >= 0); Assert(y >= 0); Assert(z >= 0);
+ Assert(x < kd->Dims[0]); Assert(y < kd->Dims[1]); Assert(z < kd->Dims[2]);
+ Assert(dir >= 0); Assert(dir < N_D3Q19);
+
+ KernelDataList * kdl = KDL(kd);
+ uint32_t * adjList = kdl->AdjList;
+
+ if (kdl->Iteration % 2 == 0) {
+ // Pdfs are stored inverse, local PDFs are located in remote nodes
+
+ // getting node index
+ uint32_t index = kdl->Grid[L_INDEX_4(kd->Dims, x, y, z)];
+
+ if (dir != D3Q19_C) {
+ #define ADJ_LIST(dir) adjList[(index - (index % VSIZE)) * N_D3Q19_IDX + (dir * VSIZE) + (index % VSIZE)]
+ *pdf = kd->PdfsActive[ADJ_LIST(D3Q19_INV[dir])];
+ #undef ADJ_LIST
+ }
+ else {
+ *pdf = kd->PdfsActive[P_INDEX_3(kdl->nCells, index, dir)];
+ }
+
+ }
+ else {
+ *pdf = kd->PdfsActive[P_INDEX_5(kdl, x, y, z, dir)];
+ }
+
+ return;
+}
+
+static void FNAME(BCSetPdf)(KernelData * kd, int x, int y, int z, int dir, PdfT pdf)
+{
+ Assert(kd != NULL);
+ Assert(kd->PdfsActive != NULL);
+ Assert(kd->PdfsActive == kd->Pdfs[0] || kd->PdfsActive == kd->Pdfs[1]);
+ Assert(x >= 0); Assert(y >= 0); Assert(z >= 0);
+ Assert(x < kd->Dims[0]); Assert(y < kd->Dims[1]); Assert(z < kd->Dims[2]);
+ Assert(dir >= 0); Assert(dir < N_D3Q19);
+
+ if (isnan(pdf)) {
+ printf("ERROR: setting nan %d %d %d %d %s\n", x, y, z, dir, D3Q19_NAMES[dir]);
+ DEBUG_BREAK_POINT();
+ exit(1);
+ }
+
+ KernelDataList * kdl = KDL(kd);
+ uint32_t * adjList = kdl->AdjList;
+
+ if (kdl->Iteration % 2 == 0) {
+ // Pdfs are stored inverse, local PDFs are located in remote nodes
+
+ // getting node index
+ uint32_t index = kdl->Grid[L_INDEX_4(kd->Dims, x, y, z)];
+
+ if (dir != D3Q19_C) {
+ #define ADJ_LIST(dir) adjList[(index - (index % VSIZE)) * N_D3Q19_IDX + (dir * VSIZE) + (index % VSIZE)]
+ kd->PdfsActive[ADJ_LIST(D3Q19_INV[dir])] = pdf;
+ #undef ADJ_LIST
+ } else {
+ kd->PdfsActive[P_INDEX_3(kdl->nCells, index, dir)] = pdf;
+ }
+
+ } else {
+ kd->PdfsActive[P_INDEX_5(kdl, x, y, z, dir)] = pdf;
+ }
+
+ return;
+}
+
+
+static void GetNode(KernelData * kd, int x, int y, int z, PdfT * pdfs)
+{
+ Assert(kd != NULL);
+ Assert(kd->PdfsActive != NULL);
+ Assert(kd->PdfsActive == kd->Pdfs[0] || kd->PdfsActive == kd->Pdfs[1]);
+ Assert(pdfs != NULL);
+ Assert(x >= 0); Assert(y >= 0); Assert(z >= 0);
+ Assert(x < kd->Dims[0]); Assert(y < kd->Dims[1]); Assert(z < kd->Dims[2]);
+
+ KernelDataList * kdl = KDL(kd);
+ uint32_t * adjList = kdl->AdjList;
+
+ if(kdl->Iteration % 2 == 0){
+
+ uint32_t index = kdl->Grid[L_INDEX_4(kdl->kd.Dims, x, y, z)];
+
+ // Load PDFs of local cell: pdf_N = src[adjList[adjListIndex + D3Q19_S]]; ...
+ #define ADJ_LIST(dir) adjList[(index - (index % VSIZE)) * N_D3Q19_IDX + (dir * VSIZE) + (index % VSIZE)]
+ #define X(name, idx, idxinv, _x, _y, _z) pdfs[idx] = kd->PdfsActive[ADJ_LIST(idxinv)];
+ D3Q19_LIST_WO_C
+ #undef X
+ #undef ADJ_LIST
+ pdfs[D3Q19_C] = kd->PdfsActive[P_INDEX_3(kdl->nCells, index, D3Q19_C)];
+
+ } else {
+
+ #define I(x, y, z, dir) P_INDEX_5(KDL(kd), (x), (y), (z), (dir))
+ #define X(name, idx, idxinv, _x, _y, _z) pdfs[idx] = kd->PdfsActive[I(x, y, z, idx)];
+ D3Q19_LIST
+ #undef X
+ #undef I
+
+ }
+
+ for (int d = 0; d < 19; ++d) {
+ if(isnan(pdfs[d]) || isinf(pdfs[d])) {
+ printf("%d %d %d %d nan! get node\n", x, y, z, d);
+ for (int d2 = 0; d2 < 19; ++d2) {
+ printf("%d: %e\n", d2, pdfs[d2]);
+ }
+ exit(1);
+ }
+ }
+
+ return;
+}
+
+
+static void SetNode(KernelData * kd, int x, int y, int z, PdfT * pdfs)
+{
+ Assert(kd != NULL);
+ Assert(kd->PdfsActive != NULL);
+ Assert(kd->PdfsActive == kd->Pdfs[0] || kd->PdfsActive == kd->Pdfs[1]);
+ Assert(pdfs != NULL);
+
+ Assert(x >= 0); Assert(y >= 0); Assert(z >= 0);
+ Assert(x < kd->Dims[0]); Assert(y < kd->Dims[1]); Assert(z < kd->Dims[2]);
+
+ for (int d = 0; d < 19; ++d) {
+ if(isnan(pdfs[d])) {
+ printf("%d %d %d %d nan! get node\n", x, y, z, d);
+ for (int d2 = 0; d2 < 19; ++d2) {
+ printf("%d: %e\n", d2, pdfs[d2]);
+ }
+ exit(1);
+ }
+ }
+
+ KernelDataList * kdl = KDL(kd);
+ uint32_t * adjList = kdl->AdjList;
+
+ if(kdl->Iteration % 2 == 0){
+
+ uint32_t index = kdl->Grid[L_INDEX_4(kdl->kd.Dims, x, y, z)];
+
+ // Load PDFs of local cell: pdf_N = src[adjList[adjListIndex + D3Q19_S]]; ...
+ kd->PdfsActive[P_INDEX_3(kdl->nCells, index, D3Q19_C)] = pdfs[D3Q19_C];
+
+ #define ADJ_LIST(dir) adjList[(index - (index % VSIZE)) * N_D3Q19_IDX + (dir * VSIZE) + (index % VSIZE)]
+ #define X(name, idx, idxinv, _x, _y, _z) kd->PdfsActive[ADJ_LIST(idxinv)] = pdfs[idx];
+ D3Q19_LIST_WO_C
+ #undef X
+ #undef ADJ_LIST
+
+ } else {
+
+ #define I(x, y, z, dir) P_INDEX_5(KDL(kd), (x), (y), (z), (dir))
+ #define X(name, idx, idxinv, _x, _y, _z) kd->PdfsActive[I(x, y, z, idx)] = pdfs[idx];
+ D3Q19_LIST
+ #undef X
+ #undef I
+
+ }
+
+ return;
+}
+
+static void ParameterUsage()
+{
+ printf("Kernel parameters:\n");
+ printf(" [-blk <n>] [-blk-[xyz] <n>]\n");
+
+ return;
+}
+
+static void ParseParameters(Parameters * params, int * blk)
+{
+ Assert(blk != NULL);
+
+ blk[0] = 0; blk[1] = 0; blk[2] = 0;
+
+ #define ARG_IS(param) (!strcmp(params->KernelArgs[i], param))
+ #define NEXT_ARG_PRESENT() \
+ do { \
+ if (i + 1 >= params->nKernelArgs) { \
+ printf("ERROR: argument %s requires a parameter.\n", params->KernelArgs[i]); \
+ exit(1); \
+ } \
+ } while (0)
+
+
+ for (int i = 0; i < params->nKernelArgs; ++i) {
+ if (ARG_IS("-blk") || ARG_IS("--blk")) {
+ NEXT_ARG_PRESENT();
+
+ int tmp = strtol(params->KernelArgs[++i], NULL, 0);
+
+ if (tmp <= 0) {
+ printf("ERROR: blocking parameter must be > 0.\n");
+ exit(1);
+ }
+
+ blk[0] = blk[1] = blk[2] = tmp;
+ }
+ else if (ARG_IS("-blk-x") || ARG_IS("--blk-x")) {
+ NEXT_ARG_PRESENT();
+
+ int tmp = strtol(params->KernelArgs[++i], NULL, 0);
+
+ if (tmp <= 0) {
+ printf("ERROR: blocking parameter must be > 0.\n");
+ exit(1);
+ }
+
+ blk[0] = tmp;
+ }
+ else if (ARG_IS("-blk-y") || ARG_IS("--blk-y")) {
+ NEXT_ARG_PRESENT();
+
+ int tmp = strtol(params->KernelArgs[++i], NULL, 0);
+
+ if (tmp <= 0) {
+ printf("ERROR: blocking parameter must be > 0.\n");
+ exit(1);
+ }
+
+ blk[1] = tmp;
+ }
+ else if (ARG_IS("-blk-z") || ARG_IS("--blk-z")) {
+ NEXT_ARG_PRESENT();
+
+ int tmp = strtol(params->KernelArgs[++i], NULL, 0);
+
+ if (tmp <= 0) {
+ printf("ERROR: blocking parameter must be > 0.\n");
+ exit(1);
+ }
+
+ blk[2] = tmp;
+ }
+ else if (ARG_IS("-h") || ARG_IS("-help") || ARG_IS("--help")) {
+ ParameterUsage();
+ exit(1);
+ }
+ else {
+ printf("ERROR: unknown kernel parameter.\n");
+ ParameterUsage();
+ exit(1);
+ }
+ }
+
+ #undef ARG_IS
+ #undef NEXT_ARG_PRESENT
+
+ return;
+}
+
+static void SetuploopStartIndices(LatticeDesc * ld, KernelDataListRia * kdlr, int nThreads)
+{
+ //#define ADJ_LIST(dir) adjList[(index - (index % VSIZE)) * N_D3Q19_IDX + (dir * VSIZE) + (index % VSIZE)]
+ Assert(ld != NULL);
+ Assert(kdlr != NULL);
+ Assert(nThreads > 0);
+
+ //uint32_t * adjList = kdlr->kdl.AdjList;
+ uint32_t * adjList = kdlr->kdl.AdjList;
+
+ uint32_t nLoopStartIndices = 0;
+ uint32_t loopStartIndex = 2;
+
+ int nFluid = kdlr->kdl.nFluid;
+
+ int * oddKernelThreadStartIndices = (int *)malloc(sizeof(int) * (nThreads + 1));
+
+ int nNodesPerThread = nFluid / nThreads;
+ //printf("nodesPerThread: %d\n", nNodesPerThread);
+
+ for (int i = 0; i < nThreads; ++i) {
+ oddKernelThreadStartIndices[i] = i * nNodesPerThread + MinI(i, nFluid % nThreads);
+ }
+
+ oddKernelThreadStartIndices[nThreads] = nFluid;
+
+ /*
+ for (int i = 0; i <= nThreads; ++i) {
+ printf("oddKernelThreadStartIndices[%d] = %d\n", i, oddKernelThreadStartIndices[i]);
+ }
+ */
+
+ int threadIndex = 1;
+
+ // We execute following code two times.
+ // - The first time to get the count of how many entries we need for the
+ // loopStartIndices array.
+ // - The second time to fill the array.
+
+ // Loop over adjacency list of all nodes.
+ // Compare if adjacent nodes within one cache line share the same access pattern.
+ // First vectorized access is assumed to be consecutive (-> may be loaded with regular load).
+
+ int lastCacheLineConsecutive = 1;
+
+ for (int fluidBaseIndex = 1; fluidBaseIndex < nFluid + 1; fluidBaseIndex += VSIZE) {
+
+ int currentCacheLineConsecutive = 1;
+
+ // Loop over all directions except the center one.
+ for(int d = 0; d < N_D3Q19 - 1; ++d) {
+ Assert(d != D3Q19_C);
+
+ // check if cache line itself has consecutive memory access pattern
+ for(int inChunkIndex = 0; (inChunkIndex < VSIZE - 1) && ((fluidBaseIndex + inChunkIndex) < nFluid); ++inChunkIndex) {
+ int index = fluidBaseIndex + inChunkIndex;
+
+ Assert(index < nFluid);
+
+ #define ADJ_LIST(idx, dir) adjList[((idx) - ((idx) % VSIZE)) * N_D3Q19_IDX + ((dir) * VSIZE) + ((idx) % VSIZE)]
+ //if (adjList[index * N_D3Q19_IDX + d] != adjList[(index - 1) * N_D3Q19_IDX + d] + 1)
+ if (ADJ_LIST(index, d) != ADJ_LIST(index-1, d) + 1) {
+ //printf("no match for index: %d\n", d);
+ //printf("ADJ_LlST(%d,%d) = %d != %d = ADJ_LlST(%d,%d) + 1\n", index, d, ADJ_LIST(index,d), ADJ_LIST(index-1,d), index-1, d);
+ // Different access pattern.
+ currentCacheLineConsecutive = 0;
+ break;
+ }
+ #undef ADJ_LIST
+
+ }
+
+ if(!currentCacheLineConsecutive){
+ break; //exit from nested loop
+ }
+
+ }
+
+ int interCacheLineConsecutive = 1;
+
+ if(currentCacheLineConsecutive && lastCacheLineConsecutive){
+ // check if cache line has consecutive memory access pattern to last entry of previous cache line
+ int lastIdxOfPreviousCacheLine = fluidBaseIndex - 2;
+ if (lastIdxOfPreviousCacheLine > 0) {
+ for(int d = 0; d < N_D3Q19 - 1; ++d) {
+ Assert(d != D3Q19_C);
+ #define ADJ_LIST(idx, dir) adjList[((idx) - ((idx) % VSIZE)) * N_D3Q19_IDX + ((dir) * VSIZE) + ((idx) % VSIZE)]
+ if (ADJ_LIST(fluidBaseIndex-1, d) != ADJ_LIST(lastIdxOfPreviousCacheLine, d) + 1) {
+ // Different access pattern.
+ //printf("not interCacheConsecutive\n");
+ interCacheLineConsecutive = 0;
+ break;
+ }
+ #undef ADJ_LIST
+
+ }
+ }
+ }
+ int threadBoundaryIndex = oddKernelThreadStartIndices[threadIndex];
+ if (fluidBaseIndex - 1 <= threadBoundaryIndex &&
+ threadBoundaryIndex < fluidBaseIndex + VSIZE - 1) {
+ // Current cache line contains thread boundary.
+ // These cache lines are treated by scalar peel and
+ // reminder loops in kernel of every thread.
+ // TODO maybe replace these loops with masked gather / scatter?!
+ if (loopStartIndex % 2 == 0) { // current index would be gather/scatter index
+ ++loopStartIndex; // reserving gather/scatter index
+ }
+ ++loopStartIndex; // reserving space for start remainder loop of thread n
+
+ if (threadIndex < nThreads){
+ ++loopStartIndex; // reserving space for starting peel loop of thread n+1
+
+ if (fluidBaseIndex - 1 == threadBoundaryIndex){
+ if(!currentCacheLineConsecutive){
+ ++loopStartIndex;
+ }
+ }
+ else {
+ currentCacheLineConsecutive = 1;
+ }
+
+ //++loopStartIndex; // reserving space for ending peel loop / starting load/store of thread n+1
+ ++loopStartIndex; // 1st load/store end == 1st gather/scatter start OR 1st gather/scatter end == 2nd load/start start
+ }
+ ++threadIndex;
+ }
+ else {
+ // We are not at a thread boundary.
+ if (currentCacheLineConsecutive) {
+ if(lastCacheLineConsecutive && !interCacheLineConsecutive){
+ loopStartIndex+=2;
+ }
+ else if(!lastCacheLineConsecutive){
+ ++loopStartIndex;
+ }
+ }
+ else {
+ if(lastCacheLineConsecutive){
+ ++loopStartIndex;
+ }
+ }
+ }
+
+ // treating special case when last thread has no remainder loop
+ if (oddKernelThreadStartIndices[nThreads] == fluidBaseIndex + VSIZE - 1) {
+ //printf("--> special case 111. loopStartIndex: %d \n", loopStartIndex);
+ if (loopStartIndex % 2 != 0) { // current index is gather/scatter end and load/store start index
+ ++loopStartIndex; //set load/store end (gather/scatter start) to same value as scalar remainder start => no more access to gather/scatter loop
+ }
+
+ ++loopStartIndex; // gather/scatter end and scalar remainder start
+ ++loopStartIndex; // scalar remainder end and scalar peel start
+
+ }
+
+ lastCacheLineConsecutive = currentCacheLineConsecutive;
+ }
+
+ if (nFluid > 0) {
+ nLoopStartIndices = loopStartIndex;
+ }
+
+ int * loopStartIndices;
+ unsigned long loopStartIndicesByte = (nLoopStartIndices + 1) * sizeof(int);
+
+ printf("# Loop Start Index Array Allocation:\n");
+ printf("# elements: \t\t%d\n", nLoopStartIndices + 1);
+ printf("# size: \t\t%e MiB\n", loopStartIndicesByte / 1024.0 / 1024.0);
+ printf("# alignment: \t\t%d b\n", PAGE_4K);
+
+ if (MemAllocAligned((void **)&loopStartIndices, loopStartIndicesByte, PAGE_4K)) {
+ printf("ERROR: allocating loopStartIndices array with MemAllocAligned failed: %lu bytes.\n", loopStartIndicesByte);
+ exit(1);
+ }
+ else {
+ printf("# allocator: \t\t\tMemAllocAligned()\n");
+ }
+
+ oddKernelThreadStartIndices[0] = 0;
+ loopStartIndices[0] = 0; //first scalar loop would start with 0
+ loopStartIndices[1] = 0; //no peel loop expected -> first load/store loop may start at index==0
+ loopStartIndices[2] = 0; //may not be set in case first access is gather/scatter -> therefore its set here
+
+ // resetting values to default
+ threadIndex = 1;
+ lastCacheLineConsecutive = 1;
+ loopStartIndex = 2;
+
+ // Loop over adjacency list of all nodes.
+ // Compare if adjacent nodes share the same access pattern.
+
+ int indexAccumulator = 0;
+
+ // for statistical reasons:
+ int gatherAccumulator = 0;
+ int loadAccumulator = 0;
+ int scalarLookups = 0;
+ int loadLookups = 0;
+
+
+ for (int fluidBaseIndex = 1; fluidBaseIndex < nFluid + 1; fluidBaseIndex += VSIZE) {
+ int currentCacheLineConsecutive = 1;
+ //printf("fluidbaseIndex: %d\n", fluidBaseIndex);
+ // Loop over all directions except the center one.
+ for(int d = 0; d < N_D3Q19 - 1; ++d) {
+ Assert(d != D3Q19_C);
+
+ // check if cache line itself has consecutive memory access pattern
+ for(int inChunkIndex = 0; (inChunkIndex < VSIZE - 1) && ((fluidBaseIndex + inChunkIndex) < nFluid); ++inChunkIndex){
+ int index = fluidBaseIndex + inChunkIndex;
+
+ Assert(index < nFluid);
+
+ #define ADJ_LIST(idx, dir) adjList[((idx) - ((idx) % VSIZE)) * N_D3Q19_IDX + ((dir) * VSIZE) + ((idx) % VSIZE)]
+ //if (adjList[index * N_D3Q19_IDX + d] != adjList[(index - 1) * N_D3Q19_IDX + d] + 1)
+ if (ADJ_LIST(index, d) != ADJ_LIST(index-1, d) + 1) {
+ // Different access pattern.
+ currentCacheLineConsecutive = 0;
+ break;
+ }
+ #undef ADJ_LIST
+ }
+
+ if(!currentCacheLineConsecutive){
+ break; //exit from nested loop
+ }
+ }
+
+ int interCacheLineConsecutive = 1;
+
+ if(currentCacheLineConsecutive && lastCacheLineConsecutive){
+ // check if cache line has consecutive memory access pattern to last entry of previous cache line
+ int lastIdxOfPreviousCacheLine = fluidBaseIndex - 2;
+ if (lastIdxOfPreviousCacheLine > 0) {
+ for(int d = 0; d < N_D3Q19 - 1; ++d) {
+ Assert(d != D3Q19_C);
+ #define ADJ_LIST(idx, dir) adjList[((idx) - ((idx) % VSIZE)) * N_D3Q19_IDX + ((dir) * VSIZE) + ((idx) % VSIZE)]
+ if (ADJ_LIST(fluidBaseIndex-1, d) != ADJ_LIST(lastIdxOfPreviousCacheLine, d) + 1) {
+ // Different access pattern.
+ interCacheLineConsecutive = 0;
+ break;
+ }
+ #undef ADJ_LIST
+
+ }
+ }
+ }
+
+ int threadBoundaryIndex = oddKernelThreadStartIndices[threadIndex];
+ //if (fluidBaseIndex > 3500)
+ // printf("threadBoundaryIndex: %d fluidBaseIndex-1: %d fluidBaseIndex + VSIZE - 1: %d\n", threadBoundaryIndex, fluidBaseIndex-1, fluidBaseIndex + VSIZE -1);
+
+ if (fluidBaseIndex - 1 <= threadBoundaryIndex &&
+ threadBoundaryIndex < fluidBaseIndex + VSIZE - 1) {
+ // Current cache line contains thread boundary.
+ // These cache lines are treated by scalar peel and
+ // reminder loops in kernel of every thread.
+ // TODO maybe replace these loops with masked gather / scatter?!
+ if (loopStartIndex % 2 == 0) { // current index would be gather/scatter index
+ //loopStartIndices[loopStartIndex] = fluidBaseIndex - 1; //same value as scalar remainder start => no more access to gather/scatter loop
+ loopStartIndices[loopStartIndex] = indexAccumulator; //same value as scalar remainder start => no more access to gather/scatter loop
+ ++loopStartIndex;
+ }
+
+ //loopStartIndices[loopStartIndex] = fluidBaseIndex - 1; // gather/scatter end and scalar remainder start
+ loopStartIndices[loopStartIndex] = indexAccumulator; // gather/scatter end and scalar remainder start
+ ++loopStartIndex;
+
+ // starting indices of thread n+1
+ loopStartIndices[loopStartIndex] = threadBoundaryIndex; // scalar remainder of thread n end and scalar peel of thread n+1 start
+ oddKernelThreadStartIndices[threadIndex] = loopStartIndex; // thread start is where scalar peel starts
+
+ if (threadIndex < nThreads){
+ indexAccumulator = ((threadBoundaryIndex + VSIZE - 1) / VSIZE ) * VSIZE; // rounding towards next multiple of VSIZE
+ ++loopStartIndex;
+ loopStartIndices[loopStartIndex] = indexAccumulator; // scalar peel end and 1st load/store start
+
+ // treating special case when there is no peel / remainder loop
+ if (fluidBaseIndex - 1 == threadBoundaryIndex){
+ if(!currentCacheLineConsecutive){
+ ++loopStartIndex;
+ loopStartIndices[loopStartIndex] = indexAccumulator; // 1st load/store end and 1st gather/scatter start
+ gatherAccumulator += VSIZE;
+ }
+ else {
+ loadLookups += VSIZE;
+ }
+ indexAccumulator += VSIZE;
+ }
+ else {
+ scalarLookups += VSIZE;
+ currentCacheLineConsecutive = 1;
+ }
+
+ ++loopStartIndex; // 1st load/store end == 1st gather/scatter start OR 1st gather/scatter end == 2nd load/start start
+ loopStartIndices[loopStartIndex] = indexAccumulator; // 1st load/store end == 1st gather/scatter start OR 1st gather/scatter end == 2nd load/start start
+ }
+ ++threadIndex;
+
+ }
+ else {
+ // We are not at a thread boundary.
+ int print = 0;
+ if (currentCacheLineConsecutive) {
+ loadAccumulator += VSIZE;
+
+ if(lastCacheLineConsecutive && !interCacheLineConsecutive){
+ loadLookups += VSIZE;
+ if (print)
+ printf("#1 loopStartIndex: %d\n", loopStartIndex);
+ // loopStartIndices[loopStartIndex] is not incremented since pointers need to be fetched again.
+ // loopStartIndices[loopStartIndex + 1] (-> start Load/Store and end Gather/Scatter)
+ // gets same value as loopStartIndices[loopStartindex] (-> start Gather/Scatter)
+ // this ensures that no gather/scatter iteration is executed
+ ++loopStartIndex;
+ loopStartIndices[loopStartIndex] = indexAccumulator;
+
+ // loopStartIndices[loopStartIndex + 2] (-> start Gather/Scatter and end Load/Store)
+ // gets set to have one Load/Store iteration
+ ++loopStartIndex;
+ indexAccumulator+=VSIZE;
+ loopStartIndices[loopStartIndex] = indexAccumulator;
+
+ }
+ else if(!lastCacheLineConsecutive){
+ loadLookups += VSIZE;
+ if (print)
+ printf("#2 loopStartIndex: %d\n", loopStartIndex);
+ ++loopStartIndex;
+ indexAccumulator+=VSIZE;
+ loopStartIndices[loopStartIndex] = indexAccumulator;
+ }
+ else { // (lastCacheLineConsecutive && interCacheLineConsecutive)
+ if (print)
+ printf("#3 loopStartIndex: %d\n", loopStartIndex);
+ indexAccumulator+=VSIZE;
+ loopStartIndices[loopStartIndex] = indexAccumulator;
+ }
+ }
+ else {
+ gatherAccumulator += VSIZE;
+ if(lastCacheLineConsecutive){
+ if (print)
+ printf("#4 loopStartIndex: %d\n", loopStartIndex);
+ ++loopStartIndex;
+ indexAccumulator+=VSIZE;
+ loopStartIndices[loopStartIndex] = indexAccumulator;
+ }
+ else { // lastCacheLine without not consecutive memory access pattern
+ if (print)
+ printf("#5 loopStartIndex: %d\n", loopStartIndex);
+ indexAccumulator+=VSIZE;
+ loopStartIndices[loopStartIndex] = indexAccumulator;
+ }
+ }
+ }
+
+ // treating special case when last thread has no remainder loop
+ if (oddKernelThreadStartIndices[nThreads] == fluidBaseIndex + VSIZE - 1) {
+ //printf("--> special case. indexAccumulator: %d\n", indexAccumulator);
+ if (loopStartIndex % 2 != 0) { // current index is gather/scatter end and load/store start index
+ ++loopStartIndex;
+ loopStartIndices[loopStartIndex] = indexAccumulator; //set load/store end (gather/scatter start) to same value as scalar remainder start => no more access to gather/scatter loop
+ }
+
+ ++loopStartIndex;
+ loopStartIndices[loopStartIndex] = indexAccumulator; // gather/scatter end and scalar remainder start
+ ++loopStartIndex;
+ loopStartIndices[loopStartIndex] = indexAccumulator; // scalar remainder end and scalar peel start
+
+ oddKernelThreadStartIndices[threadIndex] = loopStartIndex; // thread start is where scalar peel starts
+ }
+
+ lastCacheLineConsecutive = currentCacheLineConsecutive;
+
+ }
+
+ if (nLoopStartIndices != loopStartIndex){
+ printf("ERROR: nLoopStartIndices unequal loopStartIndex!\n");
+ }
+
+ /*
+ printf("loopStartIndices:\n");
+ for(int i = 0; i <= nLoopStartIndices; ++i){
+ printf("%d ", loopStartIndices[i]);
+ }
+ printf("\n");
+ printf("oddKernelThreadStartIndices:\n");
+ for(int i = 0; i <= nThreads; ++i){
+ printf("%d ", oddKernelThreadStartIndices[i]);
+ }
+ printf("\n");
+ */
+
+ kdlr->loopStartIndices = loopStartIndices;
+ kdlr->nLoopStartIndices = nLoopStartIndices;
+
+ kdlr->oddKernelThreadStartIndices = oddKernelThreadStartIndices;
+ kdlr->nOddKernelThreadStartIndices = nThreads;
+
+ printf("# vload/vstore nodes: \t% 10d \t(%3.4f %% of total fluid nodes)\n", loadAccumulator, ((double) loadAccumulator / (double) nFluid) * 100);
+ printf("# gather/scatter nodes:\t% 10d \t(%3.4f %% of total fluid nodes)\n", gatherAccumulator, ((double) gatherAccumulator / (double) nFluid) * 100.0);
+ printf("# vload/vstore lookups:\t% 10d \n", loadLookups * (N_D3Q19 - 1));
+ printf("# gather/scatter lookups:\t% 10d \n", gatherAccumulator * (N_D3Q19 - 1));
+ printf("# scalar lookups: \t% 10d \n", scalarLookups * (N_D3Q19 - 1));
+
+ double loopBalanceEven = 2.0 * 19 * sizeof(PdfT);
+ double loopBalanceOdd = 2.0 * 19 * sizeof(PdfT) /* actual PDFs */
+ + (((double)(gatherAccumulator + loadLookups + scalarLookups)) / nFluid) * sizeof(int) * (N_D3Q19 - 1) /* AdjList */
+ + (nLoopStartIndices / nFluid) * sizeof(int); // one lookup to loopStartIndices
+
+ double loopBalance = (loopBalanceEven + loopBalanceOdd) / 2.0;
+
+ kdlr->kdl.kd.LoopBalance = loopBalance;
+
+ printf("# loop balance:\n");
+ printf("# even timestep: \t\t%.2f B/FLUP\n", loopBalanceEven);
+ printf("# odd timestep: \t\t%.2f B/FLUP\n", loopBalanceOdd);
+ printf("# average: \t\t%.2f B/FLUP\n", loopBalance);
+
+ return;
+}
+
+void FNAME(D3Q19ListAaPvGatherHybridInit)(LatticeDesc * ld, KernelData ** kernelData, Parameters * params)
+{
+ KernelData * kd;
+ KernelDataList * kdl;
+ KernelDataListRia * kdlr;
+ MemAlloc((void **)&kdlr, sizeof(KernelDataListRia));
+
+ kd = (KernelData *)kdlr;
+ kdl = KDL(kdlr);
+
+ *kernelData = kd;
+
+#ifdef DEBUG
+ kd->Pdfs[0] = NULL;
+ kd->Pdfs[1] = NULL;
+ kd->PdfsActive = NULL;
+ kd->DstPdfs = NULL;
+ kd->SrcPdfs = NULL;
+ kd->Dims[0] = -1;
+ kd->Dims[1] = -1;
+ kd->Dims[2] = -1;
+ kd->GlobalDims[0] = -1;
+ kd->GlobalDims[1] = -1;
+ kd->GlobalDims[2] = -1;
+ kd->Offsets[0] = -1;
+ kd->Offsets[1] = -1;
+ kd->Offsets[2] = -1;
+
+ kd->ObstIndices = NULL;
+ kd->nObstIndices = -1;
+ kd->BounceBackPdfsSrc = NULL;
+ kd->BounceBackPdfsDst = NULL;
+ kd->nBounceBackPdfs = -1;
+
+ kdl->AdjList = NULL;
+ kdl->Coords = NULL;
+ kdl->Grid = NULL;
+ kdl->nCells = -1;
+ kdl->nFluid = -1;
+
+ kdlr->loopStartIndices = NULL;
+ kdlr->nLoopStartIndices = 0;
+ kdlr->oddKernelThreadStartIndices = NULL;
+ kdlr->nOddKernelThreadStartIndices = 0;
+#endif
+
+ kdl->Iteration = -1;
+
+ // Ajust the dimensions according to padding, if used.
+ kd->Dims[0] = kd->GlobalDims[0] = ld->Dims[0];
+ kd->Dims[1] = kd->GlobalDims[1] = ld->Dims[1];
+ kd->Dims[2] = kd->GlobalDims[2] = ld->Dims[2];
+
+ int * lDims = ld->Dims;
+
+ int lX = lDims[0];
+ int lY = lDims[1];
+ int lZ = lDims[2];
+
+ int nTotalCells = lX * lY * lZ;
+ int nCells = ld->nFluid; // TODO: + padding
+ int nFluid = ld->nFluid;
+
+ // TODO: check nCells/nFluid do not exceed 2^31. This actually has to be
+ // done during lattice setup.
+ kdl->nCells = nCells;
+ kdl->nFluid = nFluid;
+
+ PdfT * pdfs[2];
+
+ int blk[3] = { 0 };
+
+ ParseParameters(params, blk);
+
+ if (blk[0] == 0) blk[0] = lX;
+ if (blk[1] == 0) blk[1] = lY;
+ if (blk[2] == 0) blk[2] = lZ;
+
+ printf("# blocking: \t\tx: %3d y: %3d z: %3d\n", blk[0], blk[1], blk[2]);
+
+ unsigned long latByte = nCells * sizeof(PdfT) * N_D3Q19;
+ unsigned long latFluidByte = nFluid * sizeof(PdfT) * N_D3Q19;
+ unsigned long latPadByte = (nCells - nFluid) * sizeof(PdfT) * N_D3Q19;
+
+ printf("# Lattice Array Allocation:\n");
+ printf("# lattice size: \t\t%e MiB\n", latByte / 1024.0 / 1024.0);
+ printf("# fluid lattice size:\t\t%e MiB\n", latFluidByte / 1024.0 / 1024.0);
+ printf("# lattice padding: \t\t%e MiB\n", latPadByte / 1024.0 / 1024.0);
+
+
+ printf("# alignment: \t\t%d b\n", PAGE_4K);
+
+ if (PDF_ALLOCATOR((void **)&pdfs[0], latFluidByte, PAGE_4K)) {
+ printf("ERROR: allocating PDF array with %s() failed: %lu bytes.\n", STRINGIFY(PDF_ALLOCATOR), latFluidByte);
+ exit(1);
+ }
+ else {
+ printf("# allocator: \t\t\t%s()\n", STRINGIFY(PDF_ALLOCATOR));
+ }
+
+ kd->Pdfs[0] = pdfs[0];
+
+ // Initialize PDFs with some (arbitrary) data for correct NUMA placement.
+ // Here we touch only the fluid nodes as this loop is OpenMP parallel and
+ // we want the same scheduling as in the kernel.
+ #ifdef _OPENMP
+ #pragma omp parallel for
+ #endif
+ for (int i = 0; i < nFluid; ++i) { for(int d = 0; d < N_D3Q19; ++d) {
+ pdfs[0][P_INDEX_3(nCells, i, d)] = 1.0;
+ } }
+
+ // Initialize all PDFs to some standard value.
+ for (int i = 0; i < nFluid; ++i) { for(int d = 0; d < N_D3Q19; ++d) {
+ pdfs[0][P_INDEX_3(nCells, i, d)] = 0.0;
+ } }
+
+ // ----------------------------------------------------------------------
+ // create grid which will hold the index numbers of the fluid nodes
+
+ uint32_t * grid;
+
+ if (MemAlloc((void **)&grid, nTotalCells * sizeof(uint32_t))) {
+ printf("ERROR: allocating grid for numbering failed: %lu bytes.\n", nTotalCells * sizeof(uint32_t));
+ exit(1);
+ }
+ kdl->Grid = grid;
+
+ int latticeIndex;
+
+#ifdef DEBUG
+ for(int z = 0; z < lZ; ++z) {
+ for(int y = 0; y < lY; ++y) {
+ for(int x = 0; x < lX; ++x) {
+
+ latticeIndex = L_INDEX_4(ld->Dims, x, y, z);
+
+ grid[latticeIndex] = ~0;
+ }
+ }
+ }
+#endif
+
+ // ----------------------------------------------------------------------
+ // generate numbering over grid
+
+ uint32_t * coords;
+
+ if (MemAlloc((void **)&coords, nFluid * sizeof(uint32_t) * 3)) {
+ printf("ERROR: allocating coords array failed: %lu bytes.\n", nFluid * sizeof(uint32_t) * 3);
+ exit(1);
+ }
+
+ kdl->Coords = coords;
+
+ // Index for the PDF nodes can start at 0 as we distinguish solid and fluid nodes
+ // through the ld->Lattice array.
+ int counter = 0;
+
+ // Blocking is implemented via setup of the adjacency list. The kernel later will
+ // walk through the lattice blocked automatically.
+ for (int bZ = 0; bZ < lZ; bZ += blk[2]) {
+ for (int bY = 0; bY < lY; bY += blk[1]) {
+ for (int bX = 0; bX < lX; bX += blk[0]) {
+
+ int eX = MIN(bX + blk[0], lX);
+ int eY = MIN(bY + blk[1], lY);
+ int eZ = MIN(bZ + blk[2], lZ);
+
+
+ for (int z = bZ; z < eZ; ++z) {
+ for (int y = bY; y < eY; ++y) {
+ for (int x = bX; x < eX; ++x) {
+
+ latticeIndex = L_INDEX_4(lDims, x, y, z);
+
+ if (ld->Lattice[latticeIndex] != LAT_CELL_OBSTACLE) {
+ grid[latticeIndex] = counter;
+
+ coords[C_INDEX_X(counter)] = x;
+ coords[C_INDEX_Y(counter)] = y;
+ coords[C_INDEX_Z(counter)] = z;
+
+ ++counter;
+ }
+ } } }
+ } } }
+
+ Verify(counter == nFluid);
+
+ uint32_t * adjList;
+
+ // AoSoA addressing for adjList needs padding for (nFluid % VSIZE) != 0
+ int nFluid_padded = ((nFluid + VSIZE - 1) / VSIZE) * VSIZE;
+ unsigned long adjListBytes = nFluid_padded * sizeof(int) * N_D3Q19_IDX;
+
+ printf("# Adjacency List Allocation:\n");
+ printf("# size: \t\t%e MiB\n", adjListBytes / 1024.0 / 1024.0);
+ printf("# alignment: \t\t%d b\n", PAGE_4K);
+
+ // AdjList only requires 18 instead of 19 entries per node, as
+ // the center PDF needs no addressing.
+ if (ADJ_LIST_ALLOCATOR((void **)&adjList, adjListBytes, PAGE_4K)) {
+ printf("ERROR: allocating adjList array with %s() failed: %lu bytes.\n", STRINGIFY(ADJ_LIST_ALLOCATOR), adjListBytes);
+ exit(1);
+ }
+ else {
+ printf("# allocator: \t\t\t%s()\n", STRINGIFY(ADJ_LIST_ALLOCATOR));
+ }
+
+ for (int i = 0; i < nFluid_padded; ++i){
+ adjList[i] = -1;
+ }
+
+ kdl->AdjList = adjList;
+
+ int x, y, z;
+
+ uint32_t neighborIndex;
+ uint32_t dstIndex;
+
+ int nx, ny, nz, px, py, pz;
+
+ // Loop over all fluid nodes and compute the indices to the neighboring
+ // PDFs for configured data layout (AoS/SoA).
+ // Parallelized loop to ensure correct NUMA placement.
+ // #ifdef _OPENMP --> add line continuation
+ // #pragma omp parallel for default(none)
+ // shared(nFluid, nCells, coords, D3Q19_INV, D3Q19_X, D3Q19_Y, D3Q19_Z,
+ // stderr,
+ // lDims, grid, ld, lX, lY, lZ, adjList)
+ // private(x, y, z, nx, ny, nz, neighborIndex, dstIndex)
+ // #endif
+
+ for (int fluidBaseIndex = 0; fluidBaseIndex < nFluid; fluidBaseIndex+=VSIZE) {
+
+
+ // Loop over all directions except the center one.
+ for(int d = 0; d < N_D3Q19 - 1; ++d) {
+ Assert(d != D3Q19_C);
+
+ for(int inChunkIndex = 0; (inChunkIndex < VSIZE) && ((fluidBaseIndex + inChunkIndex) < nFluid); ++inChunkIndex){
+ int index = fluidBaseIndex + inChunkIndex;
+
+ Assert(index < nFluid);
+
+ x = coords[C_INDEX_X(index)];
+ y = coords[C_INDEX_Y(index)];
+ z = coords[C_INDEX_Z(index)];
+
+ Assert(x >= 0 && x < lX);
+ Assert(y >= 0 && y < lY);
+ Assert(z >= 0 && z < lZ);
+
+ Assert(ld->Lattice[L_INDEX_4(lDims, x, y, z)] != LAT_CELL_OBSTACLE);
+
+#ifdef PROP_MODEL_PUSH
+ nx = x + D3Q19_X[d];
+ ny = y + D3Q19_Y[d];
+ nz = z + D3Q19_Z[d];
+
+#elif PROP_MODEL_PULL
+ nx = x - D3Q19_X[d];
+ ny = y - D3Q19_Y[d];
+ nz = z - D3Q19_Z[d];
+#else
+ #error No implementation for this PROP_MODEL_NAME.
+#endif
+ // If the neighbor is outside the latice in X direction and we have a
+ // periodic boundary then we need to wrap around.
+ if ( ((nx < 0 || nx >= lX) && ld->PeriodicX) ||
+ ((ny < 0 || ny >= lY) && ld->PeriodicY) ||
+ ((nz < 0 || nz >= lZ) && ld->PeriodicZ)
+ ){
+ // x periodic
+
+ if (nx < 0) {
+ px = lX - 1;
+ }
+ else if (nx >= lX) {
+ px = 0;
+ } else {
+ px = nx;
+ }
+ // y periodic
+ if (ny < 0) {
+ py = lY - 1;
+ }
+ else if (ny >= lY) {
+ py = 0;
+ } else {
+ py = ny;
+ }
+
+ // z periodic
+ if (nz < 0) {
+ pz = lZ - 1;
+ }
+ else if (nz >= lZ) {
+ pz = 0;
+ } else {
+ pz = nz;
+ }
+
+ if (ld->Lattice[L_INDEX_4(lDims, px, py, pz)] == LAT_CELL_OBSTACLE) {
+ dstIndex = P_INDEX_3(nCells, index, D3Q19_INV[d]);
+ }
+ else {
+ neighborIndex = grid[L_INDEX_4(lDims, px, py, pz)];
+
+ AssertMsg(neighborIndex != ~0, "Neighbor has no Index. (%d %d %d) direction %s (%d)\n", px, py, pz, D3Q19_NAMES[d], d);
+
+ dstIndex = P_INDEX_3(nCells, neighborIndex, d);
+ }
+ }
+ else if (nx < 0 || ny < 0 || nz < 0 || nx >= lX || ny >= lY || nz >= lZ) {
+ dstIndex = P_INDEX_3(nCells, index, D3Q19_INV[d]);
+ }
+ else if (ld->Lattice[L_INDEX_4(lDims, nx, ny, nz)] == LAT_CELL_OBSTACLE) {
+ dstIndex = P_INDEX_3(nCells, index, D3Q19_INV[d]);
+ }
+ else {
+ neighborIndex = grid[L_INDEX_4(lDims, nx, ny, nz)];
+
+ Assert(neighborIndex != ~0);
+
+ dstIndex = P_INDEX_3(nCells, neighborIndex, d);
+ }
+
+ Assert(dstIndex >= 0);
+ Assert(dstIndex < nCells * N_D3Q19);
+
+ adjList[(index - (index % VSIZE)) * N_D3Q19_IDX + (d * VSIZE) + (index % VSIZE)] = dstIndex;
+ }
+ }
+ }
+
+// Sets unused adjList entries to some extreme value which triggers and SIGSEG, whenever these values are accidently accessed.
+ for(int index = nFluid; index < nFluid_padded; ++index){
+ for(int d = 0; d < N_D3Q19 - 1; ++d) {
+ adjList[(index - (index % VSIZE)) * N_D3Q19_IDX + (d * VSIZE) + (index % VSIZE)] = -10*1000*1000;
+ }
+ }
+
+/*
+ printf("============\n");
+ for (int i = 0; i < nFluid_padded * (N_D3Q19_IDX + 20);){
+ for (int j = 0; j < VSIZE; ++j){
+ printf("%d ",adjList[i]);
+ ++i;
+ }
+ printf("\n");
+ }
+ for(int dir = 0; dir < N_D3Q19; ++dir){
+ printf("dir: %d\n",dir);
+ for(int baseIndex = 0; baseIndex < nFluid + VSIZE; baseIndex+=VSIZE){
+ for(int i = 0; i < VSIZE; ++i){
+ int index = baseIndex + i;
+
+ printf("%d ", adjList[(index - (index % VSIZE)) * N_D3Q19_IDX + (dir * VSIZE) + (index % VSIZE)]);
+ }
+ printf("\n");
+ }
+ printf("\n");
+ }
+ printf("============\n");
+*/
+
+
+ int nThreads = 1;
+
+#ifdef _OPENMP
+ nThreads = omp_get_max_threads();
+#endif
+
+ SetuploopStartIndices(ld, KDLR(kd), nThreads);
+
+ // Fill remaining KernelData structures
+ kd->GetNode = GetNode;
+ kd->SetNode = SetNode;
+
+ kd->BoundaryConditionsGetPdf = FNAME(BCGetPdf);
+ kd->BoundaryConditionsSetPdf = FNAME(BCSetPdf);
+
+ kd->Kernel = FNAME(D3Q19ListAaPvGatherHybridKernel);
+
+ kd->DstPdfs = NULL;
+ kd->PdfsActive = kd->Pdfs[0];
+
+ return;
+}
+
+void FNAME(D3Q19ListAaPvGatherHybridDeinit)(LatticeDesc * ld, KernelData ** kernelData)
+{
+ KernelDataListRia ** kdlr = (KernelDataListRia **)kernelData;
+
+ MemFree((void **)&((*kdlr)->loopStartIndices));
+
+ if ((*kdlr)->oddKernelThreadStartIndices != NULL) {
+ MemFree((void **)&((*kdlr)->oddKernelThreadStartIndices));
+ }
+
+ KernelDataList ** kdl = (KernelDataList **)kernelData;
+
+ ADJ_LIST_FREE((void **)&((*kdl)->AdjList));
+
+ MemFree((void **)&((*kdl)->Coords));
+ MemFree((void **)&((*kdl)->Grid));
+
+ PDF_FREE((void **)&((*kernelData)->Pdfs[0]));
+
+ MemFree((void **)kernelData);
+ return;
+}
+#undef PAGE_4K
+#undef ADJ_LIST
projects / LbmBenchmarkKernelsPublic.git / blobdiff