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

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

Прирост от «4»-ядерного кода на С++ также невелик — 67% между 1 и 4 ядрами (оба float). Списывать можно на медленный контролер памяти и 32-битную Винду.

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


Кстати, у вас в коде на C++ есть несколько new, но ни одного delete...

Кстати, у вас в коде на C++ есть несколько new, но ни одного delete...
Как ни странно, но в данном случае освобождение памяти средствами приложения вместо средств ОС только немного замедлит программу. ;)

На самом деле списывать надо на то, что вы не специались по распараллеливанию, как и не специалист по оптимизации...
О да, Intel VTune Amplifier спасет отца русской демократии:
https://habrahabr.ru/company/intel/blog/310842/
Мне было бы интересно, поможет ли коду Интеловский компилятор, но у меня его нет :)
сам компилятор поможет несильно, а такой хороший профайлер покажет много полезного.
Да, код делался только под i7. Скорее всего, для Атлона нужно смотреть отдельно.
Cтоит попробовать вынести константные подвыражения из циклов.

Судя по ограничению скоростью памяти, вам стоит попробовать учесть иерархию памяти, чтобы повысить локальность обращений к памяти (процент попадания в кэш).
https://habrahabr.ru/post/312078/

Далее -march=native, векторизация (-sse2, avx), интринсики, профайлеры…

После чего становится грустно и взоры обращаются обратно к ANSYS ;)
Не вижу возможности вынести из циклов. Все зависит от [i][j]
-march=native проверю.
Извините, я компилирую на одной, а проверяю на разных.
Использовался вроде ‘nehalem’, проверял что от ‘sandybridge’ пользы нет (проверю, по моему у меня опции не так называются)
В интринсиках разбираться не планирую
Не вижу возможности вынести из циклов. Все зависит от [i][j]

     for (int t = 1; t <= tMax; t++) {                    // begin main loop
            float tt = Math.min(t * s + 10, ny - 1);
            //gauss
            switch (method) {
                case "cos":        

да, да конечно. Я понимаю, что это не самый внутренний цикл, но сравнение строк… Будем молиться, что компилятор сам раздвоит циклы для разных case (о чем бы было неплохо ему напомнить явно). Уж лучше лямбды/указатели на ф-ии, виртуальные методы в явном виде. Накрайняк — глобальный копипаст (как обычно, это всего быстрее:) ).

в
foreach (i)
  x=f(i)
  exp(-pow(x, 2) / w / w - (t - 1) * dt / tau)
я правильно понял, что
  • вся правя часть изменяется-таки от i и x и не может быть вынесена из цикла? Что-то я не верю, что компилятор матан учил лучше меня. Или там приколы с потерей точности?
  • pow(x, 2) это не x*x и вызов ф-ции (float,float) имеет смысл?

я Понимаю, что это вроде как мелкие придирки, но есть и более интересные замечания:
float  **Hy = new float *[nx];
— массив указателей на массивы. Итого 2 обращения к памяти. Лучше уж глобальный
float XX[X][Y]
, что сведет операцию доступа к XX[x*Y+y]. Замечание: правильно назначьте X и Y, чтобы доступ к памяти был последовательным и хорошо ложился на prefetcher. И не надо бояться, что бинарник распухнет: эти массивы будут кратко описаны в сегменте zero-data.

Потом смотрим, что весь код у вас есть несколько циклов по nx при все увеличивающимся ny. Так мы имеем, что к следующему циклу по x строка уже выпала из кэша. Возможно имеет смысл сделать
for (int t = 1; t <= tMax; t++) {
  for (i = 1; i <= nx - 1; i++) {
    //gauss //для единственного i
    // boundary conditions //для единственного i
    ....

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

Дальше. Как уже говорили ниже
Все зависит от [i][j]
, но не для всех (i,j), а некой окрестности точки (i,j), что позволяет как повысить попадания в кэш, так и подумать за уменьшение выделенной памяти и массовое распаралеливание.
Спасибо за анализ.
— Глобальный копипаст делать точно не хочу — слишком много нужно. Сделаю вызов функции или ссылку на функцию.
— exp(-pow(x, 2) — это какое-то локальное помутнение. Все квадраты расписал, а этот нет.
— exp( — (t — 1) * dt / tau), пожалуй, посчитаю заранее
— в С не шарю, именно об этом спрашивал. Но для Джавы 1-мерный массив давал проигрыш в производительности.
for (i = 1; i <= nx — 1; i++) {
//gauss //для единственного i
// boundary conditions //для единственного i

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

— порядок циклов вроде правильно ложится в кеш.
— я попробовал дополнительно делить массив еще на 4 части, имеем небольшую потерю производительности.
Еще моменты по оптимизации под 4 ядра. Я думал поделить массив пополам изначально, чтоб два потока обращались к разным учаскам памяти.

Если Вы упираетесь в скорость памяти — попробуйте разбить задачу на кусочки по-другому. Например, можно взять область (100+2x, 100+2x), и на её основе посчитать вперёд сразу на х шагов центральный кусочек размером 100*100. В идеале должно получаться, чтобы необходимая для "кусочка" память была не слишком большой и помещалась в кеш процессора — тогда у Вас на первом шаге будет чтение из ram и на следующих (х-1) шагах будет использоваться кеш.

Понял вашу мысль. Думаю, ничего хорошего не выйдет. Слишком большие накладные расходы на 2х.
это невозможно. Нельзя выделить какой-то кусок и считать его больше одного шага. Каждая точка на следующем шаге зависит от четырех соседних.
Можно только попытаться разбить циклы внутри одного шага.

нет, если взять х небольшим, то должно работать быстрее.
Например, для x=5 Вы рассчитаете квадраты размерами 110, 108,… 102. Конечно,
это немного больше вычислений по сравнению с вычислениями и чтением данных из 5 квадратов 102*102, (процентов на ~20%), но нагрузка на оперативную память упадёт почти в 5 раз — а именно она и является бутылочным горлышком.

Безотносительно самой задачи — код на С++ просто безумен. Уже много раз повторяли, в том числе и здесь, на habr-e, что нельзя сравнивать Java и С++ просто переведя алгоритм на уровне синтаксиса. Ни один вменяемый программист на С++ такую программу не напишет. Оператор new в Java и С++ не являются эквивалентами и new в с++ используется совсем по другому. Обращение с двумерным массивом поэлементно через два индекса первый из которых прыгает по таблице указателей это вообще за гранью.
На Хабре пишут, что многие научные расчеты пишут на фортране.
Даст ли выигрыш переписывание и оптимизация на него?
Кроме того, что в Фортране я совершенно не разбираюсь.
На Фортране пишут потому, что под него есть тонны математических/научных библиотек
https://benchmarksgame.alioth.debian.org/u64q/fortran.html
Фортран или на уровне чистого С, или медленнее.
Уточняю. Нашел книжку 2006года по распараллеливанью.
Там ФДТД используют на Fortran 90 с MPI
Не даст. Фортран просто имеет более высокоуровневые оптимизированные абстракции для работы с математикой. Т.е. на нем есть уже готовые примитивы, например для работы с матрицами, из «коробки». Основная причина использования фортрана для расчетов, это то, что он старше С и на нем написано куча библиотек для научных расчетов, а математики и физики очень консервативны, и, в основном, не любят программирование.
3) Проверить, могут ли лябмды/стримы ускорить основной цикл или дополнителные.
Не помогут. Лямбды и стримы — чтобы сделать код более удобным. В лучшем случае вы получите приблизительно такой же native-код с той же производительностью. В худшем случае — всё станет заметно медленнее.
Возможно некоторые циклы, в которых есть обращение к нескольким строкам, будет эффективнее разбить на несколько проходов, где идут обращения только к одной строке.

Например вместо одного прохода со строчкой
Ez[i][j] += e[i][j] * ((Hx[i][j — 1] — Hx[i][j] + Hy[i][j] — Hy[i — 1][j]));
сделать два прохода:
1:
Ez[i][j] += e[i][j] * ((Hx[i][j — 1] — Hx[i][j] + Hy[i][j]));
2:
Ez[i][j] -= e[i][j] * Hy[i — 1][j];
… а возможно даже и 3, чтобы одновременно участвовало меньше матриц:
1.
Ez[i][j] += e[i][j] * ((Hx[i][j — 1] — Hx[i][j]));
2.
Ez[i][j] += e[i][j] * Hy[i][j];
3.
Ez[i][j] -= e[i][j] * Hy[i — 1][j];
Проверил такую идею. Производительность упала в 1,6 раз
Отличная оценка, спасибо за неудачную идею.
В таком способе нужно читать эпсилон 3 раза.
Вместо 5 чтений 7+лишение записи. Плюс лишняя математика.

А в Н такого и не сделаешь
Раз производительность упала в 1.6 раз, а не в 3, значит вы именно упираетесь по памяти и вам стоит изучить варианты реструктурирования ваших циклов. В идеале — добиться векторизации всего, что только можно.
Основная идея для оптимизации в вашем коде, если не брать векторизацию и развертывание, это вынести из внутренних циклов все константные данные. Например в таком коде:
void thH (float** &Ez, float** &Hx, float** &Hy,
          int iBegin, int iEnd, float tt)
{
    for (int i = iBegin; i <= iEnd; i++)
    {
        for (int j = 1; j <= tt; j++)
        {
            Hx[i][j] += Ez[i][j] - Ez[i][j + 1];
            Hy[i][j] += Ez[i + 1][j] - Ez[i][j];
        }
    }
}

можно сделать так:
void thH (float**& Ez, float**& Hx, float**& Hy, int iBegin, int iEnd, float tt)
{
    for (int i = iBegin; i <= iEnd; i++)
    {
        float* Ez_row = Ez[i];
        float* Ez_row_plus1 = Ez[i + 1];
        float* Hx_row = Hx[i];
        float* Hy_row = Hy[i];

        for (int j = 1; j <= tt; j++)
        {
            Hx_row[j] += Ez_row[j] - Ez_row[j + 1];
            Hy_row[j] += Ez_row_plus1[j] - Ez_row[j];
        }
    }
}
После компиляции бинарник мало отличается. Разница производительности в пределах погрешности.
Замечательно, значит одно из двух:
1. Hotpoint не там и основное время тратится не в «основном» цикле
2. Компилятор сам сделал эту оптимизацию

А вы пробовали избавиться от new в циклах? И не совсем понял, вы какой профилировщик используете в С++?
1. Два основных цикла тратят до 84-85% времени (на Джаве также)
2. Да, компиляторы нынче умные
Профилировщик в С++ никакой, хотя вроде в VisualStudio смотрел.
Посмотрел в профайлере VS действительно, тормозят основные циклы, причем в них делаются простые вычисления. Развертка цикла ничего не дает у меня на машине. Значит тормозит память и нужно попытаться изменить layout матриц.
Как идея, сделать все данные более локальными. Например:
1. Сделать одну матрицу вместо 4-х где хранить коэффициенты друг за другом или как структуру.
2. Сделать одномерный массив и руками интерировать по строкам.
2.У меня и у https://habrahabr.ru/post/319216/#comment_10005284
1-мерный массив уменьшает производительность.
1. А вот хранить все локально — интересный вариант. Только все циклы придется ужасно переписывать.
Попробовал одномерный массив — подтверждаю. Ничего не дает.
Ради прикола, попробовал заменить в одномерной версии float* на char* (понимаю что считать будет не правильно) — выигрыш всего 25%. Плюс, 700 * sizeof(char) * 1024 = 700кб — это смешно по памяти (не похоже на память). Значит вам нужна векторизация, он просто не успевает считать.
К сожалению это специфично для компилятора и сами матрицы должны быть выровнены. Быстро не проверить.
Сделал еще экстремальнее — заменил 4 массива на int8_t (или это плохой тип?), с сохранением корректности всей математики. Выигрыш при 4096х700 1,26 раза. Чем больше область, тем выше.
Почему же был высокий прирост от 64бита к 32?
1024*700 — это тестовая задача, какую смог дождаться в Матлабе. С++ детально не изучал, в Джаве наибольший упор в скорость памяти при 4096х500
Еще мысли:
Все выглядит так, как-будто циклы просто очень медленно считаются, хотя тела маленькие и операции элементарны. Здесь хорошо бы помог VTune Analizer или что-то подобное. У меня его к сожалению нет, а профайлер от Microsoft очень грубый, он показывает сколько CPU затрачивает строчка кода, но при этом нельзя посмотреть ассемблерный код. VTune в таких случаях показывает более точную информацию. Может быть у вас где-то зависимость по данным и вся параллелизация убивается на корню.
Провел собственный анализ.
Много тратит запись элемента массива. Если принять время выполнения 3 циклов за 100у.е., то
каждая запись(+=) из 3х тратит 20,
запись на даблах 31 (просто для сравнения)
однократное чтение 6-7
пятикратное чтение 27
(оценки приблизительны, внимательный читатель поймет, что мы уже ниже 0)
Ох, а как вы это померили? Дело в том, что современные процессоры очень сложно устроены и при перетасовке операций (при аналогичном результате) результаты могут быть совсем другими. Идея ж простая (теоретически): есть X — исполнительных блоков, способных работать одновременно), плюс Y — стадий конвейера. Если операции без зависимостей друг от друга, то пропускная способность будет XY, если каждая операция зависит от другой, то способность будет 1 (все условно). Только вот по С — коду понять, где затык не всегда можно. А мерить время сами операции бесполезно, важно только суммарное время.
Еще вычитал:
«Согласно протоколам согласованности кеша (cache coherence protocols), если ядро пишет в кеш-строку, то строка в другом ядре, ссылающаяся на ту же память, признаётся недействительной (пробуксовка кеша, cache trashing). В результате при каждой операции записи возникают блокировки памяти
У меня такое возможно на середине массива. То есть задача влазит в кеш, но пишется с проблемами.
Все верно, но у вас такого не будет, так как, пока первое ядро дойдет до середины второе уже вытеснит все строки кеша и будет обрабатывать конец второй половины. Ваш алгоритм должен работать в 50 раз быстрее, при правильной оптимизации. Вот у меня, на той же машине есть цикл в 10 раз сложнее, а молотит в 10 раз быстрее чем ваш.
P.S. У вас нет случайно ссылки на ваш проект под студию, а-то я нечаянно похерил свой?
Я ленивый, делал только ГСС.
Под студию был 1-ядерный.
Постараюсь сегодня сделать 2-4.
Одноядерный тоже пойдет для экспериментов. Обычно если потоки абсолютно не пересекаются, то их общая производительность просто суммируется.
Мой старый проект завален какими-то левыми файлами, создал чистый 4-ядерный.
https://cloud.mail.ru/public/8K8G/VCRZSJQTt
Есть профит по сравнению с ГСС во всех случаях.
Скоро обновлю в статье
Спасибо! Попробую с главными циклами поэксперементировать. Жалко что у вас код не С++. Было бы проще и быстрее его менять, а так сплошное дублирование. Меня не покидает ощущение что я погрузился в средние века. Уж очень к С++ быстро привыкаешь.
Не обязательно выносить общие вычисления за пределы циклов — компилятор с этим замечательно справляется. Тоже самое касается констант
Смотря какой. Вот я сейчас смотрю профайлером. У меня VS 2013 c компилятором по умолчанию. Она ничего не выносит в коде автора. У меня вообще код падает, так как в главном цикле идет сравнение с tt которое float и в зависимости от причуд плавающей арифметики периодически происходит выход за границы массива.
vlanko — в С++ код при смешивании типов приводится к более общему типу, т.е. j приводится к float
да, float — это мой глюк. Пока циклы были <, а не <= все было нормально. Заменю все на инт.
Какой у вас суровый брат. А что, если поставить еще одну ОС рядышком со своим набором драйверов (CUDA)? Или там какие то другие причины?
Касаемо кода с++, по беглому просмотру, он конечно нуждается в оптимизации и, главное, в векторизации. Можно воспользоваться готовыми библиотеками с векторными классами, например eigen. Если будет время, обязательно посмотрю во что там компилируется этот код, выносятся ли условия за пределы цикла, векторизует ли код компилятор.
Кстати, а Java умеет векторизовать?
Ускорил код более чем в два раза просто развернул «горячие» циклы:
что было
void thE (float** &Ez, float** &Hx, float** &Hy, float** &e,
		  int iBegin, int iEnd, float tt)
{
	for (int i = iBegin; i <= iEnd; i++)          // main Ez
	{
		for (int j = 2; j <= tt; j++)
		{
			Ez[i][j] += e[i][j] * ((Hx[i][j - 1] - Hx[i][j] + Hy[i][j] - Hy[i - 1][j]));
		}
	}
}
void thH (float** &Ez, float** &Hx, float** &Hy,
		  int iBegin, int iEnd, float tt)
{
	for (int i = iBegin; i <= iEnd; i++)
	{
		for (int j = 1; j <= tt; j++)
		{
			Hx[i][j] += Ez[i][j] - Ez[i][j + 1];
			Hy[i][j] += Ez[i + 1][j] - Ez[i][j];
		}
	}
}

что стало
void thE (float** &Ez, float** &Hx, float** &Hy, float** &e,
		  int iBegin, int iEnd, float tt)
{
	asm("# point 2 start");
	for (int i = iBegin; i <= iEnd; i++)          // main Ez
	{
		auto ttt = tt / 4;
		for (int j = 2; j <= ttt; j+=4)
		{
			int j1 = j+0;
			int j2 = j+1;
			int j3 = j+2;
			int j4 = j+3;
			Ez[i][j1] += e[i][j1] * ((Hx[i][j1 - 1] - Hx[i][j1] + Hy[i][j1] - Hy[i - 1][j1]));
			Ez[i][j2] += e[i][j2] * ((Hx[i][j2 - 1] - Hx[i][j2] + Hy[i][j2] - Hy[i - 1][j2]));
			Ez[i][j3] += e[i][j3] * ((Hx[i][j3 - 1] - Hx[i][j3] + Hy[i][j3] - Hy[i - 1][j3]));
			Ez[i][j4] += e[i][j4] * ((Hx[i][j4 - 1] - Hx[i][j4] + Hy[i][j4] - Hy[i - 1][j4]));
		}
	}
	asm("# point 2 end");
}
void thH (float** &Ez, float** &Hx, float** &Hy,
		  int iBegin, int iEnd, float tt)
{
	asm("# point 3");
	for (int i = iBegin; i <= iEnd; i++)
	{
		auto ttt = tt / 4;
		for (int j = 1; j <= ttt; j+=4)
		{
			int j1 = j+0;
			int j2 = j+1;
			int j3 = j+2;
			int j4 = j+3;
			
			Hx[i][j1] += Ez[i][j1] - Ez[i][j1 + 1];
			Hy[i][j1] += Ez[i + 1][j1] - Ez[i][j1];
			
			Hx[i][j2] += Ez[i][j2] - Ez[i][j2 + 1];
			Hy[i][j2] += Ez[i + 1][j2] - Ez[i][j2];
			
			Hx[i][j3] += Ez[i][j3] - Ez[i][j3 + 1];
			Hy[i][j3] += Ez[i + 1][j3] - Ez[i][j3];
			
			Hx[i][j4] += Ez[i][j4] - Ez[i][j4 + 1];
			Hy[i][j4] += Ez[i + 1][j4] - Ez[i][j4];
		}
	}
}

Разумеется, граничные условия нужно дообрабатывать отдельно или указывать границы кратные количеству развертываний цикла. Но это мелочи и на производительность практически никак не скажется.
Компилирую так:
g++.exe main.cpp -std=c++11 -O3
Асм соответственно смотрю смотрю так:
g++.exe main.cpp -std=c++11 -O3 -S
Компилятор этот код не векторизовал, сейчас попробую ему в этом помочь.
Моя система: i5-3317U 1.7gHz, два ядра (+HT), ОЗУ 4gb, 1600mGz
Абсолютно ничего, это просто метка/комментарий для ассемблерного кода. Если скомпилить с флагом -S, то рядом с ехешником создастся файл с расширением .s — в нем будет ассемблерный код, то, во что преобразовывает компилятор исходники. Метка # point 2 start нужна для того, чтобы быстренько найти этот участок в асм коде. Вот, например:
Какой-то кусок кода на с++
int j3 = j+2;
int j4 = j+3;
			
asm("# point 2 sum1 start");
sum1[0] = Hx[i][j1 - 1] - Hx[i][j1];
sum1[1] = Hx[i][j2 - 1] - Hx[i][j2];
sum1[2] = Hx[i][j3 - 1] - Hx[i][j3];
sum1[3] = Hx[i][j4 - 1] - Hx[i][j4];
asm("# point 2 sum1 end");
			
sum2[0] = Hy[i][j1] - Hy[i - 1][j1];
sum2[1] = Hy[i][j2] - Hy[i - 1][j2];

Этот же код полученный через флаг -S
movl	(%ecx,%ebp), %ecx
	movl	(%edx,%ebp), %edx
	.p2align 4,,10
L7:
/APP
 # 28 "main.cpp" 1
	# point 2 sum1 start
 # 0 "" 2
/NO_APP
	vmovss	8(%ecx), %xmm0
	vmovss	12(%ecx), %xmm3
	vmovss	16(%ecx), %xmm2
/APP
 # 33 "main.cpp" 1
	# point 2 sum1 end
 # 0 "" 2
/NO_APP
	vaddss	4(%edi,%eax,4), %xmm0, %xmm4
	addl	$16, %ecx

То есть зная метку через ctrl+F легко находится нужный участок в асм коде
Ахтунг!
Применил ваш код, и вышла фигня.
Потом понял:
ttt = tt / 4;
Вы просто уменьшили объем вычислений в 4 раза
Я не уменьшал количество вычислений, я уменьшил количество итераций цикла в четыре раза, а в самом теле цикла увеличил в четыре раза. А фигня у вас получилась скорее всего изза того, что вы не обрабатываете данные на границах — я об этом писал.
Вы уменьшили в 4 раза.
Я же посмотрел визуализацию на Яве. Фигня была именно в том, что обрабатывается только первая четверть картинки.
В вашем коде:
auto ttt = tt / 4;
		for (int j = 2; j <= ttt; j+=4)
		{
			int j1 = j+0;
			int j2 = j+1;
			int j3 = j+2;
			int j4 = j+3;

Нужно:
а) int j1 = j * 4+0;
или б) j <= tt- 4
или в) j <= ttt * 4. (Но эти два варианта глупые)
В вашем коде максимальное значение переменных tt / 4+3 << tt
j <= ttt
ttt = tt / 4;
да, действительно, мой косяк. В таком случае да, раз вы упираетесь в память, что либо с вычислениями делать бесполезно
сравните два варианта кода:
for (int i = 0; i = 10; i++)
    a[i] = b[i] + c[i];

for (int i = 0; i = 10/4; i+=4)
{
    a[i] = b[i] + c[i];
    a[i+1] = b[i+1] + c[i+1];
    a[i+2] = b[i+2] + c[i+2];
    a[i+3] = b[i+3] + c[i+3];
}

Во втором случае я развернул цикл, но так как 10 итераций не делятся на цело на 4, то получается оставшиеся две итерации нужно обработать отдельно:
for (int i = 0; i = 10/4; i+=4)
{
    a[i] = b[i] + c[i];
    a[i+1] = b[i+1] + c[i+1];
    a[i+2] = b[i+2] + c[i+2];
    a[i+3] = b[i+3] + c[i+3];
}
a[8] = b[8] + c[8];
a[9] = b[9] + c[9];

Этого я не сделал в вашем алгоритме, о чем сразу же и написал:
Разумеется, граничные условия нужно дообрабатывать отдельно или указывать границы кратные количеству развертываний цикла. Но это мелочи и на производительность практически никак не скажется.

Советую вам корректно развернуть циклы, так как потенциально это может дать сильный прирост производительности на больших итерациях
В моем случае я вообще не вижу смысла в таком разворачивании циклов. У меня запас производительности процессора примерно в 48 раз выше, чем скорость памяти.
Пока на Джаве разворачивание только Е дало падение производительности на 3 % (в пределах погрешности)
vlanko, как выяснилось векторизовать данные циклы невозможно по причине смещения данных (Hx[i][j — 1] и Hx[i][j], например).
Но развертывание цикла, как я показал выше, вполне имеет место быть, у меня время просчета сократилось с 4500мс до 1900мс
Может быть возможно изменить алгоритм так, чтобы данные не пересекались — тогда можно будет векторизовать.
По идее любая разностная схема считает производные именно как разницу между соседними элементами и обойти это никак не выйдет.
Кстати, можно воспользоваться встроенными векторными классами
С ними так же просто, как в шейдерном коде. Если вам не удастся векторизовать ваш алгоритм, то и на CUDA особого прироста не ждите.
Проверить скорость чистого С (не ++). По benchmarksgame он всегда быстрее.

У вас код и так написан на почти чистом Си. Так что никакого прироста не получится.

Извините, но почему не meep? Или любой другой коммерческий вариант (rsoft, lumeric)? Я понимаю что лицензирование коммерческих продуктов стоит денег, но было бы интересно сравнить ваш код с meepом.

Я понимаю что писать свой FDTD это занимательно (лично знаю человека который из принципа написал свой 2D FDTD на CUDA), но для исследовательских целей почему бы не использовать готовый продукт и не строить свой велосипед?

И еще по поводу Матлаба. Циклы и условные ветвления (а так же вызовы анонимных функций) в нем действительно медленные, но вот встроенные алгоритмы работы с массивами (сложение, умножение матриц, нахождение определителей, собственных значений и т.д.) давольно быстры. На мой первый взгляд в матлабе возможно реализовать FDTD используя только лишь один внешний цикл для времени, все остальное можно сделать используя матричные операции.
Когда я начинал разбираться в вопросе в 2010м, Meep был непонятно чем. Сейчас, наверное, использовать его оптимальнее, чем самому писать на С++.
Да, я писал. Матлаб очень быстро умножает матрицы и с нормальной скоростью берет функции от массивов. Возможно, получится такое сделать. Но это потребует создание 4 новых матриц для вычитания. Боюсь, со скоростью памяти будет хуже.
Я переписал в Матлабе 2-мерный цикл полностью на вычитание матриц:
код
Ez0(2:end-1, 2:tt) = Ez(2:end-1, 2:tt)+e(2:end-1, 2:tt).*(Hx(2:end, 1:tt-1)-Hx(2:end, 2:tt)+Hy(2:end, 2:tt)-Hy(1:end-1, 2:tt));
Hx0(1:end, 1:tt) = Hx(1:end, 1:tt) + dH *(Ez(1:end-1, 1:tt) - Ez(1:end-1, 2:tt+1));
Hy0(1:end, 1:tt) = Hy(1:end, 1:tt) + dH *(Ez(2:end, 1:tt) - Ez(1:end-1, 1:tt));


Выигрыш от 3,1х для малых до 3,8х до больших матриц. Если оптимизировать идеально остальное, может выжмется 4,5х.
В любом случае, спасибо за идею.
Еще можно посмотреть в матлабовском профайлере что происходит (profile viewer). Я не пытаюсь сказать что маталб может провести все вычисления так же быстро как и си/ява, плюс возможно в матлабе не получится так просто сделать распараллеливание (если честно кроме parfor я не знаю других способов распаралелливания вычислений, хотя можно попробовать запускать ява код из матлаба). Но возможности матлаба по визуализации подкупают.

Если вы в лс опишите структуру и скажете что посчитать, я могу сделать .ctl файл для meep и прогнать вычисления.
Предлагаю начать с того, что параллелизация на с++ выполнена с ошибкой.

Видимо автором ожидается, что цикл по t ожидает завершения стартованных по join() нитей, но это не так.
Теоретически да. Но я проверяю, и у меня полный бардак в порядке исполнения при использовании mingw-gcc
Уточняю. В исходном коде (и в Яве и в С++), заявленном как 4-х поточный, выполняется следующее
— запуск 2х потоков
— ожидание их завершения
— запуск еще 2х потоков
— ожидание их завершения

Какое может быть ускорение? Это обычный 2-поточный код
Нет, не двухпоточный.
— запуск 2х потоков
— в каждом из двух потоков еще запуск по 2 потока
То есть два треда, в каждом запускаются сотни раз еще по паре тредов, в итоге одновременно запущено в среднем 4 потока.
Принято, похоже сплю.

В mingw действительно что то не так в потоках, некоторые версии mingw с потоками не работают совсем, некоторые абы как. В общем, принудительно поставив пару mutex'ов, заставил отработать как надо.
На результат не повлияло.

Простая замена на линейный массив [y + maxx *x] привела к большому падению скорости.
Дальше надо лезть в ассемблер.

И главный вопрос — почему при распараллеливании не растет пропорционально скорость?

FFT хорошая задача, еще в каком то старом тесте была
Потоки — стандартные из C++11.
Именно в mingw нормально работают.
Да, одномерный массив у меня тоже давал падение скорости.
Имхо, ускорение слабое, потому что даже 2-поточный код близок к скорости памяти.
Да, все именно так и планировалось. Вопрос у меня в том, можно ли на С++ эффективнее порождать потоки (с меньшими затратами)
Можно использовать пул потоков
Поскольку в быстром Фурье не разбираюсь, взял первый попавшийся. Минус — ширина только степень двойки.


Это стандартная особенность БПФ.
Вообще, последние i7 имеют от десяти до двадцати мег кэша L3 — хотя он тормознутее, конечно — но ограничение скорости по памяти должно быть слабее, чем для других.
А вообще, получается необходимо дробить алгоритм на отдельные «островки», хорошо попадающие в кэш.
А на БПФ я наезжаю потому, что Матлаб умеет как-угодно без существенных потерь скорости.
Да, но у моих процей только 6-8 МБ кеша.
Малые задачи у меня выполняются быстрее — эквивалент 30 ГБ/сек против 17 ГБ/сек, (память 20около).
То есть лезут в кеш.
Но тупое деление большой области на 4 куска по і делает немного хуже.
Буду проверять развернутые циклы.
Если «как угодно» — то это либо классическое дополнение нулями области значений до требуемой размерности, либо какие другие способы дополнения.

В Матлабе много чего наворочено «внутри» — но, как всегда в таких случаях — вы берёте готовое.
И «под капот» вас не пускают.

Кстати, дурацкая идея — я гляжу, что во всех научных расчётах есть тенденция разбивать всё на составляющие — Ez, Ex, Ey и т.д.
А если поступить наоборот?
Ну т.е. у нас фиксируется пространство «точек», в каждой точке имеется набор свойств — запись либо массив, Ex, Ey, Ez, Hx, Hy, Hz — и, возможно — локальные производные.
Тогда локальные расчёты автоматически локализуются ;)) — для двумерного случая у нас будет три «слайса» из массива по три элемента в каждом. Или даже два.
Возможно алгоритм будет работать «компактнее».
В таком случае размер структуры становится некратен ^2 и плохеет префетчеру и кратности размеру строки кэша.
Intel VTune Amplifier в https://habrahabr.ru/company/intel/blog/310842/ предлагает обратно разбить на несколько массивов
Кхе.
Падами забить до требуемого размера?
(Двумерный случай, кстати — должен быть кратен, трёхмерный придётся добивать до четырёх измерений. %)
Думаю, имелось в виду, что массив из 3-6 компонент не кратен 4м.
Но я мог бы сложить в одну матрицу Ez, Hx, Hy, e — как раз 2^2
Кстати про преобразование фурье, в матлаб используется библиотека fftw. Сама библиотека написана на си. Там не совсем «как угодно», производительность возрастает если длина факторизуется на малые простые числа.
Ну, в принципе, известно, что БПФ предъявляет требования к размеру набора исходных данных.
Наиболе распространённый вариант — Кули-Тьюки — предполагает размер датасета в степень двойки, другие варианты прореживания — хотя-бы факторизуемость.
Кстати, свести этот метод к методу конечных элементов нельзя?
А то под МКЭ существуют свои пакеты — от коммерческих до бесплатной Salome.
В общем, я решил задачку.

Посмотрев ассемблерный код на отличном сайте https://godbolt.org, выяснилось, что gcc в 32-битной версии не использует SSE и подобное.

В итоге нужно брать нормальную 64-битную версию и компилировать примерно так
>-O2 -lstdc++ -march=core-avx2 -m64

Для матрицы 4906х1000 у меня результат изменился со 154с на 29.5с, т.е грубо в 5(пять) раз.
Исходник автора.
Уточнение, для 32-бит просто нужно дополнительно указывать -mfpmath=sse

Т.е так >-lstdc++ -O2 -march=core-avx2 -m32 -o sie32.exe -mfpmath=sse

Это тоже решает проблему
У меня 32-битный gcc 4.92
Включение опции -mfpmath=sse не изменяет ничего.
Полный набор ваших опций увеличивает время на 5%
Если что, мои опции mingw32-g++.exe -Wall -fexceptions -O2 -march=corei7-avx -O3 -std=c++11
Буду пытаться 64-битный компилятор. Под винду сложно найти нормальное.
http://netassist.dl.sourceforge.net/project/mingw-w64/Toolchains%20targetting%20Win32/Personal%20Builds/mingw-builds/installer/mingw-w64-install.exe

Это загрузчик. Нужно брать вариант с posix-threads. Оба варианта версии 6.2 i686 и x86-64 работают корректно.

-O3 не рекомендую использовать — глючная
Собственно, вы не могли бы выложить готовое приложение, уточнив размер массива (4096*1000 ?)
И добавив в конце FDTDmainFunction
код типа
double sum = 0;
    for (i = 0; i < nx+1; i++)
    {
        for (j=0; j<ny+1; j++)
        {
            sum+=abs(Ez[i][j]);
        }
    }
    cout << " " << sum << " ";


Чтоб я мог сравнить скорость программы и результат
Не вопрос https://yadi.sk/d/87Kyv5eQ38ZYHf

Там исходник, необходимый рантайм, батник для компиляции под AVX, и две версии — SIMD и обычная O2 — см батник
Требует libgcc_s_dw2-1.dll
После установки гсс оно появится?
Сорри, недоглядел. Перевыложил по той же ссылке

Конечно. Рантайм GCC лежит в папке bin компилятора. К сожалению с ним затруднительно слинковаться статически.
У вас что-то не так с компилятором. Возможно, О3 рулит.
i7-4771
Моя версия выполняется 22,3 секунд
ваша в моем компиляторе 22,4 (конечно, это тоже самое)
Ваш АВХ 23,2
Ваша без инструкций 93,5 секунд
Все так. Это всего лишь значит, что у вас по умолчанию используется sse математика.

Для этого процессора архитектура -march=haswell
В этом тесте -O3 вроде ничего не ломает.

У меня сейчас i7-4702mq, той же архитектуры.

Итого
gcc (i686-posix-dwarf-rev1, Built by MinGW-W64 project) 6.2.0
gcc.exe aut-res.cpp -lstdc++ -O3 -march=haswell -mfpmath=sse -o aut-avx.exe
4096х1000 результат 22.8с

Java(TM) SE Runtime Environment (build 1.8.0_112-b15)
Java HotSpot(TM) 64-Bit Server VM (build 25.112-b15, mixed mode)
4096х1000 результат Full calculate time 37.323 sec

Компилировалось javac.exe -encoding UTF-8 BasicEx.java
Потребление памяти 132Мб против 222Мб у Явы.
131 МБайт — это и есть объем массивов двух компонент. А Джава-машина сама жрет много.
Но главное — я проверил gcc 64 6.2 -march=native. На широких задачах (16384 — лучше для 4 ядер) на 35% быстрее старого.
Но VisualStudio еще на 3-10% быстрее.
for (j = 2; j <= tt; j++) {//можно ли тут оптимизировать порядок доступа к ячейкам?

Можно уменьшить количество обращений к памяти. Примерно так:
final int lastNx = nx - 1;
for (int i = 1; i <= lastNx; i++) {
    double hx0 = Hx[i][1];
    double hy0 = Hy[i][1];
    double hx, hy;
    for (int j = 2; j <= tt; j++) {
        hx = Hx[i][j];
        hy = Hy[i][j];
        Ez[i][j] += e[i][j] * (hx0 - hx + hy - hy0);
	
        hx0 = hx;
        hy0 = hy;
    }
}

Осталось только посчитать магнитное поле

Аналогично оптимизируется.

co.E[i][j] = co.E[i][j] * co.E[i][j] + si.E[i][j] * si.E[i][j];

Опять избыточные обращения к памяти. Лучше так:
double co = co.E[i][j];
double si = si.E[i][j];
co.E[i][j] = co * co + si * si;

Компилятор может такие вещи оптимизировать. Но всегда полагаться на оптимизатор я бы не советовал.

К тому же:
  • Объявления счетчиков циклов (i,j) нельзя выносить за пределы циклов.
  • Циклы с условиями вида «i < nx + 1» лучше разворачивать в обратную сторону либо вычислять nx+1 заранее.
  • В цикле условия на счетчик типа «j < begin» и «j — begin» следует избегать. Стоит разбить на 2 цикла: до «begin» и после.
  • Всё, что выносится за пределы цикла, надо таки вынести. Например: «dx / period», "(t — 1) * dt / tau", "(t — 1) * dt) * omega", "(double) nx / 2 + 0.5" и т.д.
  • Где посмотреть полный код на java? Желательно с юнит-тестами… Простор для оптимизации тут явно есть. Ну и параметры запуска надо бы озвучить.
Можно ли этот код использовать для расчёта акустических волн?
Акустикой я никогда не занимался. Знаю, что впервые такие сетки использовались в гидродинамике, для векторов давления и скорости.
Нужно смотреть конкретные уравнения. Возможно, в акустике будет проще.
По сути, это уравнения Максвелла — то, что до волн. Если нужны просто волны — там более быстрые методы.
А можно ли увидеть код на MATLAB?
Закинул сюда
https://github.com/VolodymyrKoliadenko/FDTD-2threads/tree/master/Matlab
TE TM версии и матричную.
Но там весьма неаккуратно.
Если говорить об оптимизации распараллеливанием, то логично обратится к OpenCL, да и готовые реализации уже существуют.
Кроме того, что в OpenCL мне разобраться показалось намного сложнее, чем в Куде,
в первом же результате:
«Although the portability resulted in a slightly reduced performance (10–35%) of the OpenCL-accelerated FDTD simulations compared with… Open Multiprocessing implementations»
Я сейчас занимаюсь подобной задачей, только для акустики. OpenCL показался мне предпочтительнее CUDA, потому что даже без GPU он способен распараллелить задачу по ядрам.
(а взялся за самостоятельную реализацию, потому что интересует решение на гексагональной сетке).
Владимир, очень странные вещи вы пишите в адрес производительности MATLAB.

Дело в том, что MATLAB, в основе своей, работает на коммерческих библиотеках, таких как BLAS, Lapack и т.д. Кроме того, производительность и оптимальность работы алгоритмов проверена временем и реальными проектами.

Возможно, некоторые алгоритмы в MATLAB будут работать медленнее чем С++, Java и т.д. Например пузырьковый алгоритм работает в 2 раза медленнее в MATLAB по сравнению с С++ (а не в 6 раз, как вы пишите). Но это исключение, а не правило.

Прежде всего, надо учитывать, что MATLAB это матричная лаборатория (MATrix LABaratory). Поэтому, максимально эффективно использовать его при векторизованных операциях. Та же векторизованная (встроенная) функция сортировки sort работает в 20 раз быстрее чем аналог C++ (библиотека boost).

Еще несколько тестов MATLAB / C++ (библиотека boost):


Размер Time MATLAB Time С++ C++/MATLAB
Умножение матриц 500х500 0,00794 C = A*B; 6,144 C = prod(A, B); 773,80
Умножение матрицы на вектор 10000х10000 0,033 c = A*b; 3,73 c = prod(A,b); 113,03
Транспонирование матрицы 10000х10000 0,226 B = A'; 6,44 B = trans(A); 28,50


Здесь лежит хорошая статья которая обобщает все способы ускорения вычисления в MALAB. Возможно она поможет вам разобраться с узкими местами кода и получить необходимый результат.
В комментариях таблицы прикладываются криво, поэтому размещу картинкой:
2017_02_27_18_42_41
«Да, я писал. Матлаб очень быстро умножает матрицы и с нормальной скоростью берет функции от массивов.»
Матрицы — намного быстрее наивной реализации на Джаве.
Функции — примерно со скоростью Джавы.
Речь о том, что Матлаб ужасно медленно крутит сложные самописные циклы.
Зарегистрируйтесь на Хабре, чтобы оставить комментарий

Публикации