Некоторая процедура рисования Мандельброта от c до sse2

Я хочу переписать такую ​​простую подпрограмму в код SSE2 (желательно в nasm), и я не совсем уверен, как это сделать, две неясные вещи (как выразить вычисления (внутренний цикл и те из внешнего цикла) и как вызвать c функция кода "SetPixelInDibInt(i,j, palette[n]);" из-под статически связанного кода asm

    void DrawMandelbrotD(double ox, double oy, double lx, int N_ITER)
    {
     double ly = lx * double(CLIENT_Y)/double(CLIENT_X);
     double dx = lx / CLIENT_X;
     double dy = ly / CLIENT_Y;
     double ax = ox - lx * 0.5 + dx * 0.5;
     double ay = oy - ly * 0.5 + dy * 0.5;
    static  double re, im, re_n, im_n, c_re, c_im, rere, imim, int n;

    for(int j=0; j<CLIENT_Y; j+=1)
    {
     for(int i=0; i<CLIENT_X; i+=1)
     {
      c_re = ax + i * dx;
      c_im = ay + j * dy;
      re = c_re;
      im = c_im;
      rere=re*re;
      imim=im*im;
      n=1;

      for(int k=0;k<N_ITER;k++)
      {
        im =  (re+re)*im    + c_im;
        re =   rere - imim  + c_re;
        rere=re*re;
        imim=im*im;
        if ( (rere + imim) > 4.0 ) break;
        n++;
       }
        SetPixelInDibInt(i ,j, palette[n]);
      }
     }
    }

Может ли кто-то помочь, я хотел бы не видеть других реализаций кода, а просто nasm-sse-перевод этих выше - это было бы наиболее полезно в моем случае - кто-то может помочь с этим?

1 ответ

Intel имеет полную реализацию в качестве примера AVX. Увидеть ниже.

Что делает Мандельброта хитрым, так это то, что раннее условие для каждой точки в наборе (т. Е. Пикселя) отличается. Вы можете оставить пару или четверку пикселей итерацией, пока величина обоих не превысит 2,0 (или вы достигнете максимума итераций). Чтобы поступить иначе, потребуется отслеживать, какие точки пикселя были в каком элементе вектора.

В любом случае, упрощенная реализация для работы с вектором 2 (или 4 с AVX) удваивается за раз, а ее пропускная способность ограничивается задержкой цепочек зависимостей. Вам нужно было бы сделать несколько цепочек зависимостей параллельно, чтобы поддерживать оба модуля FMA в Haswell. Таким образом, вы будете дублировать свои переменные и чередовать операции для двух итераций внешнего цикла внутри внутреннего цикла.

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

Я подозреваю, что проверка состояния выхода из цикла каждые 4 итерации может быть выигрышной. Получение кода для ветвления, основанного на сравнении упакованных векторов, немного дороже, чем в скалярном случае. Требуется дополнительное добавление FP, это также дорого. (Haswell может выполнять два FMA за цикл (задержка = 5). Одиночный модуль добавления FP является одним и тем же портом, что и один из модулей FMA. Два модуля FP mul находятся на тех же портах, которые могут работать FMA.)

Состояние цикла можно проверить с помощью упакованного сравнения, чтобы сгенерировать маску нулей и единиц, а также (V)PTEST этого регистра с самим собой, чтобы увидеть, если это все ноль. (редактировать: movmskps затем test+jcc меньше мопов, но, возможно, более высокая задержка.) Тогда, очевидно, je или же jne в зависимости от обстоятельств, в зависимости от того, выполняли ли вы сравнение FP, которое оставляет нули, когда вы должны выйти, или нули, когда вы не должны. NAN не должен быть возможен, но нет причин не выбирать опцию сравнения так, чтобы NAN приводил к выполнению условия выхода.

const __mm256d const_four = _mm256_set1_pd(4.0);  // outside the loop

__m256i cmp_result = _mm256_cmp_pd(mag_squared, const_four, _CMP_LE_OQ);  // vcmppd.  result is non-zero if at least one element < 4.0
if (_mm256_testz_si256(cmp_result, cmp_result))
    break;

Там может быть какой-то способ использовать PTEST непосредственно в упакованном двойнике, с некоторой битовой взломанной AND-маской, которая выберет биты, которые будут установлены, если значение FP> 4.0. Как, может быть, некоторые биты в показателе степени? Возможно стоит задуматься. Я нашел сообщение на форуме об этом, но не попробовал это.

Хм, дерьмо, это не записывает, КОГДА условие петли не удалось, для каждого элемента вектора отдельно, с целью раскраски точек за пределами набора Мандельброта. Может быть, проверить любой элемент, попадающий в условие (вместо всех), записать результат, а затем установить этот элемент (и c для этого элемента) до 0,0, поэтому он не будет вызывать условие выхода снова. Или, может быть, планирование пикселей в векторные элементы - это путь, который все-таки нужен. Этот код вполне может работать на многопоточном ЦП, так как будет много ошибочных прогнозов ветвлений, когда каждый элемент отдельно инициирует условие раннего завершения.

Это может потратить большую часть вашей пропускной способности, и учитывая, что выполнимо 4 мопа за цикл, но только 2 из них могут быть FP mul/add/FMA, есть место для значительного количества целочисленного кода для планирования точек в векторные элементы. (В Sandybridge/Ivybrideg без FMA пропускная способность FP ниже. Но есть только 3 порта, которые могут обрабатывать целочисленные операции, и 2 из них являются портами для модулей FP mul и FP add.)

Поскольку вам не нужно читать какие-либо исходные данные, для каждой цепочки dep есть только 1 поток доступа к памяти, и это поток записи. (И это низкая пропускная способность, так как большинству точек требуется много итераций, прежде чем вы будете готовы записать однопиксельное значение.) Таким образом, количество потоков аппаратной предварительной выборки не является ограничивающим фактором для количества параллельных цепочек деп для выполнения., Кэширование задержек должно быть скрыто буферами записи.

Я могу написать некоторый код, если кто-то все еще заинтересован в этом (просто оставьте комментарий). Я остановился на этапе проектирования высокого уровня, так как это старый вопрос.

==============

Я также обнаружил, что Intel уже использовала набор Мандельброта в качестве примера для одного из своих учебных пособий по AVX. Они используют метод mask-off-vector-elements для условия цикла. (используя маску, сгенерированную непосредственно vcmpps в И). Их результаты показывают, что AVX (одинарная точность) обеспечил 7-кратное ускорение по сравнению со скалярным плавающим числом, поэтому очевидно, что соседние пиксели не достигают условия раннего выхода при очень различном количестве итераций. (по крайней мере, для увеличения / панорамирования, с которыми они тестировали.)

Они просто позволяют накапливать результаты FP для элементов, которые не соответствуют условию раннего запуска. Они просто прекращают увеличивать счетчик для этого элемента. Надеемся, что большинство систем по умолчанию установят управляющее слово для обнуления денормали, если денормали все еще требуют дополнительных циклов.

Однако их код глуп с одной стороны: они отслеживают количество итераций для каждого элемента вектора с вектором с плавающей запятой, а затем преобразуют его в int в конце перед использованием. Было бы быстрее, и не занимать исполняющую единицу FP, использовать для этого упакованные целые числа. О, я знаю, почему они это делают: AVX (без AVX2) не поддерживает 256-битные целочисленные векторные операции. Они могли бы использовать упакованные 16-битные счетчики цикла int, но это могло бы привести к переполнению. (И им придется сжать маску с 256b до 128b).

Они также тестируют для всех элементов, являющихся> 4.0 с movmskps а затем проверить это, вместо того, чтобы использовать ptest, Я думаю, test / jcc может слиться с макросами и работать на другом исполнительном модуле, чем векторные операции FP, так что, возможно, это даже не медленнее. О, и, конечно, AVX (без AVX2) не имеет 256 бит PTEST, Также, PTEST это 2 моп, так что на самом деле movmskps + test / jcc меньше мопсов, чем ptest + jcc, (PTEST 1 Snop-домен UOP на SnB, но все еще 2 неиспользованных UOP для портов выполнения. На IvB/HSW 2 мопа даже в слитой области.) Так это выглядит movmskps является оптимальным способом, если вы не можете воспользоваться побитовым И это является частью PTEST или нужно проверять не только старший бит каждого элемента. Если ветка непредсказуема, ptest может быть меньше задержки, и, следовательно, стоит того, чтобы поймать ошибочно предсказывает цикл раньше.

Другие вопросы по тегам