Как стать автором
Обновить

Комментарии 22

Хм, вообще относительная ошибка порядка 1E-7 (как в вашем случае) при работе с float'ами — это нормально. Она такая примерно всегда получается. Просто если у вас вещественные вычисления, то воспроизводимость результатов надо мерить не тупым равенством, а приближённым.
Это НЕ нормально, урезание количества знаков после запятой работать должно по одному и тому же алгоритму и выдавать всегда один и тот же (floor, round, ceil) результат, тем более, что они тестировали с тем же количеством нитей, выполняющих задачу. Фортран я не знаю, а вот код на C в подобных случаях всегда имело бы смысл сначало транслировать в ассемблерный код (gcc -S), его анализировать, а потом уже в бинарник, ну или для чистоты вообще реверсить с objdump, чтобы совсем точно.
Алгоритм «урезания» заложен в стандарте IEEE 754, и компилятор его поддерживает. Но если вы посмотрите пост, на который я ссылался в самом начале, то увидите пример, когда всё в соответствии с стандартом, а результат зависит от порядка суммирования, потому что критичны случаи, когда складывается большое и очень маленькое число.
Ну в Вашем случае оптимизация в компиляторе наотсебятила. Я вот тут еще нашел про выравнивание и –fp-model (or /fp) flag в «Using the Intel® Math Kernel Library (Intel® MKL) and Intel® Compilers to Obtain Run-to-Run Numerical Reproducible Results by Todd Rosenquist, Technical Consulting Engineer, Intel® Math Kernal Library and Shane Story, Manager of Intel® MKL Technical Strategy, Intel Corporation»

The Intel compiler also
provides threading via the OpenMP* model. With the latest compiler,
the reduction stage in threaded parallel sections can be forced to
provide reproducible results by setting
KMP_DETERMINISTIC_REDUCTION=yes.
Let’s consider the case where you would like the resulting executable
to produce the same results on the following three systems:
> One 4-core processor supporting SSE 4.2 instructions
> One 4-core processor supporting AVX instructions
> Two 4-core processors supporting AVX instructions
First, to ensure that Intel MKL functions provide reproducible results,
we must make sure all arrays are aligned on 16-byte boundaries (the
mkl_alloc() function suits this purpose). Then, ensure that the number
of threads used on each system remains fixed and does not vary
during the run. By default, Intel MKL uses as many threads as there
are cores, so you should fix the number of threads at four using the
mkl_set_num_threads() function. Finally, because, Intel MKL dispatches
an optimized code path based on the instruction set available, you
need to use the new reproducibility control to configure this: mkl_
cbwr_set(MKL_CBWR_SSE4_2. You will then also need to use an
appropriate –fp-model (or /fp) flag to ensure the compiler returns
consistent results. Applications threaded using OpenMP should also
specify a fixed number of threads (omp_set_num_threads(4)) and set
KMP_DETERMINISTIC_REDUCTION=yes.
П.С. Извиняюсь, в посте по ссылке все это есть. Прочитал пост после того, как предыдущий комментарий отправил.
НЛО прилетело и опубликовало эту надпись здесь
У меня сейчас gcc на Celeron (Bay Trail) флоаты для деления и умножения в XMM-регистры кладет и с SSE считает, а не в FPU.
Это правильно, сейчас FPU почти не используется.
Почему? Неужели так быстрее? Не вижу логики.
Классический FPU это стековая архитектура (привет наследию x87), SSE — регистровая
С регистрами часто работать гораздо быстрее чем со стеком
В FPU регистровый же стек, а не RAM.
Так загружать-выгружать данные то из этого стека куда-то надо. А поскольку старый набор инструкций FPU был создан когда сопроцессор был отдельным чипом, то FST/FIST выгружает данные из стека сугубо в память. И загрузкой та же история — данные можно брать только из памяти. Ну и в целом, как ни крути, неудобно манипулировать со стеком. Даже с учетом FXCH все равно нужно туда-сюда данные крутить, поскольку почти все инструкции работают сугубо с вершиной стека. Кроме того в суперскалярах то что практически все инструкции используют один и тот же регистр порядком убивает производительность, поскольку две операции с FPU одновременно выполнить уже становится проблемой.
Всё просто — так быстрее!
Понятно, что быстрее. Непонятно, почему векторная подсистема при использовании не по назначению быстрее, чем FPU при использовании по назначению. Там же такие сложности с загрузкой и чтением XMM-регистров — латентность, требования по выравниванию памяти и т. д.
Нет там сложностей с загрузкой и чтением XMM регистров. Многие SSE команды имеют векторный и скалярный варианты и скалярный имеет меньше ограничений чем векторный

Применительно к чтению XMM, чтение одного значения в SSE регистр — это MOVSS/MOVSD и там нет требований к выравниванию. Чтение вектора требующее выравнивания — это другая инструкция MOVAPS/MOVAPD и у нее есть и невыровненный эквивалент MOVUPS/MOVUPD, просто он более медленный
Не, ну тут действительно не будет один и тот же резульатат, много факторов, которые вне контроля программиста: состояния в L-кешах, какие из тредов ядра/процессоры немного раньше или немного позже посчитают. В данном случае важно знать допустимую погрешность при вычислениях, чтобы ее учитывать при разработке алгоритма.
Я не силен в фортране но при помощи подручной браузерной консоли удалось получить ситуацию как у вас (и без параллельного программирования) — не выполнения ассоциативности сложения.

var a = 1000; var b = 0.0004; console.log(((a + b) + b) == (a + (b + b)));
var a = 10000000; var b = 0.0004; console.log(((a + b) + b) == (a + (b + b)));

image

В свое время когда писал магистерскую на CUDA, я тоже наткнулся на эту проблему. Эта была программа по вычислению гауссова размытия изображения и тоже заметил что свертка матрицы на GPU и CPU немного отличались. На тот момент я решил эту проблему. Нужно было чтобы порядок операций над числами на GPU и CPU совпадал. Условно: ((a + b) + (a + b)) если я выполню это на CPU или если я выполню на GPU отдельно сложения в скобках а потом просуммирую результат на CPU то результат будет тот же.

Поправьте, если я не прав
Всё правильно — ассоциативность не выполняется. Но вы должны учитывать, что кроме написанного вами кода, есть ещё компилятор, который может его изменить, и библиотеки, которые могут не учитывать вашу необходимость в «правильном» порядке суммирования. Про это и разговор.
Вообще как мне кажется основная тема вашей публикации — это как раз не выполнение этой ассоциативности. Вы вроде указали даже команды которые отключают оптимизацию кода компилятора но ни слова не сказали о главной сути проблемы. Решение проблемы не понимая ее причины обычно не приводит к успеху. В обычном кольце всегда выполняет равенство 2 + 2 = 4, а в программировании мы можем сделать и так как у вас на рисунке. Чисто алгоритмически процессоры на GPU и CPU считают одинаково. При решении своей задачи, у меня сверки матриц были идентичными во всех случаях, но стоило заменить последовательность вычисления как результат получался уже немного другой.
В самом начале я давал ссылку на статью, где проблема описана намного шире и детальнее. Здесь речь в основном о частном случае, а именно использовании OpenMP и возможностей по контролю воспроизводимости там.
Зарегистрируйтесь на Хабре , чтобы оставить комментарий