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 // Michael Hussnaetter, 2017-2018
12 // University of Erlangen-Nuremberg, Germany
13 // michael.hussnaetter -at- fau.de
15 // This file is part of the Lattice Boltzmann Benchmark Kernels (LbmBenchKernels).
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.
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.
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/>.
30 // --------------------------------------------------------------------------
31 #include "BenchKernelD3Q19ListAaPvGatherHybridCommon.h"
45 #if ALLOC_ADJ_LIST_IN_HBM == 1
46 #define ADJ_LIST_ALLOCATOR HbwAllocAligned
47 #define ADJ_LIST_FREE HbwFree
49 #define ADJ_LIST_ALLOCATOR MemAllocAligned
50 #define ADJ_LIST_FREE MemFree
53 #if ALLOC_PDF_IN_HBM == 1
54 #define PDF_ALLOCATOR HbwAllocAligned
55 #define PDF_FREE HbwFree
57 #define PDF_ALLOCATOR MemAllocAligned
58 #define PDF_FREE MemFree
61 // Forward definition.
62 void FNAME(D3Q19ListAaPvGatherHybridKernel)(LatticeDesc * ld, struct KernelData_ * kd, CaseData * cd);
66 // -----------------------------------------------------------------------
67 // Functions which are used as callback by the kernel to read or write
70 static void FNAME(BCGetPdf)(KernelData * kd, int x, int y, int z, int dir, PdfT * pdf)
73 Assert(kd->PdfsActive != NULL);
74 Assert(kd->PdfsActive == kd->Pdfs[0] || kd->PdfsActive == kd->Pdfs[1]);
77 Assert(x >= 0); Assert(y >= 0); Assert(z >= 0);
78 Assert(x < kd->Dims[0]); Assert(y < kd->Dims[1]); Assert(z < kd->Dims[2]);
79 Assert(dir >= 0); Assert(dir < N_D3Q19);
81 KernelDataList * kdl = KDL(kd);
82 uint32_t * adjList = kdl->AdjList;
84 if (kdl->Iteration % 2 == 0) {
85 // Pdfs are stored inverse, local PDFs are located in remote nodes
88 uint32_t index = kdl->Grid[L_INDEX_4(kd->Dims, x, y, z)];
91 #define ADJ_LIST(dir) adjList[(index - (index % VSIZE)) * N_D3Q19_IDX + (dir * VSIZE) + (index % VSIZE)]
92 *pdf = kd->PdfsActive[ADJ_LIST(D3Q19_INV[dir])];
96 *pdf = kd->PdfsActive[P_INDEX_3(kdl->nCells, index, dir)];
101 *pdf = kd->PdfsActive[P_INDEX_5(kdl, x, y, z, dir)];
107 static void FNAME(BCSetPdf)(KernelData * kd, int x, int y, int z, int dir, PdfT pdf)
110 Assert(kd->PdfsActive != NULL);
111 Assert(kd->PdfsActive == kd->Pdfs[0] || kd->PdfsActive == kd->Pdfs[1]);
112 Assert(x >= 0); Assert(y >= 0); Assert(z >= 0);
113 Assert(x < kd->Dims[0]); Assert(y < kd->Dims[1]); Assert(z < kd->Dims[2]);
114 Assert(dir >= 0); Assert(dir < N_D3Q19);
117 printf("ERROR: setting nan %d %d %d %d %s\n", x, y, z, dir, D3Q19_NAMES[dir]);
122 KernelDataList * kdl = KDL(kd);
123 uint32_t * adjList = kdl->AdjList;
125 if (kdl->Iteration % 2 == 0) {
126 // Pdfs are stored inverse, local PDFs are located in remote nodes
128 // getting node index
129 uint32_t index = kdl->Grid[L_INDEX_4(kd->Dims, x, y, z)];
131 if (dir != D3Q19_C) {
132 #define ADJ_LIST(dir) adjList[(index - (index % VSIZE)) * N_D3Q19_IDX + (dir * VSIZE) + (index % VSIZE)]
133 kd->PdfsActive[ADJ_LIST(D3Q19_INV[dir])] = pdf;
136 kd->PdfsActive[P_INDEX_3(kdl->nCells, index, dir)] = pdf;
140 kd->PdfsActive[P_INDEX_5(kdl, x, y, z, dir)] = pdf;
147 static void GetNode(KernelData * kd, int x, int y, int z, PdfT * pdfs)
150 Assert(kd->PdfsActive != NULL);
151 Assert(kd->PdfsActive == kd->Pdfs[0] || kd->PdfsActive == kd->Pdfs[1]);
152 Assert(pdfs != NULL);
153 Assert(x >= 0); Assert(y >= 0); Assert(z >= 0);
154 Assert(x < kd->Dims[0]); Assert(y < kd->Dims[1]); Assert(z < kd->Dims[2]);
156 KernelDataList * kdl = KDL(kd);
157 uint32_t * adjList = kdl->AdjList;
159 if(kdl->Iteration % 2 == 0){
161 uint32_t index = kdl->Grid[L_INDEX_4(kdl->kd.Dims, x, y, z)];
163 // Load PDFs of local cell: pdf_N = src[adjList[adjListIndex + D3Q19_S]]; ...
164 #define ADJ_LIST(dir) adjList[(index - (index % VSIZE)) * N_D3Q19_IDX + (dir * VSIZE) + (index % VSIZE)]
165 #define X(name, idx, idxinv, _x, _y, _z) pdfs[idx] = kd->PdfsActive[ADJ_LIST(idxinv)];
169 pdfs[D3Q19_C] = kd->PdfsActive[P_INDEX_3(kdl->nCells, index, D3Q19_C)];
173 #define I(x, y, z, dir) P_INDEX_5(KDL(kd), (x), (y), (z), (dir))
174 #define X(name, idx, idxinv, _x, _y, _z) pdfs[idx] = kd->PdfsActive[I(x, y, z, idx)];
181 for (int d = 0; d < 19; ++d) {
182 if(isnan(pdfs[d]) || isinf(pdfs[d])) {
183 printf("%d %d %d %d nan! get node\n", x, y, z, d);
184 for (int d2 = 0; d2 < 19; ++d2) {
185 printf("%d: %e\n", d2, pdfs[d2]);
195 static void SetNode(KernelData * kd, int x, int y, int z, PdfT * pdfs)
198 Assert(kd->PdfsActive != NULL);
199 Assert(kd->PdfsActive == kd->Pdfs[0] || kd->PdfsActive == kd->Pdfs[1]);
200 Assert(pdfs != NULL);
202 Assert(x >= 0); Assert(y >= 0); Assert(z >= 0);
203 Assert(x < kd->Dims[0]); Assert(y < kd->Dims[1]); Assert(z < kd->Dims[2]);
205 for (int d = 0; d < 19; ++d) {
207 printf("%d %d %d %d nan! get node\n", x, y, z, d);
208 for (int d2 = 0; d2 < 19; ++d2) {
209 printf("%d: %e\n", d2, pdfs[d2]);
215 KernelDataList * kdl = KDL(kd);
216 uint32_t * adjList = kdl->AdjList;
218 if(kdl->Iteration % 2 == 0){
220 uint32_t index = kdl->Grid[L_INDEX_4(kdl->kd.Dims, x, y, z)];
222 // Load PDFs of local cell: pdf_N = src[adjList[adjListIndex + D3Q19_S]]; ...
223 kd->PdfsActive[P_INDEX_3(kdl->nCells, index, D3Q19_C)] = pdfs[D3Q19_C];
225 #define ADJ_LIST(dir) adjList[(index - (index % VSIZE)) * N_D3Q19_IDX + (dir * VSIZE) + (index % VSIZE)]
226 #define X(name, idx, idxinv, _x, _y, _z) kd->PdfsActive[ADJ_LIST(idxinv)] = pdfs[idx];
233 #define I(x, y, z, dir) P_INDEX_5(KDL(kd), (x), (y), (z), (dir))
234 #define X(name, idx, idxinv, _x, _y, _z) kd->PdfsActive[I(x, y, z, idx)] = pdfs[idx];
244 static void ParameterUsage()
246 printf("Kernel parameters:\n");
247 printf(" [-blk <n>] [-blk-[xyz] <n>]\n");
252 static void ParseParameters(Parameters * params, int * blk)
256 blk[0] = 0; blk[1] = 0; blk[2] = 0;
258 #define ARG_IS(param) (!strcmp(params->KernelArgs[i], param))
259 #define NEXT_ARG_PRESENT() \
261 if (i + 1 >= params->nKernelArgs) { \
262 printf("ERROR: argument %s requires a parameter.\n", params->KernelArgs[i]); \
268 for (int i = 0; i < params->nKernelArgs; ++i) {
269 if (ARG_IS("-blk") || ARG_IS("--blk")) {
272 int tmp = strtol(params->KernelArgs[++i], NULL, 0);
275 printf("ERROR: blocking parameter must be > 0.\n");
279 blk[0] = blk[1] = blk[2] = tmp;
281 else if (ARG_IS("-blk-x") || ARG_IS("--blk-x")) {
284 int tmp = strtol(params->KernelArgs[++i], NULL, 0);
287 printf("ERROR: blocking parameter must be > 0.\n");
293 else if (ARG_IS("-blk-y") || ARG_IS("--blk-y")) {
296 int tmp = strtol(params->KernelArgs[++i], NULL, 0);
299 printf("ERROR: blocking parameter must be > 0.\n");
305 else if (ARG_IS("-blk-z") || ARG_IS("--blk-z")) {
308 int tmp = strtol(params->KernelArgs[++i], NULL, 0);
311 printf("ERROR: blocking parameter must be > 0.\n");
317 else if (ARG_IS("-h") || ARG_IS("-help") || ARG_IS("--help")) {
322 printf("ERROR: unknown kernel parameter.\n");
329 #undef NEXT_ARG_PRESENT
334 static void SetuploopStartIndices(LatticeDesc * ld, KernelDataListRia * kdlr, int nThreads)
336 //#define ADJ_LIST(dir) adjList[(index - (index % VSIZE)) * N_D3Q19_IDX + (dir * VSIZE) + (index % VSIZE)]
338 Assert(kdlr != NULL);
339 Assert(nThreads > 0);
341 //uint32_t * adjList = kdlr->kdl.AdjList;
342 uint32_t * adjList = kdlr->kdl.AdjList;
344 uint32_t nLoopStartIndices = 0;
345 uint32_t loopStartIndex = 2;
347 int nFluid = kdlr->kdl.nFluid;
349 int * oddKernelThreadStartIndices = (int *)malloc(sizeof(int) * (nThreads + 1));
351 int nNodesPerThread = nFluid / nThreads;
352 //printf("nodesPerThread: %d\n", nNodesPerThread);
354 for (int i = 0; i < nThreads; ++i) {
355 oddKernelThreadStartIndices[i] = i * nNodesPerThread + MinI(i, nFluid % nThreads);
358 oddKernelThreadStartIndices[nThreads] = nFluid;
361 for (int i = 0; i <= nThreads; ++i) {
362 printf("oddKernelThreadStartIndices[%d] = %d\n", i, oddKernelThreadStartIndices[i]);
368 // We execute following code two times.
369 // - The first time to get the count of how many entries we need for the
370 // loopStartIndices array.
371 // - The second time to fill the array.
373 // Loop over adjacency list of all nodes.
374 // Compare if adjacent nodes within one cache line share the same access pattern.
375 // First vectorized access is assumed to be consecutive (-> may be loaded with regular load).
377 int lastCacheLineConsecutive = 1;
379 for (int fluidBaseIndex = 1; fluidBaseIndex < nFluid + 1; fluidBaseIndex += VSIZE) {
381 int currentCacheLineConsecutive = 1;
383 // Loop over all directions except the center one.
384 for(int d = 0; d < N_D3Q19 - 1; ++d) {
385 Assert(d != D3Q19_C);
387 // check if cache line itself has consecutive memory access pattern
388 for(int inChunkIndex = 0; (inChunkIndex < VSIZE - 1) && ((fluidBaseIndex + inChunkIndex) < nFluid); ++inChunkIndex) {
389 int index = fluidBaseIndex + inChunkIndex;
391 Assert(index < nFluid);
393 #define ADJ_LIST(idx, dir) adjList[((idx) - ((idx) % VSIZE)) * N_D3Q19_IDX + ((dir) * VSIZE) + ((idx) % VSIZE)]
394 //if (adjList[index * N_D3Q19_IDX + d] != adjList[(index - 1) * N_D3Q19_IDX + d] + 1)
395 if (ADJ_LIST(index, d) != ADJ_LIST(index-1, d) + 1) {
396 //printf("no match for index: %d\n", d);
397 //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);
398 // Different access pattern.
399 currentCacheLineConsecutive = 0;
406 if(!currentCacheLineConsecutive){
407 break; //exit from nested loop
412 int interCacheLineConsecutive = 1;
414 if(currentCacheLineConsecutive && lastCacheLineConsecutive){
415 // check if cache line has consecutive memory access pattern to last entry of previous cache line
416 int lastIdxOfPreviousCacheLine = fluidBaseIndex - 2;
417 if (lastIdxOfPreviousCacheLine > 0) {
418 for(int d = 0; d < N_D3Q19 - 1; ++d) {
419 Assert(d != D3Q19_C);
420 #define ADJ_LIST(idx, dir) adjList[((idx) - ((idx) % VSIZE)) * N_D3Q19_IDX + ((dir) * VSIZE) + ((idx) % VSIZE)]
421 if (ADJ_LIST(fluidBaseIndex-1, d) != ADJ_LIST(lastIdxOfPreviousCacheLine, d) + 1) {
422 // Different access pattern.
423 //printf("not interCacheConsecutive\n");
424 interCacheLineConsecutive = 0;
432 int threadBoundaryIndex = oddKernelThreadStartIndices[threadIndex];
433 if (fluidBaseIndex - 1 <= threadBoundaryIndex &&
434 threadBoundaryIndex < fluidBaseIndex + VSIZE - 1) {
435 // Current cache line contains thread boundary.
436 // These cache lines are treated by scalar peel and
437 // reminder loops in kernel of every thread.
438 // TODO maybe replace these loops with masked gather / scatter?!
439 if (loopStartIndex % 2 == 0) { // current index would be gather/scatter index
440 ++loopStartIndex; // reserving gather/scatter index
442 ++loopStartIndex; // reserving space for start remainder loop of thread n
444 if (threadIndex < nThreads){
445 ++loopStartIndex; // reserving space for starting peel loop of thread n+1
447 if (fluidBaseIndex - 1 == threadBoundaryIndex){
448 if(!currentCacheLineConsecutive){
453 currentCacheLineConsecutive = 1;
456 //++loopStartIndex; // reserving space for ending peel loop / starting load/store of thread n+1
457 ++loopStartIndex; // 1st load/store end == 1st gather/scatter start OR 1st gather/scatter end == 2nd load/start start
462 // We are not at a thread boundary.
463 if (currentCacheLineConsecutive) {
464 if(lastCacheLineConsecutive && !interCacheLineConsecutive){
467 else if(!lastCacheLineConsecutive){
472 if(lastCacheLineConsecutive){
478 // treating special case when last thread has no remainder loop
479 if (oddKernelThreadStartIndices[nThreads] == fluidBaseIndex + VSIZE - 1) {
480 //printf("--> special case 111. loopStartIndex: %d \n", loopStartIndex);
481 if (loopStartIndex % 2 != 0) { // current index is gather/scatter end and load/store start index
482 ++loopStartIndex; //set load/store end (gather/scatter start) to same value as scalar remainder start => no more access to gather/scatter loop
485 ++loopStartIndex; // gather/scatter end and scalar remainder start
486 ++loopStartIndex; // scalar remainder end and scalar peel start
490 lastCacheLineConsecutive = currentCacheLineConsecutive;
494 nLoopStartIndices = loopStartIndex;
497 int * loopStartIndices;
498 unsigned long loopStartIndicesByte = (nLoopStartIndices + 1) * sizeof(int);
500 printf("# Loop Start Index Array Allocation:\n");
501 printf("# elements: \t\t%d\n", nLoopStartIndices + 1);
502 printf("# size: \t\t%e MiB\n", loopStartIndicesByte / 1024.0 / 1024.0);
503 printf("# alignment: \t\t%d b\n", PAGE_4K);
505 if (MemAllocAligned((void **)&loopStartIndices, loopStartIndicesByte, PAGE_4K)) {
506 printf("ERROR: allocating loopStartIndices array with MemAllocAligned failed: %lu bytes.\n", loopStartIndicesByte);
510 printf("# allocator: \t\t\tMemAllocAligned()\n");
513 oddKernelThreadStartIndices[0] = 0;
514 loopStartIndices[0] = 0; //first scalar loop would start with 0
515 loopStartIndices[1] = 0; //no peel loop expected -> first load/store loop may start at index==0
516 loopStartIndices[2] = 0; //may not be set in case first access is gather/scatter -> therefore its set here
518 // resetting values to default
520 lastCacheLineConsecutive = 1;
523 // Loop over adjacency list of all nodes.
524 // Compare if adjacent nodes share the same access pattern.
526 int indexAccumulator = 0;
528 // for statistical reasons:
529 int gatherAccumulator = 0;
530 int loadAccumulator = 0;
531 int scalarLookups = 0;
535 for (int fluidBaseIndex = 1; fluidBaseIndex < nFluid + 1; fluidBaseIndex += VSIZE) {
536 int currentCacheLineConsecutive = 1;
537 //printf("fluidbaseIndex: %d\n", fluidBaseIndex);
538 // Loop over all directions except the center one.
539 for(int d = 0; d < N_D3Q19 - 1; ++d) {
540 Assert(d != D3Q19_C);
542 // check if cache line itself has consecutive memory access pattern
543 for(int inChunkIndex = 0; (inChunkIndex < VSIZE - 1) && ((fluidBaseIndex + inChunkIndex) < nFluid); ++inChunkIndex){
544 int index = fluidBaseIndex + inChunkIndex;
546 Assert(index < nFluid);
548 #define ADJ_LIST(idx, dir) adjList[((idx) - ((idx) % VSIZE)) * N_D3Q19_IDX + ((dir) * VSIZE) + ((idx) % VSIZE)]
549 //if (adjList[index * N_D3Q19_IDX + d] != adjList[(index - 1) * N_D3Q19_IDX + d] + 1)
550 if (ADJ_LIST(index, d) != ADJ_LIST(index-1, d) + 1) {
551 // Different access pattern.
552 currentCacheLineConsecutive = 0;
558 if(!currentCacheLineConsecutive){
559 break; //exit from nested loop
563 int interCacheLineConsecutive = 1;
565 if(currentCacheLineConsecutive && lastCacheLineConsecutive){
566 // check if cache line has consecutive memory access pattern to last entry of previous cache line
567 int lastIdxOfPreviousCacheLine = fluidBaseIndex - 2;
568 if (lastIdxOfPreviousCacheLine > 0) {
569 for(int d = 0; d < N_D3Q19 - 1; ++d) {
570 Assert(d != D3Q19_C);
571 #define ADJ_LIST(idx, dir) adjList[((idx) - ((idx) % VSIZE)) * N_D3Q19_IDX + ((dir) * VSIZE) + ((idx) % VSIZE)]
572 if (ADJ_LIST(fluidBaseIndex-1, d) != ADJ_LIST(lastIdxOfPreviousCacheLine, d) + 1) {
573 // Different access pattern.
574 interCacheLineConsecutive = 0;
583 int threadBoundaryIndex = oddKernelThreadStartIndices[threadIndex];
584 //if (fluidBaseIndex > 3500)
585 // printf("threadBoundaryIndex: %d fluidBaseIndex-1: %d fluidBaseIndex + VSIZE - 1: %d\n", threadBoundaryIndex, fluidBaseIndex-1, fluidBaseIndex + VSIZE -1);
587 if (fluidBaseIndex - 1 <= threadBoundaryIndex &&
588 threadBoundaryIndex < fluidBaseIndex + VSIZE - 1) {
589 // Current cache line contains thread boundary.
590 // These cache lines are treated by scalar peel and
591 // reminder loops in kernel of every thread.
592 // TODO maybe replace these loops with masked gather / scatter?!
593 if (loopStartIndex % 2 == 0) { // current index would be gather/scatter index
594 //loopStartIndices[loopStartIndex] = fluidBaseIndex - 1; //same value as scalar remainder start => no more access to gather/scatter loop
595 loopStartIndices[loopStartIndex] = indexAccumulator; //same value as scalar remainder start => no more access to gather/scatter loop
599 //loopStartIndices[loopStartIndex] = fluidBaseIndex - 1; // gather/scatter end and scalar remainder start
600 loopStartIndices[loopStartIndex] = indexAccumulator; // gather/scatter end and scalar remainder start
603 // starting indices of thread n+1
604 loopStartIndices[loopStartIndex] = threadBoundaryIndex; // scalar remainder of thread n end and scalar peel of thread n+1 start
605 oddKernelThreadStartIndices[threadIndex] = loopStartIndex; // thread start is where scalar peel starts
607 if (threadIndex < nThreads){
608 indexAccumulator = ((threadBoundaryIndex + VSIZE - 1) / VSIZE ) * VSIZE; // rounding towards next multiple of VSIZE
610 loopStartIndices[loopStartIndex] = indexAccumulator; // scalar peel end and 1st load/store start
612 // treating special case when there is no peel / remainder loop
613 if (fluidBaseIndex - 1 == threadBoundaryIndex){
614 if(!currentCacheLineConsecutive){
616 loopStartIndices[loopStartIndex] = indexAccumulator; // 1st load/store end and 1st gather/scatter start
617 gatherAccumulator += VSIZE;
620 loadLookups += VSIZE;
622 indexAccumulator += VSIZE;
625 scalarLookups += VSIZE;
626 currentCacheLineConsecutive = 1;
629 ++loopStartIndex; // 1st load/store end == 1st gather/scatter start OR 1st gather/scatter end == 2nd load/start start
630 loopStartIndices[loopStartIndex] = indexAccumulator; // 1st load/store end == 1st gather/scatter start OR 1st gather/scatter end == 2nd load/start start
636 // We are not at a thread boundary.
638 if (currentCacheLineConsecutive) {
639 loadAccumulator += VSIZE;
641 if(lastCacheLineConsecutive && !interCacheLineConsecutive){
642 loadLookups += VSIZE;
644 printf("#1 loopStartIndex: %d\n", loopStartIndex);
645 // loopStartIndices[loopStartIndex] is not incremented since pointers need to be fetched again.
646 // loopStartIndices[loopStartIndex + 1] (-> start Load/Store and end Gather/Scatter)
647 // gets same value as loopStartIndices[loopStartindex] (-> start Gather/Scatter)
648 // this ensures that no gather/scatter iteration is executed
650 loopStartIndices[loopStartIndex] = indexAccumulator;
652 // loopStartIndices[loopStartIndex + 2] (-> start Gather/Scatter and end Load/Store)
653 // gets set to have one Load/Store iteration
655 indexAccumulator+=VSIZE;
656 loopStartIndices[loopStartIndex] = indexAccumulator;
659 else if(!lastCacheLineConsecutive){
660 loadLookups += VSIZE;
662 printf("#2 loopStartIndex: %d\n", loopStartIndex);
664 indexAccumulator+=VSIZE;
665 loopStartIndices[loopStartIndex] = indexAccumulator;
667 else { // (lastCacheLineConsecutive && interCacheLineConsecutive)
669 printf("#3 loopStartIndex: %d\n", loopStartIndex);
670 indexAccumulator+=VSIZE;
671 loopStartIndices[loopStartIndex] = indexAccumulator;
675 gatherAccumulator += VSIZE;
676 if(lastCacheLineConsecutive){
678 printf("#4 loopStartIndex: %d\n", loopStartIndex);
680 indexAccumulator+=VSIZE;
681 loopStartIndices[loopStartIndex] = indexAccumulator;
683 else { // lastCacheLine without not consecutive memory access pattern
685 printf("#5 loopStartIndex: %d\n", loopStartIndex);
686 indexAccumulator+=VSIZE;
687 loopStartIndices[loopStartIndex] = indexAccumulator;
692 // treating special case when last thread has no remainder loop
693 if (oddKernelThreadStartIndices[nThreads] == fluidBaseIndex + VSIZE - 1) {
694 //printf("--> special case. indexAccumulator: %d\n", indexAccumulator);
695 if (loopStartIndex % 2 != 0) { // current index is gather/scatter end and load/store start index
697 loopStartIndices[loopStartIndex] = indexAccumulator; //set load/store end (gather/scatter start) to same value as scalar remainder start => no more access to gather/scatter loop
701 loopStartIndices[loopStartIndex] = indexAccumulator; // gather/scatter end and scalar remainder start
703 loopStartIndices[loopStartIndex] = indexAccumulator; // scalar remainder end and scalar peel start
705 oddKernelThreadStartIndices[threadIndex] = loopStartIndex; // thread start is where scalar peel starts
708 lastCacheLineConsecutive = currentCacheLineConsecutive;
712 if (nLoopStartIndices != loopStartIndex){
713 printf("ERROR: nLoopStartIndices unequal loopStartIndex!\n");
717 printf("loopStartIndices:\n");
718 for(int i = 0; i <= nLoopStartIndices; ++i){
719 printf("%d ", loopStartIndices[i]);
722 printf("oddKernelThreadStartIndices:\n");
723 for(int i = 0; i <= nThreads; ++i){
724 printf("%d ", oddKernelThreadStartIndices[i]);
729 kdlr->loopStartIndices = loopStartIndices;
730 kdlr->nLoopStartIndices = nLoopStartIndices;
732 kdlr->oddKernelThreadStartIndices = oddKernelThreadStartIndices;
733 kdlr->nOddKernelThreadStartIndices = nThreads;
735 printf("# vload/vstore nodes: \t% 10d \t(%3.4f %% of total fluid nodes)\n", loadAccumulator, ((double) loadAccumulator / (double) nFluid) * 100);
736 printf("# gather/scatter nodes:\t% 10d \t(%3.4f %% of total fluid nodes)\n", gatherAccumulator, ((double) gatherAccumulator / (double) nFluid) * 100.0);
737 printf("# vload/vstore lookups:\t% 10d \n", loadLookups * (N_D3Q19 - 1));
738 printf("# gather/scatter lookups:\t% 10d \n", gatherAccumulator * (N_D3Q19 - 1));
739 printf("# scalar lookups: \t% 10d \n", scalarLookups * (N_D3Q19 - 1));
741 double loopBalanceEven = 2.0 * 19 * sizeof(PdfT);
742 double loopBalanceOdd = 2.0 * 19 * sizeof(PdfT) /* actual PDFs */
743 + (((double)(gatherAccumulator + loadLookups + scalarLookups)) / nFluid) * sizeof(int) * (N_D3Q19 - 1) /* AdjList */
744 + (nLoopStartIndices / nFluid) * sizeof(int); // one lookup to loopStartIndices
746 double loopBalance = (loopBalanceEven + loopBalanceOdd) / 2.0;
748 kdlr->kdl.kd.LoopBalance = loopBalance;
750 printf("# loop balance:\n");
751 printf("# even timestep: \t\t%.2f B/FLUP\n", loopBalanceEven);
752 printf("# odd timestep: \t\t%.2f B/FLUP\n", loopBalanceOdd);
753 printf("# average: \t\t%.2f B/FLUP\n", loopBalance);
758 void FNAME(D3Q19ListAaPvGatherHybridInit)(LatticeDesc * ld, KernelData ** kernelData, Parameters * params)
761 KernelDataList * kdl;
762 KernelDataListRia * kdlr;
763 MemAlloc((void **)&kdlr, sizeof(KernelDataListRia));
765 kd = (KernelData *)kdlr;
773 kd->PdfsActive = NULL;
779 kd->GlobalDims[0] = -1;
780 kd->GlobalDims[1] = -1;
781 kd->GlobalDims[2] = -1;
786 kd->ObstIndices = NULL;
787 kd->nObstIndices = -1;
788 kd->BounceBackPdfsSrc = NULL;
789 kd->BounceBackPdfsDst = NULL;
790 kd->nBounceBackPdfs = -1;
798 kdlr->loopStartIndices = NULL;
799 kdlr->nLoopStartIndices = 0;
800 kdlr->oddKernelThreadStartIndices = NULL;
801 kdlr->nOddKernelThreadStartIndices = 0;
806 // Ajust the dimensions according to padding, if used.
807 kd->Dims[0] = kd->GlobalDims[0] = ld->Dims[0];
808 kd->Dims[1] = kd->GlobalDims[1] = ld->Dims[1];
809 kd->Dims[2] = kd->GlobalDims[2] = ld->Dims[2];
811 int * lDims = ld->Dims;
817 int nTotalCells = lX * lY * lZ;
818 int nCells = ld->nFluid; // TODO: + padding
819 int nFluid = ld->nFluid;
821 // TODO: check nCells/nFluid do not exceed 2^31. This actually has to be
822 // done during lattice setup.
823 kdl->nCells = nCells;
824 kdl->nFluid = nFluid;
830 ParseParameters(params, blk);
832 if (blk[0] == 0) blk[0] = lX;
833 if (blk[1] == 0) blk[1] = lY;
834 if (blk[2] == 0) blk[2] = lZ;
836 printf("# blocking: \t\tx: %3d y: %3d z: %3d\n", blk[0], blk[1], blk[2]);
838 unsigned long latByte = nCells * sizeof(PdfT) * N_D3Q19;
839 unsigned long latFluidByte = nFluid * sizeof(PdfT) * N_D3Q19;
840 unsigned long latPadByte = (nCells - nFluid) * sizeof(PdfT) * N_D3Q19;
842 printf("# Lattice Array Allocation:\n");
843 printf("# lattice size: \t\t%e MiB\n", latByte / 1024.0 / 1024.0);
844 printf("# fluid lattice size:\t\t%e MiB\n", latFluidByte / 1024.0 / 1024.0);
845 printf("# lattice padding: \t\t%e MiB\n", latPadByte / 1024.0 / 1024.0);
848 printf("# alignment: \t\t%d b\n", PAGE_4K);
850 if (PDF_ALLOCATOR((void **)&pdfs[0], latFluidByte, PAGE_4K)) {
851 printf("ERROR: allocating PDF array with %s() failed: %lu bytes.\n", STRINGIFY(PDF_ALLOCATOR), latFluidByte);
855 printf("# allocator: \t\t\t%s()\n", STRINGIFY(PDF_ALLOCATOR));
858 kd->Pdfs[0] = pdfs[0];
860 // Initialize PDFs with some (arbitrary) data for correct NUMA placement.
861 // Here we touch only the fluid nodes as this loop is OpenMP parallel and
862 // we want the same scheduling as in the kernel.
864 #pragma omp parallel for
866 for (int i = 0; i < nFluid; ++i) { for(int d = 0; d < N_D3Q19; ++d) {
867 pdfs[0][P_INDEX_3(nCells, i, d)] = 1.0;
870 // Initialize all PDFs to some standard value.
871 for (int i = 0; i < nFluid; ++i) { for(int d = 0; d < N_D3Q19; ++d) {
872 pdfs[0][P_INDEX_3(nCells, i, d)] = 0.0;
875 // ----------------------------------------------------------------------
876 // create grid which will hold the index numbers of the fluid nodes
880 if (MemAlloc((void **)&grid, nTotalCells * sizeof(uint32_t))) {
881 printf("ERROR: allocating grid for numbering failed: %lu bytes.\n", nTotalCells * sizeof(uint32_t));
889 for(int z = 0; z < lZ; ++z) {
890 for(int y = 0; y < lY; ++y) {
891 for(int x = 0; x < lX; ++x) {
893 latticeIndex = L_INDEX_4(ld->Dims, x, y, z);
895 grid[latticeIndex] = ~0;
901 // ----------------------------------------------------------------------
902 // generate numbering over grid
906 if (MemAlloc((void **)&coords, nFluid * sizeof(uint32_t) * 3)) {
907 printf("ERROR: allocating coords array failed: %lu bytes.\n", nFluid * sizeof(uint32_t) * 3);
911 kdl->Coords = coords;
913 // Index for the PDF nodes can start at 0 as we distinguish solid and fluid nodes
914 // through the ld->Lattice array.
917 // Blocking is implemented via setup of the adjacency list. The kernel later will
918 // walk through the lattice blocked automatically.
919 for (int bZ = 0; bZ < lZ; bZ += blk[2]) {
920 for (int bY = 0; bY < lY; bY += blk[1]) {
921 for (int bX = 0; bX < lX; bX += blk[0]) {
923 int eX = MIN(bX + blk[0], lX);
924 int eY = MIN(bY + blk[1], lY);
925 int eZ = MIN(bZ + blk[2], lZ);
928 for (int z = bZ; z < eZ; ++z) {
929 for (int y = bY; y < eY; ++y) {
930 for (int x = bX; x < eX; ++x) {
932 latticeIndex = L_INDEX_4(lDims, x, y, z);
934 if (ld->Lattice[latticeIndex] != LAT_CELL_OBSTACLE) {
935 grid[latticeIndex] = counter;
937 coords[C_INDEX_X(counter)] = x;
938 coords[C_INDEX_Y(counter)] = y;
939 coords[C_INDEX_Z(counter)] = z;
946 Verify(counter == nFluid);
950 // AoSoA addressing for adjList needs padding for (nFluid % VSIZE) != 0
951 int nFluid_padded = ((nFluid + VSIZE - 1) / VSIZE) * VSIZE;
952 unsigned long adjListBytes = nFluid_padded * sizeof(int) * N_D3Q19_IDX;
954 printf("# Adjacency List Allocation:\n");
955 printf("# size: \t\t%e MiB\n", adjListBytes / 1024.0 / 1024.0);
956 printf("# alignment: \t\t%d b\n", PAGE_4K);
958 // AdjList only requires 18 instead of 19 entries per node, as
959 // the center PDF needs no addressing.
960 if (ADJ_LIST_ALLOCATOR((void **)&adjList, adjListBytes, PAGE_4K)) {
961 printf("ERROR: allocating adjList array with %s() failed: %lu bytes.\n", STRINGIFY(ADJ_LIST_ALLOCATOR), adjListBytes);
965 printf("# allocator: \t\t\t%s()\n", STRINGIFY(ADJ_LIST_ALLOCATOR));
968 for (int i = 0; i < nFluid_padded; ++i){
972 kdl->AdjList = adjList;
976 uint32_t neighborIndex;
979 int nx, ny, nz, px, py, pz;
981 // Loop over all fluid nodes and compute the indices to the neighboring
982 // PDFs for configured data layout (AoS/SoA).
983 // Parallelized loop to ensure correct NUMA placement.
984 // #ifdef _OPENMP --> add line continuation
985 // #pragma omp parallel for default(none)
986 // shared(nFluid, nCells, coords, D3Q19_INV, D3Q19_X, D3Q19_Y, D3Q19_Z,
988 // lDims, grid, ld, lX, lY, lZ, adjList)
989 // private(x, y, z, nx, ny, nz, neighborIndex, dstIndex)
992 for (int fluidBaseIndex = 0; fluidBaseIndex < nFluid; fluidBaseIndex+=VSIZE) {
995 // Loop over all directions except the center one.
996 for(int d = 0; d < N_D3Q19 - 1; ++d) {
997 Assert(d != D3Q19_C);
999 for(int inChunkIndex = 0; (inChunkIndex < VSIZE) && ((fluidBaseIndex + inChunkIndex) < nFluid); ++inChunkIndex){
1000 int index = fluidBaseIndex + inChunkIndex;
1002 Assert(index < nFluid);
1004 x = coords[C_INDEX_X(index)];
1005 y = coords[C_INDEX_Y(index)];
1006 z = coords[C_INDEX_Z(index)];
1008 Assert(x >= 0 && x < lX);
1009 Assert(y >= 0 && y < lY);
1010 Assert(z >= 0 && z < lZ);
1012 Assert(ld->Lattice[L_INDEX_4(lDims, x, y, z)] != LAT_CELL_OBSTACLE);
1014 #ifdef PROP_MODEL_PUSH
1015 nx = x + D3Q19_X[d];
1016 ny = y + D3Q19_Y[d];
1017 nz = z + D3Q19_Z[d];
1019 #elif PROP_MODEL_PULL
1020 nx = x - D3Q19_X[d];
1021 ny = y - D3Q19_Y[d];
1022 nz = z - D3Q19_Z[d];
1024 #error No implementation for this PROP_MODEL_NAME.
1026 // If the neighbor is outside the latice in X direction and we have a
1027 // periodic boundary then we need to wrap around.
1028 if ( ((nx < 0 || nx >= lX) && ld->PeriodicX) ||
1029 ((ny < 0 || ny >= lY) && ld->PeriodicY) ||
1030 ((nz < 0 || nz >= lZ) && ld->PeriodicZ)
1037 else if (nx >= lX) {
1046 else if (ny >= lY) {
1056 else if (nz >= lZ) {
1062 if (ld->Lattice[L_INDEX_4(lDims, px, py, pz)] == LAT_CELL_OBSTACLE) {
1063 dstIndex = P_INDEX_3(nCells, index, D3Q19_INV[d]);
1066 neighborIndex = grid[L_INDEX_4(lDims, px, py, pz)];
1068 AssertMsg(neighborIndex != ~0, "Neighbor has no Index. (%d %d %d) direction %s (%d)\n", px, py, pz, D3Q19_NAMES[d], d);
1070 dstIndex = P_INDEX_3(nCells, neighborIndex, d);
1073 else if (nx < 0 || ny < 0 || nz < 0 || nx >= lX || ny >= lY || nz >= lZ) {
1074 dstIndex = P_INDEX_3(nCells, index, D3Q19_INV[d]);
1076 else if (ld->Lattice[L_INDEX_4(lDims, nx, ny, nz)] == LAT_CELL_OBSTACLE) {
1077 dstIndex = P_INDEX_3(nCells, index, D3Q19_INV[d]);
1080 neighborIndex = grid[L_INDEX_4(lDims, nx, ny, nz)];
1082 Assert(neighborIndex != ~0);
1084 dstIndex = P_INDEX_3(nCells, neighborIndex, d);
1087 Assert(dstIndex >= 0);
1088 Assert(dstIndex < nCells * N_D3Q19);
1090 adjList[(index - (index % VSIZE)) * N_D3Q19_IDX + (d * VSIZE) + (index % VSIZE)] = dstIndex;
1095 // Sets unused adjList entries to some extreme value which triggers and SIGSEG, whenever these values are accidently accessed.
1096 for(int index = nFluid; index < nFluid_padded; ++index){
1097 for(int d = 0; d < N_D3Q19 - 1; ++d) {
1098 adjList[(index - (index % VSIZE)) * N_D3Q19_IDX + (d * VSIZE) + (index % VSIZE)] = -10*1000*1000;
1103 printf("============\n");
1104 for (int i = 0; i < nFluid_padded * (N_D3Q19_IDX + 20);){
1105 for (int j = 0; j < VSIZE; ++j){
1106 printf("%d ",adjList[i]);
1111 for(int dir = 0; dir < N_D3Q19; ++dir){
1112 printf("dir: %d\n",dir);
1113 for(int baseIndex = 0; baseIndex < nFluid + VSIZE; baseIndex+=VSIZE){
1114 for(int i = 0; i < VSIZE; ++i){
1115 int index = baseIndex + i;
1117 printf("%d ", adjList[(index - (index % VSIZE)) * N_D3Q19_IDX + (dir * VSIZE) + (index % VSIZE)]);
1123 printf("============\n");
1130 nThreads = omp_get_max_threads();
1133 SetuploopStartIndices(ld, KDLR(kd), nThreads);
1135 // Fill remaining KernelData structures
1136 kd->GetNode = GetNode;
1137 kd->SetNode = SetNode;
1139 kd->BoundaryConditionsGetPdf = FNAME(BCGetPdf);
1140 kd->BoundaryConditionsSetPdf = FNAME(BCSetPdf);
1142 kd->Kernel = FNAME(D3Q19ListAaPvGatherHybridKernel);
1145 kd->PdfsActive = kd->Pdfs[0];
1150 void FNAME(D3Q19ListAaPvGatherHybridDeinit)(LatticeDesc * ld, KernelData ** kernelData)
1152 KernelDataListRia ** kdlr = (KernelDataListRia **)kernelData;
1154 MemFree((void **)&((*kdlr)->loopStartIndices));
1156 if ((*kdlr)->oddKernelThreadStartIndices != NULL) {
1157 MemFree((void **)&((*kdlr)->oddKernelThreadStartIndices));
1160 KernelDataList ** kdl = (KernelDataList **)kernelData;
1162 ADJ_LIST_FREE((void **)&((*kdl)->AdjList));
1164 MemFree((void **)&((*kdl)->Coords));
1165 MemFree((void **)&((*kdl)->Grid));
1167 PDF_FREE((void **)&((*kernelData)->Pdfs[0]));
1169 MemFree((void **)kernelData);