【发布时间】：2017-06-12 16:37:56
【问题描述】：
我有一个表示长方体中点的数组。它是一个一维数组，它使用以下索引函数来实现3维：
int getCellIndex(int ix, int iy, int iz) {
return ix + (iy * numCellsX) + (iz * numCellsX * numCellsY);
}
域中的单元格数为：
numCells = (numX + 2) * (numY + 2) * (numZ + 2)
其中 numX/numY/numZ 是 X/Y/Z 方向的单元格数。每个方向的 +2 是在域外部创建填充单元。每个方向的单元格数由下式给出：
numX = 5 * numY
numZ = numY/2
numY = userInput
对于每个单元格，我想根据它的邻居值（即模板）计算该单元格的新值，其中邻居在上方、下方、左侧、右侧、正面和背面。但是，我只想对还不错的单元格进行此计算。我有一个布尔数组来跟踪单元格是否坏。这是当前计算的样子：
for(int z = 1; z < numZ+1; z++) {
for(int y = 1; y < numY+1; y++) {
for(int x = 1; x < numX+1; x++) {
if(!isBadCell[ getCellIndex(x,y,z) ] {
// Do stencil Computation
}
}
}
}
这不是很好的表现。我希望能够对循环进行矢量化以提高性能，但是由于 if 语句我不能。我提前知道细胞是否坏了，这在整个计算过程中都不会改变。我想将域分成块，最好是 4x4x4 块，这样我可以先验地计算每个块是否包含坏细胞，如果是这样，照常处理它，或者如果不是，使用一个优化的函数，可以采取矢量化的优势例如
for(block : blocks) {
if(isBadBlock[block]) {
slowProcessBlock(block) // As above
} else {
fastVectorizedProcessBlock(block)
}
}
注意：块不需要物理存在，即这可以通过更改索引函数并使用不同的索引循环数组来实现。我愿意接受任何最有效的方法。
fastVectorizedProcessBlock() 函数看起来类似于 slowProcessBlock() 函数，但使用 if 语句删除（因为我们知道它不包含坏单元格）和矢量化 pragma。
如何将我的域拆分为多个块以便完成此任务？这似乎很棘手，因为 a) 每个方向上的单元格数量不相等，b) 我们需要考虑填充单元格，因为我们绝不能尝试计算它们的值，因为这会导致内存访问失效的界限。
如何在不使用 if 语句的情况下处理不包含坏单元格的块？
编辑：
这是我最初的想法：
for(int i = 0; i < numBlocks; i++) { // use blocks of 4x4x4 = 64
if(!isBadBlock[i]) {
// vectorization pragma here
for(int z = 0; z < 4; z++) {
for(int y = 0; y < 4; y++) {
for(int x = 0; x < 4; x++) {
// calculate stencil using getCellIndex(x,y,z)*i
}
}
}
} else {
for(int z = 0; z < 4; z++) {
for(int y = 0; y < 4; y++) {
for(int x = 0; x < 4; x++) {
if(!isBadCell[i*getCellIndex(x,y,z)]) {
// calculate stencil using getCellIndex(x,y,z)*i
}
}
}
}
}
现在单元格将存储在块中，即第一个 4x4x4 块中的所有单元格将存储在 pos 0-63 中，然后第二个块中的所有单元格将存储在 pos 64-127 等中。
但是，如果 numX/numY/numZ 值不正确，我认为不会起作用。例如，如果 numY = 2、numZ = 1 和 numX = 10 会怎样？ for 循环期望 z 方向至少有 4 个单元格深。有什么好办法可以解决吗？
更新 2 - 以下是模板计算的样子：
if ( isBadCell[ getCellIndex(x,y,z) ] ) {
double temp = someOtherArray[ getCellIndex(x,y,z) ] +
1.0/CONSTANT/CONSTANT*
(
- 1.0 * cells[ getCellIndex(x-1,y,z) ]
- 1.0 * cells[ getCellIndex(x+1,y,z) ]
- 1.0 * cells[ getCellIndex(x,y-1,z) ]
- 1.0 * cells[ getCellIndex(x,y+1,z) ]
- 1.0 * cells[ getCellIndex(x,y,z-1) ]
- 1.0 * cells[ getCellIndex(x,y,z+1) ]
+ 6.0 * cells[ getCellIndex(x,y,z) ]
);
globalTemp += temp * temp;
cells[ getCellIndex(x,y,z) ] += -omega * temp / 6.0 * CONSTANT * CONSTANT;
}
【问题讨论】：
-
小贴士：而不是
ix + (iy * numCellsX) + (iz * numCellsX * numCellsY)，使用链式计算ix + numCellsX*(iy + iz * numCellsY)，少一个*。 -
如果所有块都不是 4x4x4，那么你如何确定哪些不是？您是分配一个精确大小的数组
X*Y*Z还是将其四舍五入以便最终得到精确的(X*Y*Z) / (4*4*4)立方体？我认为您需要先弄清楚这一点，然后才能回答问题。 -
@Lundin 可以填充不是 4x4x4 的块以使其成为 4x4x4，但是当我遍历每个块时，我需要知道哪些单元格是填充单元格，这需要一个 if 语句停止矢量化。任何想法如何解决这个问题？
-
我假设这些设置类似于
int cells[numCellsZ][numCellsY][numCellsX]... 在这种情况下，您不能通过执行某些操作从单元格[0][0][0] 中获取偏移量喜欢(int *) &cells[z][y][x] - (int *) &cells[0][0][0] -
@technosaurus 这些单元格被声明为单元格[numCells]，其中 numCells = (numX+2) + (numY+2) + (numZ+2)。请参阅布局的索引功能。你能解释一下你的解决方案吗？我认为它只会给我一个单元格的全局索引，这就是 get index 函数的作用。
标签： c arrays optimization multidimensional-array tiling