У меня есть массив, который представляет точки в кубоиде. Это одномерный массив, который использует следующую функцию индексирования для реализации трех измерений:
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 remove (поскольку мы знаем, что он не содержит плохих ячеек) и прагмой векторизации.
Как я могу разделить свой домен на блоки, чтобы я мог это сделать? Кажется сложным, потому что: а) количество ячеек в каждом направлении не равно, б) нам нужно учитывать ячейки заполнения, так как мы никогда не должны пытаться вычислить их значение, так как это приведет к отсутствию доступа к памяти границ.
Как я могу обработать блоки, которые не содержат плохие ячейки, не используя оператор if?
EDIT:
Это идея, которую я изначально имел:
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 - Вот как выглядит трафаретное вычисление:
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;
}
