Комментарии 104
Прирост от «4»-ядерного кода на С++ также невелик — 67% между 1 и 4 ядрами (оба float). Списывать можно на медленный контролер памяти и 32-битную Винду.
На самом деле списывать надо на то, что вы не специались по распараллеливанию, как и не специалист по оптимизации...
Кстати, у вас в коде на C++ есть несколько new, но ни одного delete...
Кстати, у вас в коде на C++ есть несколько new, но ни одного delete...Как ни странно, но в данном случае освобождение памяти средствами приложения вместо средств ОС только немного замедлит программу. ;)
На самом деле списывать надо на то, что вы не специались по распараллеливанию, как и не специалист по оптимизации...О да, Intel VTune Amplifier спасет отца русской демократии:
https://habrahabr.ru/company/intel/blog/310842/
Судя по ограничению скоростью памяти, вам стоит попробовать учесть иерархию памяти, чтобы повысить локальность обращений к памяти (процент попадания в кэш).
https://habrahabr.ru/post/312078/
Далее -march=native, векторизация (-sse2, avx), интринсики, профайлеры…
После чего становится грустно и взоры обращаются обратно к ANSYS ;)
-march=native проверю.
Извините, я компилирую на одной, а проверяю на разных.
Использовался вроде ‘nehalem’, проверял что от ‘sandybridge’ пользы нет (проверю, по моему у меня опции не так называются)
В интринсиках разбираться не планирую
Не вижу возможности вынести из циклов. Все зависит от [i][j]
да, да конечно. Я понимаю, что это не самый внутренний цикл, но сравнение строк… Будем молиться, что компилятор сам раздвоит циклы для разных case (о чем бы было неплохо ему напомнить явно). Уж лучше лямбды/указатели на ф-ии, виртуальные методы в явном виде. Накрайняк — глобальный копипаст (как обычно, это всего быстрее:) ).for (int t = 1; t <= tMax; t++) { // begin main loop float tt = Math.min(t * s + 10, ny - 1); //gauss switch (method) { case "cos":
в
я правильно понял, что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 части, имеем небольшую потерю производительности.
Если Вы упираетесь в скорость памяти — попробуйте разбить задачу на кусочки по-другому. Например, можно взять область (100+2x, 100+2x), и на её основе посчитать вперёд сразу на х шагов центральный кусочек размером 100*100. В идеале должно получаться, чтобы необходимая для "кусочка" память была не слишком большой и помещалась в кеш процессора — тогда у Вас на первом шаге будет чтение из ram и на следующих (х-1) шагах будет использоваться кеш.
Можно только попытаться разбить циклы внутри одного шага.
нет, если взять х небольшим, то должно работать быстрее.
Например, для x=5 Вы рассчитаете квадраты размерами 110, 108,… 102. Конечно,
это немного больше вычислений по сравнению с вычислениями и чтением данных из 5 квадратов 102*102, (процентов на ~20%), но нагрузка на оперативную память упадёт почти в 5 раз — а именно она и является бутылочным горлышком.
Даст ли выигрыш переписывание и оптимизация на него?
На Фортране пишут потому, что под него есть тонны математических/научных библиотек
https://benchmarksgame.alioth.debian.org/u64q/fortran.html
Фортран или на уровне чистого С, или медленнее.
Там ФДТД используют на 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];
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];
Отличная оценка, спасибо за неудачную идею.
В таком способе нужно читать эпсилон 3 раза.
Вместо 5 чтений 7+лишение записи. Плюс лишняя математика.
А в Н такого и не сделаешь
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 в циклах? И не совсем понял, вы какой профилировщик используете в С++?
2. Да, компиляторы нынче умные
Профилировщик в С++ никакой, хотя вроде в VisualStudio смотрел.
Как идея, сделать все данные более локальными. Например:
1. Сделать одну матрицу вместо 4-х где хранить коэффициенты друг за другом или как структуру.
2. Сделать одномерный массив и руками интерировать по строкам.
1-мерный массив уменьшает производительность.
1. А вот хранить все локально — интересный вариант. Только все циклы придется ужасно переписывать.
Ради прикола, попробовал заменить в одномерной версии float* на char* (понимаю что считать будет не правильно) — выигрыш всего 25%. Плюс, 700 * sizeof(char) * 1024 = 700кб — это смешно по памяти (не похоже на память). Значит вам нужна векторизация, он просто не успевает считать.
К сожалению это специфично для компилятора и сами матрицы должны быть выровнены. Быстро не проверить.
Почему же был высокий прирост от 64бита к 32?
1024*700 — это тестовая задача, какую смог дождаться в Матлабе. С++ детально не изучал, в Джаве наибольший упор в скорость памяти при 4096х500
Все выглядит так, как-будто циклы просто очень медленно считаются, хотя тела маленькие и операции элементарны. Здесь хорошо бы помог VTune Analizer или что-то подобное. У меня его к сожалению нет, а профайлер от Microsoft очень грубый, он показывает сколько CPU затрачивает строчка кода, но при этом нельзя посмотреть ассемблерный код. VTune в таких случаях показывает более точную информацию. Может быть у вас где-то зависимость по данным и вся параллелизация убивается на корню.
Много тратит запись элемента массива. Если принять время выполнения 3 циклов за 100у.е., то
каждая запись(+=) из 3х тратит 20,
запись на даблах 31 (просто для сравнения)
однократное чтение 6-7
пятикратное чтение 27
(оценки приблизительны, внимательный читатель поймет, что мы уже ниже 0)
«Согласно протоколам согласованности кеша (cache coherence protocols), если ядро пишет в кеш-строку, то строка в другом ядре, ссылающаяся на ту же память, признаётся недействительной (пробуксовка кеша, cache trashing). В результате при каждой операции записи возникают блокировки памяти.»
У меня такое возможно на середине массива. То есть задача влазит в кеш, но пишется с проблемами.
P.S. У вас нет случайно ссылки на ваш проект под студию, а-то я нечаянно похерил свой?
Под студию был 1-ядерный.
Постараюсь сегодня сделать 2-4.
https://cloud.mail.ru/public/8K8G/VCRZSJQTt
Есть профит по сравнению с ГСС во всех случаях.
Скоро обновлю в статье
vlanko — в С++ код при смешивании типов приводится к более общему типу, т.е. j приводится к float
Касаемо кода с++, по беглому просмотру, он конечно нуждается в оптимизации и, главное, в векторизации. Можно воспользоваться готовыми библиотеками с векторными классами, например 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
Скажите нубу, что делает это
asm("# 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];
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 раза
Я же посмотрел визуализацию на Яве. Фигня была именно в том, что обрабатывается только первая четверть картинки.
В вашем коде:
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];
Этого я не сделал в вашем алгоритме, о чем сразу же и написал:
Разумеется, граничные условия нужно дообрабатывать отдельно или указывать границы кратные количеству развертываний цикла. Но это мелочи и на производительность практически никак не скажется.
Советую вам корректно развернуть циклы, так как потенциально это может дать сильный прирост производительности на больших итерациях
Но развертывание цикла, как я показал выше, вполне имеет место быть, у меня время просчета сократилось с 4500мс до 1900мс
Может быть возможно изменить алгоритм так, чтобы данные не пересекались — тогда можно будет векторизовать.
С ними так же просто, как в шейдерном коде. Если вам не удастся векторизовать ваш алгоритм, то и на CUDA особого прироста не ждите.
Проверить скорость чистого С (не ++). По benchmarksgame он всегда быстрее.
У вас код и так написан на почти чистом Си. Так что никакого прироста не получится.
Я понимаю что писать свой FDTD это занимательно (лично знаю человека который из принципа написал свой 2D FDTD на CUDA), но для исследовательских целей почему бы не использовать готовый продукт и не строить свой велосипед?
И еще по поводу Матлаба. Циклы и условные ветвления (а так же вызовы анонимных функций) в нем действительно медленные, но вот встроенные алгоритмы работы с массивами (сложение, умножение матриц, нахождение определителей, собственных значений и т.д.) давольно быстры. На мой первый взгляд в матлабе возможно реализовать FDTD используя только лишь один внешний цикл для времени, все остальное можно сделать используя матричные операции.
Да, я писал. Матлаб очень быстро умножает матрицы и с нормальной скоростью берет функции от массивов. Возможно, получится такое сделать. Но это потребует создание 4 новых матриц для вычитания. Боюсь, со скоростью памяти будет хуже.
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х.
В любом случае, спасибо за идею.
Если вы в лс опишите структуру и скажете что посчитать, я могу сделать .ctl файл для meep и прогнать вычисления.
Видимо автором ожидается, что цикл по t ожидает завершения стартованных по join() нитей, но это не так.
— запуск 2х потоков
— ожидание их завершения
— запуск еще 2х потоков
— ожидание их завершения
Какое может быть ускорение? Это обычный 2-поточный код
— запуск 2х потоков
— в каждом из двух потоков еще запуск по 2 потока
То есть два треда, в каждом запускаются сотни раз еще по паре тредов, в итоге одновременно запущено в среднем 4 потока.
В mingw действительно что то не так в потоках, некоторые версии mingw с потоками не работают совсем, некоторые абы как. В общем, принудительно поставив пару mutex'ов, заставил отработать как надо.
На результат не повлияло.
Простая замена на линейный массив [y + maxx *x] привела к большому падению скорости.
Дальше надо лезть в ассемблер.
И главный вопрос — почему при распараллеливании не растет пропорционально скорость?
FFT хорошая задача, еще в каком то старом тесте была
Поскольку в быстром Фурье не разбираюсь, взял первый попавшийся. Минус — ширина только степень двойки.
Это стандартная особенность БПФ.
Вообще, последние i7 имеют от десяти до двадцати мег кэша L3 — хотя он тормознутее, конечно — но ограничение скорости по памяти должно быть слабее, чем для других.
А вообще, получается необходимо дробить алгоритм на отдельные «островки», хорошо попадающие в кэш.
Да, но у моих процей только 6-8 МБ кеша.
Малые задачи у меня выполняются быстрее — эквивалент 30 ГБ/сек против 17 ГБ/сек, (память 20около).
То есть лезут в кеш.
Но тупое деление большой области на 4 куска по і делает немного хуже.
Буду проверять развернутые циклы.
В Матлабе много чего наворочено «внутри» — но, как всегда в таких случаях — вы берёте готовое.
И «под капот» вас не пускают.
Кстати, дурацкая идея — я гляжу, что во всех научных расчётах есть тенденция разбивать всё на составляющие — Ez, Ex, Ey и т.д.
А если поступить наоборот?
Ну т.е. у нас фиксируется пространство «точек», в каждой точке имеется набор свойств — запись либо массив, Ex, Ey, Ez, Hx, Hy, Hz — и, возможно — локальные производные.
Тогда локальные расчёты автоматически локализуются ;)) — для двумерного случая у нас будет три «слайса» из массива по три элемента в каждом. Или даже два.
Возможно алгоритм будет работать «компактнее».
Intel VTune Amplifier в https://habrahabr.ru/company/intel/blog/310842/ предлагает обратно разбить на несколько массивов
А то под МКЭ существуют свои пакеты — от коммерческих до бесплатной Salome.
Посмотрев ассемблерный код на отличном сайте https://godbolt.org, выяснилось, что gcc в 32-битной версии не использует SSE и подобное.
В итоге нужно брать нормальную 64-битную версию и компилировать примерно так
>-O2 -lstdc++ -march=core-avx2 -m64
Для матрицы 4906х1000 у меня результат изменился со 154с на 29.5с, т.е грубо в 5(пять) раз.
Исходник автора.
Т.е так >-lstdc++ -O2 -march=core-avx2 -m32 -o sie32.exe -mfpmath=sse
Это тоже решает проблему
Включение опции -mfpmath=sse не изменяет ничего.
Полный набор ваших опций увеличивает время на 5%
Если что, мои опции mingw32-g++.exe -Wall -fexceptions -O2 -march=corei7-avx -O3 -std=c++11
Буду пытаться 64-битный компилятор. Под винду сложно найти нормальное.
Это загрузчик. Нужно брать вариант с posix-threads. Оба варианта версии 6.2 i686 и x86-64 работают корректно.
-O3 не рекомендую использовать — глючная
И добавив в конце 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 << " ";
Чтоб я мог сравнить скорость программы и результат
Там исходник, необходимый рантайм, батник для компиляции под AVX, и две версии — SIMD и обычная O2 — см батник
После установки гсс оно появится?
Конечно. Рантайм GCC лежит в папке bin компилятора. К сожалению с ним затруднительно слинковаться статически.
i7-4771
Моя версия выполняется 22,3 секунд
ваша в моем компиляторе 22,4 (конечно, это тоже самое)
Ваш АВХ 23,2
Ваша без инструкций 93,5 секунд
Для этого процессора архитектура -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
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? Желательно с юнит-тестами… Простор для оптимизации тут явно есть. Ну и параметры запуска надо бы озвучить.
Нужно смотреть конкретные уравнения. Возможно, в акустике будет проще.
По сути, это уравнения Максвелла — то, что до волн. Если нужны просто волны — там более быстрые методы.
в первом же результате:
«Although the portability resulted in a slightly reduced performance (10–35%) of the OpenCL-accelerated FDTD simulations compared with… Open Multiprocessing implementations»
Дело в том, что 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. Возможно она поможет вам разобраться с узкими местами кода и получить необходимый результат.
Простая реализация FDTD на Java