c - 如何将数组分成 block

Question

我有一个代表长方体中点的数组。它是一个一维数组，使用如下索引函数实现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

对于每个单元格，我想根据它的邻居值(即模板)为该单元格计算一个新值，它的邻居位于上方、下方、左侧、右侧、前方和后方。但是，我只想对还不错的单元格进行此计算。我有一个 bool 数组，用于跟踪单元格是否损坏。这是当前计算的样子:

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 语句我不能。我知道细胞是否提前坏了，这在整个计算过程中都不会改变。我想将域分成 block ，最好是 4x4x4 block ，这样我可以计算每个 block 的先验，如果它包含坏单元格，如果这样处理它，或者如果没有，使用可以采取的优化函数矢量化的优势，例如

for(block : blocks) {
    if(isBadBlock[block]) {
        slowProcessBlock(block) // As above
    } else {
        fastVectorizedProcessBlock(block)
    }
}

注意: block 不需要物理存在，即这可以通过更改索引函数并使用不同的索引遍历数组来实现。我愿意接受最有效的方法。

fastVectorizedProcessBlock() 函数看起来类似于 slowProcessBlock() 函数，但带有 if 语句 remove(因为我们知道它不包含坏单元格)和矢量化编译指示。

如何将我的域拆分成多个 block 以便完成此操作？这看起来很棘手，因为 a) 每个方向的单元格数量不相等，b) 我们需要考虑填充单元格，因为我们绝不能尝试计算它们的值，因为这会导致内存访问失败的界限。

然后如何在不使用 if 语句的情况下处理不包含坏单元格的 block ？

编辑:

这是我最初的想法:

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
                }
             }
         }
     }
 }

单元格现在将存储在 block 中，即第一个 4x4x4 block 中的所有单元格将存储在 pos 0-63 中，然后第二个 block 中的所有单元格将存储在 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;
}

Answer 1

getCellIndex() 在哪里检索 numCellX 和 numCellY 的值？最好将它们作为参数传递而不是依赖全局变量，并使此函数static inline 以允许编译器优化。

static line int getCellIndex(int ix, int iy, int iz, int numCellsX, numCellsY) {
    return ix + (iy * numCellsX) + (iz * numCellsX * numCellsY);
}

for (int z = 1; z <= numZ; z++) {
    for (int y = 1; y <= numY; y++) {
        for (int x = 1; x <= numX; x++) {
            if (!isBadCell[getCellIndex(x, y, z, numX + 2, numY + 2)] {
                // Do stencil Computation
            }
        }
    }
}

您还可以删除所有与一些局部变量的乘法:

int index = (numY + 2) * (numX + 2);  // skip top padding plane
for (int z = 1; z <= numZ; z++) {
    index += numX + 2;  // skip first padding row
    for (int y = 1; y <= numY; y++) {
        index += 1;   // skip first padding col
        for (int x = 1; x <= numX; x++, index++) {
            if (!isBadCell[index] {
                // Do stencil Computation
            }
        }
        index += 1;   // skip last padding col
    }
    index += numX + 2;   // skip last padding row
}

这些方向是否可行在很大程度上取决于为获得模板值而执行的实际计算。你也应该发布它。

如果您可以更改坏单元格的 bool 数组格式，将行填充为 8 的倍数并使用 8 列的水平填充来改进对齐将很有用。使 bool 数组成为位数组允许一次检查 8、16、32 甚至 64 个单元格。

您可以调整数组指针以使用基于 0 的坐标。

这是它的工作原理:

int numCellsX = 8 + ((numX + 7) & ~7) + 8;
int numCellsY = 1 + numY + 1;
int numCellsXY = numCellsX * numCellsY;
// adjusted array_pointer
array_pointer = allocated_pointer + 8 + numCellsX + numCellsXY;
// assuming the isBadCell array is 0 based too.
for (int z = 0, indexZ = 0; z < numZ; z++, indexZ += numCellsXY) {
    for (int y = 0, indexY = indexZ; y < numY; y++, indexY += numCellsX) {
        for (int x = 0, index = indexY; x <= numX - 8; x += 8, index += 8) {
            int mask = isBadCell[index >> 3];
            if (mask == 0) {
                // let the compiler unroll computation for 8 pixels with
                for (int i = 0; i < 8; i++) {
                   // compute stencil value for x+i,y,z at index+i
                }
            } else {
                for (int i = 0; i < 8; i++, mask >>= 1) {
                    if (!(mask & 1)) {
                       // compute stencil value for x+i,y,z at index+i
                    }
                }
            }
        }
        int mask = isBadCell[index >> 3];
        for (; x < numX; x++, index++, mask >>= 1) {
            if (!(mask & 1)) {
                // compute stencil value for x,y,z at index
            }
        }
    }
}

编辑:

模板函数对 getCellIndex 的调用过多。以下是如何使用上面代码中计算的索引值对其进行优化:

// index is the offset of cell x,y,z
// numCellsX, numCellsY are the dimensions of the plane
// numCellsXY is the offset between planes: numCellsX * numCellsY

if (isBadCell[index]) {
    double temp = someOtherArray[index] +
                1.0 / CONSTANT / CONSTANT *
                ( - 1.0 * cells[index - 1]
                  - 1.0 * cells[index + 1]
                  - 1.0 * cells[index - numCellsX]
                  - 1.0 * cells[index + numCellsX]
                  - 1.0 * cells[index - numCellsXY]
                  - 1.0 * cells[index + numCellsXY]
                  + 6.0 * cells[index]
                );
    cells[index] += -omega * temp / 6.0 * CONSTANT * CONSTANT;
    globalTemp += temp * temp;
}

预先计算 &cells[index] 作为指针可能会改进代码，但编译应该能够检测到这个公共(public)子表达式并生成高效代码。

编辑 2:

这是一种平铺方法:您可以添加缺少的参数，大多数大小都假定为全局大小，但您可能应该传递一个指向包含所有这些值的上下文结构的指针。它使用 isBadTile[] 和 isGoodTile[]: bool 数组分别判断给定的图 block 是否包含所有单元格都坏和所有单元格都好。

void handle_tile(int x, int y, int z, int nx, int ny, int nz) {
    int index0 = x + y * numCellsX + z * numCellsXY;
    // skipping a tile with all cells bad.
    if (isBadTile[index0] && nx == 4 && ny == 4 && nz == 4)
        return;
    // handling a 4x4x4 tile with all cells OK.
    if (isGoodTile[index0] && nx == 4 && ny == 4 && nz == 4) {
        for (int iz = 0; iz < 4; iz++) {
            for (int iy = 0; iy < 4; iy++) {
                for (int ix = 0; ix < 4; ix++) {
                    int index = index0 + ix + iy * numCellsX + iz + numCellsXY;
                    // Do stencil computation using `index`
                }
            }
        }
    } else {
        for (int iz = 0; iz < nz; iz++) {
            for (int iy = 0; iy < ny; iy++) {
                for (int ix = 0; ix < nx; ix++) {
                    int index = index0 + ix + iy * numCellsX + iz + numCellsXY;
                    if (!isBadCell[index] {
                        // Do stencil computation using `index`
                }
            }
        }
    }
}

void handle_cells() {
    int x, y, z;
    for (z = 1; z <= numZ; z += 4) {
        int nz = min(numZ + 1 - z, 4);
        for (y = 1; y <= numY; y += 4) {
            int ny = min(numY + 1 - y, 4);
            for (x = 1; x <= numX; x += 4) {
                int nx = min(numX + 1 - x, 4);
                handle_tile(x, y, z, nx, ny, nz);
            }
        }
    }
}

这是一个计算 isGoodTile[] 数组的函数。唯一正确计算的偏移对应于 x 的倍数 4 + 1，y 和 z 的最大值小于 3。

这个实现不是最优的，因为可以计算的元素更少。不完整的边界 block (距离边缘少于 4 个)可能被标记为不好，以跳过单个案例的好案例。如果为边缘图 block 正确计算了 isBadTile 数组，则对坏图 block 的测试可能适用于这些边缘图 block ，但目前情况并非如此。

void computeGoodTiles() {
    int start = 1 + numCellsX + numCellsXY;
    int stop = numCellsXY * numCellsZ - 1 - numCellsX - numCellsXY;

    memset(isGoodTile, 0, sizeof(*isGoodTile) * numCellsXY * numCellsZ);
    for (int i = start; i < stop; i += 4) {
        isGoodTile[i] = (isBadCell[i + 0] | isBadCell[i + 1] |
                         isBadCell[i + 2] | isBadCell[i + 3]) ^ 1;
    }
    for (int i = start; i < stop - 3 * numCellsX; i += 4) {
        isGoodTile[i] = isGoodTile[i + 0 * numCellsX] &
                        isGoodTile[i + 1 * numCellsX] &
                        isGoodTile[i + 2 * numCellsX] &
                        isGoodTile[i + 3 * numCellsX];
    }
    for (int i = start; i < stop - 3 * numCellsXY; i += 4) {
        isGoodTile[i] = isGoodTile[i + 0 * numCellsXY] &
                        isGoodTile[i + 1 * numCellsXY] &
                        isGoodTile[i + 2 * numCellsXY] &
                        isGoodTile[i + 3 * numCellsXY];
    }
}

void computeBadTiles() {
    int start = 1 + numCellsX + numCellsXY;
    int stop = numCellsXY * numCellsZ - 1 - numCellsX - numCellsXY;

    memset(isBadTile, 0, sizeof(*isBadTile) * numCellsXY * numCellsZ);
    for (int i = start; i < stop; i += 4) {
        isBadTile[i] = isBadCell[i + 0] & isBadCell[i + 1] &
                       isBadCell[i + 2] & isBadCell[i + 3];
    }
    for (int i = start; i < stop - 3 * numCellsX; i += 4) {
        isBadTile[i] = isBadTile[i + 0 * numCellsX] &
                       isBadTile[i + 1 * numCellsX] &
                       isBadTile[i + 2 * numCellsX] &
                       isBadTile[i + 3 * numCellsX];
    }
    for (int i = start; i < stop - 3 * numCellsXY; i += 4) {
        isBadTile[i] = isBadTile[i + 0 * numCellsXY] &
                       isBadTile[i + 1 * numCellsXY] &
                       isBadTile[i + 2 * numCellsXY] &
                       isBadTile[i + 3 * numCellsXY];
    }
}