Уменьшение количества инструкций, удаленных после развертывания цикла

Непонятно, где gcc будет сокращать количество инструкций. Возможно, что повышенное давление на регистры может побудить gcc использовать инструкции load+operate (то есть такое же количество примитивных операций, но меньше инструкций). Индекс для f_L будет увеличиваться только один раз за самый внутренний цикл, но это сэкономит только 6,2 ГБ (3 * 64 * 256 * 256 * 496) инструкций. (Кстати, накладные расходы цикла должны составлять только три инструкции, поскольку i должно оставаться в регистре.)

Следующая псевдосборка (для RISC-подобной ISA) с использованием двухстороннего развертывания показывает, как можно сохранить приращение:

// the address of f_L[k * L * L + j * L + i] is in r1
// (float)(1.0 / (w_n * w_n) * pix_1) results are in f1 and f2
load-single f9 [r1];    // load float at address in r1 to register f9
add-single f9 f9 f1;    // f9 = f9 + f1
store-single [r1] f9;   // store float in f9 to address in r1
load-single f10 4[r1];  // load float at address of r1+4 to f10
add-single f10 f10 f2;  // f10 = f10 + f2
store-single 4[r1] f10; // store float in f10 to address of r1+4
add r1 r1 #8;           // increase the address by 8 bytes

Трассировка двух итераций неразвернутой версии будет выглядеть примерно так:

load-single f9 [r1];  // load float at address of r1 to f9
add-single f9 f9 f2;  // f9 = f9 + f2
store-single [r1] f9; // store float in f9 to address of r1
add r1 r1 #4;         // increase the address by 4 bytes
...
load-single f9 [r1];  // load float at address of r1 to f9
add-single f9 f9 f2;  // f9 = f9 + f2
store-single [r1] f9; // store float in f9 to address of r1
add r1 r1 #4;         // increase the address by 4 bytes

Поскольку инструкции по адресации памяти обычно включают добавление непосредственного смещения (Itanium является необычным исключением), а конвейеры обычно не реализуются для оптимизации случая, когда непосредственное значение равно нулю, использование ненулевого непосредственного смещения обычно «бесплатно». (Конечно, это уменьшает количество инструкций — 7 вместо 8 в данном случае — но в целом также повышает производительность.)

Что касается прогнозирования ветвлений, согласно Микроархитектура процессоров Intel, AMD и VIA Агнера Фога: руководство по оптимизации для программистов на ассемблере и разработчиков компиляторов(PDF) предиктор ветвлений микроархитектуры Core2 использует 8-битную глобальную историю. Это означает, что он отслеживает результаты последних 8 переходов и использует эти 8 бит (вместе с битами из адреса инструкции) для индексации таблицы. Это позволяет распознавать корреляции между соседними ветвями.

Для вашего кода ветвь, соответствующая, например, 8-й предыдущей ветви, не является одной и той же ветвью на каждой итерации (поскольку используется короткое замыкание), поэтому непросто понять, насколько хорошо будут распознаваться корреляции.

Некоторые корреляции в ветвях очевидны. Если px_x[0]>=0 истинно, то px_x[0]+1>=0 также будет истинным. Если px_x[0] <(int)threadCopy[0].S_x верно, то px_x[0]+1 <(int)threadCopy[0].S_x скорее всего будет правдой.

Если развертывание выполнено таким образом, что px_x[n] проверяется для всех четырех значений n, тогда эти корреляции будут отодвигаться дальше, чтобы результаты не использовались предсказателем ветвления.

Некоторые возможности оптимизации

Хотя вы не спрашивали о каких-либо возможностях оптимизации, я собираюсь предложить некоторые возможности для исследования.

Во-первых, для ветвей, если не подходит строгая переносимость, тест x>=0 && x<y можно упростить до (unsigned)x<(unsigned)y. Это не является строго переносимым, потому что, например, машина теоретически может представлять отрицательные числа в формате знака-величины со старшим значащим битом в качестве знака и отрицательным значением, обозначаемым нулевым битом. Для обычных представлений целых чисел со знаком такое повторное приведение будет работать до тех пор, пока y является положительным целым числом со знаком, поскольку отрицательное значение x будет иметь установленный старший бит и, следовательно, будет больше, чем y, интерпретируемое как целое число без знака.

Во-вторых, количество ветвей можно значительно уменьшить, используя 100% корреляции для px_x или px_y:

if ((unsigned) px_y[0]<(unsigned int)threadCopy[0].S_y)
{
    if ((unsigned)px_x[0]<(unsigned int)threadCopy[0].S_x)
        pixel_1[0] = threadCopy[0].I_n[px_y[0] * threadCopy[0].S_x + px_x[0]];
    else
        pixel_1[0] = 0.0;
    if ((unsigned)px_x[0]+1<(unsigned int)threadCopy[0].S_x)
        pixel_1[2] = threadCopy[0].I_n[px_y[0] * threadCopy[0].S_x + (px_x[0]+1)];
    else
        pixel_1[2] = 0.0;
}
if ((unsigned)px_y[0]+1<(unsigned int)threadCopy[0].S_y)
{
    if ((unsigned)px_x[0]<(unsigned int)threadCopy[0].S_x)
        pixel_1[1] = threadCopy[0].I_n[(px_y[0]+1) * threadCopy[0].S_x + px_x[0]];
    else
        pixel_1[1] = 0.0;
    if ((unsigned)px_x[0]+1<(unsigned int)threadCopy[0].S_x)
        pixel_1[3] = threadCopy[0].I_n[(px_y[0]+1) * threadCopy[0].S_x + (px_x[0]+1)];
    else
        pixel_1[3] = 0.0;
}

(Если приведенный выше фрагмент кода реплицируется для развертывания, его, вероятно, следует реплицировать как блок, а не чередовать тесты для разных значений px_x и px_y, чтобы ветвь px_y находилась рядом с ветвью px_y+1, а первая ветвь px_x была рядом с ветвью px_x. другой филиал px_x и филиалы px_x+1.)

Другая возможная оптимизация — изменение вычисления w_n на вычисление его обратной величины. Это изменит умножение и три деления на четыре умножения и одно деление. Деление намного дороже умножения. Кроме того, вычисление приблизительной обратной величины намного более удобно для SIMD, поскольку обычно существуют инструкции по обратной оценке, которые обеспечивают отправную точку, которую можно уточнить с помощью метода Ньютона-Рафсона.

Если приемлема даже худшая обфускация кода, вы можете рассмотреть возможность изменения кода, такого как double y = O_L + (double)j * R_L;, на double y = O_L; ... y += R_L;. (Я провел тест, и gcc, похоже, не распознал эту оптимизацию, вероятно, из-за использования плавающей запятой и приведения к двойному.) Таким образом:

for(int imageNo =0; imageNo<496;imageNo++){

double z = O_L;
for (unsigned int k=0; k<256; k++)
{

    double y = O_L;
    for (unsigned int j=0; j<256; j++)
    {
        double x[1]; x[0] = O_L;
        for (unsigned int i=0; i<256; i++)
        {
            ...
            x[0] +=  R_L ;
        } // end of i loop
        y += R_L;
    }  // end of j loop
    z += R_L;
} // end of k loop
    } // end of imageNo loop

Я предполагаю, что это лишь незначительно улучшит производительность, поэтому стоимость запутывания будет выше по сравнению с выгодой.

Еще одно изменение, которое, возможно, стоит попробовать, — это включение некоторых вычислений pix_1 в разделы, устанавливающие условные значения pixel_1[]. Это значительно запутало бы код и не принесло бы большой пользы. Кроме того, это может затруднить автовекторизацию компилятором. (При условной установке значений на соответствующий I_n или на ноль сравнение SIMD может установить каждый элемент на -1 или 0, а простое and со значением I_n даст правильное значение. В этом случае накладные расходы на формирование I_n vector, вероятно, не будет иметь смысла, учитывая, что Core2 поддерживает только SIMD двойной точности с двойной точностью, но с поддержкой сбора или даже просто с более длинным вектором компромиссы могут измениться.)

Однако это изменение увеличило бы размер базовых блоков и уменьшило объем вычислений, когда какие-либо из px_x и px_y выходят за пределы допустимого диапазона (я предполагаю, что это редкость, поэтому преимущество будет очень небольшим при Лучший).

double pix_1 = 0.0;
double alpha_diff = 1.0 - alpha;
if ((unsigned) px_y[0]<(unsigned int)threadCopy[0].S_y)
{
    double beta_diff = 1.0 - beta;
    if ((unsigned)px_x[0]<(unsigned int)threadCopy[0].S_x)
        pix1 += alpha_diff * beta_diff 
             * threadCopy[0].I_n[px_y[0] * threadCopy[0].S_x + px_x[0]];
    // no need for else statement since pix1 is already zeroed and not 
    // adding the pixel_1[0] factor is the same as zeroing pixel_1[0]
    if ((unsigned)px_x[0]+1<(unsigned int)threadCopy[0].S_x)
        pix1 += alpha * beta_diff 
             * threadCopy[0].I_n[px_y[0] * threadCopy[0].S_x + (px_x[0]+1)];
}
if ((unsigned)px_y[0]+1<(unsigned int)threadCopy[0].S_y)
{
    if ((unsigned)px_x[0]<(unsigned int)threadCopy[0].S_x)
        pix1 += alpha_diff * beta 
             * threadCopy[0].I_n[(px_y[0]+1) * threadCopy[0].S_x + px_x[0]];
    if ((unsigned)px_x[0]+1<(unsigned int)threadCopy[0].S_x)
        pix1 += alpha * beta 
             * threadCopy[0].I_n[(px_y[0]+1) * threadCopy[0].S_x + (px_x[0]+1)];
}

В идеале такой код, как ваш, должен быть векторизован, но я не знаю, как заставить gcc распознавать возможности, как выражать возможности с помощью встроенных функций, а также стоит ли прилагать значительные усилия по ручной векторизации этого кода при ширине SIMD, равной всего двум. .

Я не программист (просто тот, кто любит изучать и думать об архитектуре компьютера), и у меня есть значительная склонность к микрооптимизации (как ясно из вышеизложенного), поэтому приведенные выше предложения следует рассматривать в этом свете.