我有一个表示长方体中点的数组。它是一个一维数组，使用以下索引函数来实现三个维度：
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语句（因为我们知道它不包含坏单元），并添加了一个向量化预处理器指令。
如何将我的域分成块以便我可以完成这个任务？这似乎很棘手，因为 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块中的所有细胞将存储在位置0-63，然后第二个块中的所有细胞将存储在位置64-127等。
然而，如果numX/numY/numZ值不是很好，我认为这种方法行不通。例如，如果numY = 2，numZ = 1，numX = 10怎么办？for循环将期望z方向至少深入4个单元格。有没有好的方法来解决这个问题？
更新2 - 这是stencil计算的样子：
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 + numCellsX*(iy + iz * numCellsY)代替ix + (iy * numCellsX) + (iz * numCellsX * numCellsY)，这样可以少用一个乘号。 - chux - Reinstate Monica