Assembler VS C/C++

use FORTRAN

На ассемблере имхо лучше, только придётся под все target архитектуры писать код. Впрочем это, имхо, не большая проблема, программа обычно работает на одной архитектуре.

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

А пацаны-то с http://www.netlib.org и не знают :{

Системы моделирования на Fortranе пишут только извращенцы, тут вопрос о грамотной реализации пакетных SIMD расширений, которые в процах кстати не зря появились, и увы из "обычного" языка будь то Pascal, Fortran, C, C++, ADA НЕПРОИЗРАСТАЮТ ибо компилеры ДАЖЕ интеловский ПАКЕТНЫЕ инструкции генерировать НЕУМЕЮТ, и уж если и писать на языке высокого уровня то нужно пользоваться intristicами а для них простите нативным языком является С и заниматься извращенствами с фортраном уже слишком, когда стоит вопрос о том НЕ ИЗВРАТОМ ли является ВООБЩЕ язык высокого уровня перед ассемблером. Так что харе флудить господа.

2Legioner: Помне так переносимость падает как только я решаю использовать нечто этакое что выходит за рамки общих методик программирования. Ну нет в соседних архитектурах SIMD в том виде как он есть у интеля и конечно в любом случае под другую архитектуру придется переписывать - это ежику ясно (даже чисто сишный код придется пеерписывать!)

> Системы моделирования на Fortranе пишут только извращенцы

Да что вы говорите, мужики-то и не знают :-(

Мне не очень понятно зачем вы хотите переписывать на ассемблере. Не устраивает скорость обработки? Пробовали распараллеливать вычисления на несколько процессоров, на несколько компов? Пробовали спец пакеты оптимизированные для матричных вычислений?

> вопрос о том НЕ ИЗВРАТОМ ли является ВООБЩЕ язык высокого уровня перед ассемблером

Как по мне, так вопрос в том не является ли спешное переписывание задачи на ассемблере извратом, когда не все высокоуровневые возможности испробованы.

Жалко когда другие судят о том что ты пробывал а что нет... Знаете аноним когда N=10000 то хвататься будешь за каждую мелоч лишбы поделить коичество операций на цифру побольше. Какая блин быгода от нескольких процессов когда аппаратно вычислительное ядро ОДНО (ну макчимум 4, ну совсем максимум 8) Больше процессов чем аппаратных вычислителей делать бесполезно, а вот обсчитывать 2 или 4 потока данных на одном вычислительном ядре очень даже приятно и ПРАВИЛЬНЕЙ чем тыкать пальцем в ЧУЖИЕ хоть и хорошие либы но написанные НЕ ДЛЯ ТВОЕЙ задачи. Переписывание не спешное а продуманное. Оптимизировать количество разностных вычислений на нессиметричной неэрмитовой матрице увы нельзя (факт медицинский). А нескольких компов у меня нет, как нет и сервачка с двумя четырехядерными ксеонами. А распаралелив потоки данных (если вам такой факт известен) можно много денег сэкономить... только тут не надо фортранолюбителей слушать, которые сос своим языком на перевес бегают и тыкают им куда не надо тыкать. Или фортран умеет генерить SIMD вычисления? Код и команда компиляции после которой в дизасемблере я найду addpd в студию!

Ты, собственно, от использования SIMD какого прироста ожидаешь?

>Так что лучше (с точки зрения оптимизации, простоты и удобства)?

По-моему без разницы: ты свой код уже накрепко привязал к конкретному типу процессора.

В принципе более предпочтителен второй вариант: возникают некоторое подобие API, которое можно реализовать и на ЯВУ, если возникнет жесткая необходимость. С другой стороны - лишняя межмодульная связанность и необходимость документировать интерфейсы.

Вообще-то есть какая-то низкоуровневая интеловская библиотека IPP, хотя закрытая и платная (для коммерческого использования), но цена ее не очень высокая. Может тебе подойдет?

Еще можно посмотреть в сторону генераторов кода, т.е. что-то типа DSL, хотя тут придется потратиться на разработку компилятора для целевой платформы.

Фича в том что например один addpd выполняется ровно столькоже сколько один addsd, по тем же причинам два идущих подряд (или почти подряд) addsd выполняются столькоже...

Могу пояснить - при конвеерной архитектуре на корках есть 4 потока обработчиков мопсов которые вообще говоря пахают параллельно и получается так что параллельно может абрабатываться либо одна задача сложения упакованных чисел (2 или 4 сложения зависит от размера складываемых чисел) либо параллельно столькоже эквивалентных задач сложения скалярных (неупакованных) чисел. Причина в том что они бъются в почти одинаковый набор мопсов. Под "или почти подряд" имелось ввиду допустимость наличия посредине операций загрузки выгрузки из паимяти в/из регистров. Эти операции могут быть упорядочены на ранних стадиях конвеерной обработки и распараллелены (а при наличии данных в кэше становятся вообще бесплатными).

Однако как понимаете весь этот трюкотаж работает на уровне конвеерной обработки и следовательно зыбок. Сложение упакованных чисел при вызове процедуры в цикле всегда будет быстрее двух сложений в двух последовательных шагах цикла (ибо гарантии попадания двух сложений в один конвеерный цикл уже просто нет да и не попадут они туда так как есть еще инкремент счетчика который сбивает конвеер). Итого в сухом остатке - гипотетически выигрыш должен быть 100% минус потери на упаковку/распаковку, что экспериментально дает прирост на 75-88% в зависимости от модели проца, а это согласитесь не мало (Это померено на модельной программе - куске проги реальной).

По поводу интеловских либ. Ну имея разностные выражения с простыми алгебраическими операциями глупо ожидать что они будут в сторонних либах да еще и оптимизированы под SIMD (у меня в одной задачи вариаций трехэтажных дробей 43 штуки и что все я найду в либе?). Речь идет об жесткой оптимизации по времени и памяти "школьной формулы больших размкров" а для этого либы не пишут.

ПО теме. Возникает подозрение что оптимизацию под SIMD проще делать руками чим интристиками, гарантий больше а гемор тот же.

Даже у intel в optimisation reference manual для simd архитектур почти все примеры на ассемблере - я бы даже не задумывался над выбором языка :)

из личного опыта (со-автор библиотеки по обработки изображений) ..
Мы тоже , по-возможности, везде старались использовать в коде MMX,
(сам код на C), но с некоторых пор оказалось, что приемущество его
использования стала почти нулевая (это был некий шок),
т.е. современные компиляторы умудряются так оптимизировать код,
что сравнимо с использованием SIMD (и, возможно они это делают на
уровне MMX).
Конечно, многое зависит от характера задачи, но если использовать MMX
и не терять кроссплатфор/процессоро-менность -> никуда не деться от
#ifdefов

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

> Знаете аноним когда N=10000

ну, это еще цветочки ... тут у нас задачка с N=1000000 :))

Экстраполяция нашей задачи на трехмерный случай приведет не к N^4 а к N^6 так что и n=10000 уже весело.

Про MMX не знаю, а вот SIMD оптимизацию на уровне пакетной арифметики с плавающей точкой компилеры делать не умеют...

Посмотри на BLAS библиотеки, их оптимизируют производители процессоров AMD, Indel etc, и ещё ATLAS, она сама потом выдаст оптимированную под конкретную систему BLAS.

Ну и вообще при сложности О(N^4), незнаю насколько имеет смысл оптимиривать на ниском уровне. Какого типа хоть матрицы (spars?) и что за алгоритмы? Может лучше попытатся снизить сложность?

К сожелению снизить сложность нельзя. Задача - моделирование кинетики квантовых систем. Матрицы - это результат сеточного расчета вероятностей переходов частиц в "К" пространстве при рассеении. Типичные вычисляемые величины это:
sqrt(x^2+y^2), x^2+y^2, 1/(x^2+y^2), (x+y)^2, 1/sqrt(x^2+y^2)

При этом матрицы симметрией не обладают (никакой) увы уж такова специфика задачи.

Касательно стоит-ли оптимизировать на низком уровне - ответ ДА стоит. SIMD оптимизацию на сколько я понят моно сделать только ручками в ассемблере (увы в гцц интристики работают кривовато - периодически сегволтятся, хотя может и к лучшему). Так вот при использовании float можно теоретически выйиграть в 4 раза (или на 300%) при использовании SIMD по отношению к линейному коду на С/С++/FORTRAN, при использовании double в 2 раза (или на 100%).

Так вот эксперименты по расчету одинакового количества контрольных выражений (приведены выше) показали что выигрыш при
float - +202% (+67% для AMD)
double - +53% (-37% для AMD)
В реальных вычислениях из контрольный выражений составляются трехэтажные дроби, так что эксперимент отражает реальность.

Можно также сказать что выигрыш (или проигрыш для AMD) тем больше чем сложнее выражение (при этом * и / сложнее + и - а sqrt сложнее * и /). И полагаю что 50 строк ассемблерного кода вполне стоят свечь, особенно когда вычисления идут на float (И расчеты идут на машине которую не скупердяи покупали - в смысле проц интеловский стоит).

PS: Отмечу что данный результат наглядно показывает (по крайне мере мне) какой гадостью являются АМДшные процы. Их работа никак не корреллируется с теоретическими ожиданиями, а SIMD блоки вообще работают через одно место... Единственное где они еще с интелями могут потягяться так это на операциях сложения и умножения (с вычитанием хуже) и то если есть много операций доступа к памяти (интел тут проигрывает ибо контроллер памяти внешний и частота доступа поменьше...)

конечно, у нас матрицы специфичные
(задачи из области томографии, image fidelity и т.д.)
по-поводу библиотек, которые работают с большими
матрицами (все они используют, в какой-то мерe, модифицированный BLAS):
http://trilinos.sandia.gov/
http://crd.lbl.gov/%7Exiaoye/SuperLU/#superlu?
я их смотрел, но пока не пробовал (руки не доходят).
Было бы интересно узнать о них мнение от тех, кто с ними работал.

> Про MMX не знаю, а вот SIMD
MMX - это и есть SIMD http://ru.wikipedia.org/wiki/MMX
но не "плавающей точкой"

Запор как и Белаз тоже машина. Вот вы например какой ММХ применяли? (Да их вообще говоря два, над mm и xmm регистрами). И как применяли? А потом если смотреть в доки от интеля например то MMX у них SIMDом не называется, хотя формально им и является. И блоки реализующие разную SIMD арифметику разные и обладают разной эффективностью

> вы например какой ММХ применяли?
пример кода, с коментарием :)

/* MMX implementaion of alpha-blending below turns out to be
30% slower then the original integer math implementation under it
I'm probably stupid or something.
*/
__m64 va = _mm_set_pi16 (ca, a, ca, a);
__m64 vd = _mm_set_pi16 (br[i],tr[i],ba[i],ta[i]);

/* b=(b*ca + t*a)>>8 */
vd = _mm_srli_pi16( vd, 8 );
vd = _mm_madd_pi16( va, vd );
ba[i] = _mm_cvtsi64_si32( vd );
vd = _mm_srli_si64( vd, 32 );
br[i] = _mm_cvtsi64_si32( vd );

vd = _mm_set_pi16 (bb[i],tb[i],bg[i],tg[i]);
vd = _mm_srli_pi16( vd, 8 );
vd = _mm_madd_pi16( va, vd );
bg[i] = _mm_cvtsi64_si32( vd );
vd = _mm_srli_si64( vd, 32 );
bb[i] = _mm_cvtsi64_si32( vd );
_mm_empty();

код выше делает:

ba[i] = ((ba[i]*ca)>>8)+ta[i] ;
br[i] = (br[i]*ca+tr[i]*a)>>8 ;
bg[i] = (bg[i]*ca+tg[i]*a)>>8 ;
bb[i] = (bb[i]*ca+tb[i]*a)>>8 ;

А вы попробуйте перейти на с _mm_..._pi.. на _mm_..._epi.. и во вторых никогда не транслировали этот код а ассемблер? Если нет и думаете что каждая команда _mm равна соответствующей инструкции то вы СИЛЬНО неправы!

GCC способен генерировать такое:
movapd %xmm0, -16(%rbp)
movapd -16(%rbp), %xmm0
movapd %xmm0, -32(%rbp)
movapd -32(%rbp), %xmm1
movapd %xmm1, -16(%rbp)
movapd -16(%rbp), %xmm0
movapd -64(%rbp), %xmm1
movapd %xmm1, -128(%rbp)
movapd -128(%rbp), %xmm1
movapd %xmm1, -134(%rbp)
movapd -134(%rbp), %xmm1
addpd %xmm1, %xmm0
movapd %xmm0, -16(%rbp)
movapd -16(%rbp), %xmm0
movapd %xmm0, -32(%rbp)
movapd %xmm1, -128(%rbp)
movapd -128(%rbp), %xmm1
movapd %xmm1, -134(%rbp)
Прикольно да? Правда при оптимизации -О2 большинство (но не вся!) подобная хня пропадает, но чем больше интристиков тем прикольней код... Так что увы вашей оптимизацией вы могли только испохабить код...