Продолжаю оптимизировать

Точнее, пытаюсь закончить оптимизировать Scalar SSE в тех местах, которые до конца не проверил, или проверил не правильно. 

1. Замена вычитания сложением. В чем суть этой попытки оптимизации? У нас есть вот такое выражение:
  S[j] := S[j]-R[j]*S[index];
На ассемблере я его транслировал в такой набор команд
  movapd XMM1, XMM0;
  mulsd XMM1, [R11+RDX*8];
  movsd XMM2, [RCX+RDX*8];
  subsd XMM2, XMM1;
  movsd [RCX+RDX*8], XMM2;
в двух операторной форме команд, принятой в x86, я вынужден сначала загрузить S[j] в регистр, а уже потом выполнять вычитание.
Если же у S[index] перед внутренним циклом поменять знак, то операцию вычитания можно заменить сложением без предварительной загрузки из памяти в регистр XMM2. В результате получится вот такой код:
  movapd XMM1, XMM0;
  mulsd XMM1, [R11+RDX*8];
  addsd XMM1, [RCX+RDX*8];
  movsd [RCX+RDX*8], XMM1;
Как видим, он на одну команду короче, и, возможно, быстрее. Не стоит забывать, что нам пришлось во внешней цикле добавить команду изменения знака S[index].
Но будет ли этот код реально быстрее? Зависит от того, какая команда выполняется быстрее: умножение XMM1 на значение в память или загрузка в XMM2 значения из памяти. Дело в том, что команда addsd не может начаться, пока не закончится умножение XMM1.
Поэтому, если умножение медленное, то XMM2 успеет загрузиться из памяти (так как эти команды независимы) и затем можно будет вычесть из XMM2 XMM1 без потери производительности. То есть в этом случае замена вычитания сложением не нужна, так как ничего не даст.
В том случае, если умножение выполняется быстрее, чем загрузка XMM2, то ситуация будет обратная: выгоднее будет второй вариант со сложением.
На мой взгляд, разницы тут не должно быть много, но давайте посмотрим на результаты.
На типе single (float) ускорение от замены вычитания на сложение составило от сотых долей процента на системах из 3 уравнений, до 4-5% на больших системах до 48 переменных включительно. А вот на больших размерностях результат оказался парадоксальный: незначительное снижение производительности на системах из 96 уравнений и потери до 3.5% и выше на больших размерностях.
На типе double ускорение не превышает 3%, достигая максимума на 24 уравнениях, затем постепенно снижается и заменяется замедлением на системах в 1536 уравнений.
В общем, я не увидел особого смысла в такой оптимизации.

2. Развертывание цикла. Я прошлый раз ограничился четырьмя операциями, но Иван предложил попробовать увеличить их до восьми.
Результат моих размышлений на эту тему.
Ускорение от развертывания можно получить за счет снижения доли служебных операций на организацию цикла. Но это работает только в том случае, когда тело цикла (рабочая нагрузка) очень проста и по времени сравнима с работой по организации самого цикла.
В моем случае сам цикл организуется всего тремя командами: двумя короткими регистровыми операциями и условным переходом. А полезной нагрузкой является 5 команд, в том числе три с памятью. В принципе, доля служебных команд достаточна велика и такой цикл имеет смысл разложить на две итерации. Дальнейшее разложение приведет лишь к незначительному росту производительности, тем более, что хорошо работает блок предсказаний переходов, команды организации цикла и полезной нагрузки исполняются на разных блоках, так что частично могут выполняться параллельно.
Второй источник ускорения связан с отсутствием сбросов конвейера. Чем длиннее линейная часть программы, без условных переходов, тем полнее заполнен конвейер и выше производительность процессора. Наличие условных переходов сбрасывает конвейер,что приводит к сильному снижению производительности.
Вышеописанная страшилка более актуальная для очень старых процессоров. Современные версии процессоров имеют очень хороший блок предсказания переходов, который работает совместно с процессором и остановка конвейера происходит только в том случае, если блок предсказаний ошибся. В частности, для for-циклов блок предсказаний ошибается очень редко, в идеальном случае лишь один раз ‒ когда цикл заканчивается. Думаю, что именно этот источник ускорения при решении моей задачи имеет очень низкое значение.
Третий источник ускорения ‒ суперскалярная архитектура, причем именно она, на мой взгляд, дает наибольший эффект в моем случае. Посмотрим на полезную нагрузку цикла. В случае одной операции на итерацию у нас требуется выполнить одно умножение и одно сложение, причем сложение зависит от умножение. Фактически здесь нет операций, которые могли бы выполнятся параллельно.
Если выполнять две операции за итерацию, то нужно выполнить 2 умножения и 2 сложения, тут уже есть что распараллелить. 2 умножения можно выполнять параллельно, но со сложениям из-за зависимости от умножений ситуация сложнее, скорее всего параллельно умножениям удастся выполнить лишь одну команду сложения.
При увеличении количества операций на итерацию ситуация с распараллеливанием становится все лучше и лучше. Теоретически могут выполняться 4 инструкции за такт, но практически, есть подозрение, их количество лежит в пределах от 2 до 3, причем ближе к двум, по всей видимости. Это связано с тем, что умножение все же выполняется медленнее сложения, поэтому для сумматоров данные всегда поставляются с задержкой.
Тем не менее, попробуем оценить, насколько увеличивается производительность при переходе от 4 операций на итерацию к 8.
На типе single (float) ускорение составило от 0% (3 переменных) до 11% (192 переменных), затем уменьшилось до 4% на 384 переменных.
На типе double ускорение максимум составило около 7% на 192 переменных, на всех остальных размерах ускорение составляет около 1%, что ниже погрешности измерений.
В результате решил оставить развертывание на восемь операций для типа single (float), а для типа double остановиться на четырех.

3. Решил побаловаться с командами PREFETCH. Эффект для моей задачи равен нулю, так как, видимо, механизм автоматической предзагрузки кэша в этом случае справляется на ура.