В данной заметке я расскажу о своем опыте распараллеливания программы с использованием OpenMP, написанной одновременно и на С++, и на Fortran 90, причем вызов фортрановской части кода осуществляется параллельно в цикле из C++ части.  Остановлюсь в основном на тех деталях и тонкостях, которые мне показались настоящими сюрпризами. Суть программы достаточно проста: есть некая цилиндрическая структура, которая модельно разбивается на аксиальные ячейки, каждая аксиальная ячейка с использованием методов математического моделирования обсчитывается независимо. Основная часть кода написана на C++, но вот то, что нужно рассчитать для каждой аксиальной ячейки, написано на Fortran 90, и надо сказать, что эта фортрановская часть достаточно внушительная. Код испокон веков обсчитывал последовательно каждую аксиальную ячейку, и в виду того, что фортрановская часть делает объемные вычисления, код считал долго. И тут была поставлена задача – распараллелить код, т.е. считать каждую аксиальную ячейку параллельно, дабы ускорить время расчета всей программы. Была принята следующая идея: выделение памяти под массивы и т.п., необходимые для фортрановских расчетов, оставить как и прежде, т.е. там же в фортране, удобно перегруппировав их в массив объектов структуры, описывающей аксиальную ячейку, и плюсом расширив данные, а вот вызов главной функции, делающей вычисления в фортрановской части, делать в цикле по аксиальным ячейкам параллельно в C++ части, т.е. примерно так

# pragma omp parallel for
for (iaxialmesh = 0; iaxialmesh<Nmeshes; iaxialmesh++) {
	CallFortranCodeForAxMesh(iaxialmesh);
}

Здесь опущены все детали с private и т.п., необходимые в строке c # pragma. Таким образом, мы вызываем параллельно из C++ части функцию, которая написана на фортране. Эта фортрановская функция в свою очередь зовет массу других функций, в общей сложности которых несколько сотен.

Итак, главной модификацией фортрановской части было избавление от глобальных переменных, которые в большом количестве были помещены в модули. Именно в такие моменты приходит понимание, что глобальные переменные – мины замедленного действия. Фортрановская часть программы до этого работала на глобальных переменных, каждый раз перезаписывая эти глобальные переменные перед выполнением расчета для текущей аксиальной ячейки, а затем выкачивала нужные насчитанные данные. Теперь же была создана структура tAxmesh. Массив объектов данной структуры tAxmesh axmeshes(NtotAxMesh) благополучно хранит все данные для каждой аксиальной ячейки изолированно в памяти (разумеется, объем потребляемой программой оперативной памяти увеличился) и позволяет выполнять вычисления для каждой аксиальной ячейки параллельно. Дело в том, что многие массивы нужны только для выполнения расчета на текущем временном шаге, потом их можно без проблем затереть. Назовем эти массивы «рабочими». Как раз они-то и лежали в глобальной области. Теперь, если мы хотим проводить расчет сразу для всех ячеек, эти рабочие массивы должны быть созданы для каждой аксиальной ячейки, потому и вырос объем потребляемой памяти.

Теперь встает вопрос, как передавать во внутренние фортрановские функции информацию о местонахождении в памяти массивов, относящихся к той или иной аксиальной ячейке. Как вариант, можно передавать только один параметр - объект структуры для текущей ячейки axmeshes(iCurrentAxMesh).  Но тогда в большом количестве мест фортрановского кода пришлось бы дописывать к переменным somevariable приставку axmeshes(iCurrentAxMesh), т.е. axmeshes(iCurrentAxMesh)%somevariable. Для того чтобы это не делать, было принято решение передавать напрямую нужные переменные через интерфейсы. Объем интерфейсов функций в фортрановской части программы в результате этого изрядно подрос, т.к. раньше функции просто обращались за переменные в модули use SomeModule, теперь же нужно передавать указатели на массивы в тех местах памяти, которые соответствуют данной обсчитываемой аксиальной ячейке.

Первое опасение заключалось в том, что нет ли каких либо лимитов по число параметров, передаваемых в интерфейсе функции в Fortran 90. Опасения не оправдались, можно ставить хоть несколько сотен, проблем нет. Второе опасение заключалось в том, что не вылезем ли мы за пределы памяти, доступной для текущего вычислительного потока. Здесь также, проблем не возникло, памяти, используемой потоком по умолчанию, с лихвой хватало для вычислений достаточно объемной программы.

Настоящим сюрпризом оказалась следующая деталь. Когда я захожу в управляющую функцию в фортране с определенным id аксиальной ячейки CallFortranCodeForAxMesh(iaxialmesh), далее мне необходимо передать информацию о месте памяти, где лежит все относящееся к данной ячейке. Логично сделать это через указатель, но такой код работает некорректно!

! код, работающий некорректно
subroutine CallFortranCodeForAxMesh(iaxialmesh) bind (C, name = “CallFortranCodeForAxMesh”)
use AxMeshModule ! модуль, где дано описание структуры tAxmesh
integer(c_int), intent(in) :: iaxialmesh
type (tAxmesh), pointer :: paxmesh => NULL()
paxmesh => axmeshes(iaxialmesh)
call AxMeshComputation(paxmexh)
end CallFortranCodeForAxMesh

Проблемы с такой реализацией начались с того, что компилятор начал ругаться на повторную деаллокацию массивов во внутренних функциях фортрана. Речь идет о таких массивах, которые объявляются, аллоцируются и деаллоцируются внутри какой-либо подпрограммы для совершенно локальных задач. В корректности кода не было никаких сомнений. В итоге удалось выяснить, что при такой реализации два различных потока начинают работать над одной и той же областью памяти, относящейся к одной и той же аксиальной ячейки. Почему это именно так, я так и не понял, приписав это к багу фортрановского компилятора при параллельной работе кода.

Эта проблема была решена прямой передачей объекта массива:

! корректный код
subroutine CallFortranCodeForAxMesh(iaxialmesh) bind (C, name = “CallFortranCodeForAxMesh”)
use AxMeshModule ! модуль, где дано описание структуры tAxmesh
integer(c_int), intent(in) :: iaxialmesh
call AxMeshComputation(axmeshes(iaxialmesh))
end CallFortranCodeForAxMesh

С такой реализацией потоки разбирают четко отдельные области памяти, относящиеся к различным аксиальным ячейкам.

Следует отметить, что когда мы находимся во внутренних участках фортрановского кода и считаем конкретную аксиальную ячейку, уже находясь в изолированной области памяти, относящейся к заданной аксиальной ячейке, здесь уже передача по указателю, каких-либо объектов частных структур axmeshes(iaxialmesh)%someStructureObj не вызывает проблем.

Следующий сюрприз поджидал при инициализации данных. Для того чтобы зачитать входной файл был создан массив объектов некоторой частной структуры structObj, лежащий в глобальной области, а затем эти данные раскидывались по соответствующим местам, относящимся к различным аксиальным ячейкам axmeshes(i)%structObj. Изначально это было так:

! копируются только указатели!
axmeshes(i)%structObj(j) = structObj(j)

Компилятор вновь начал выдавать конфликт по памяти при параллельной работе потоков, та самая гонка, или race condition. Оказалось, что данное копирование происходит только для указателей, т.е. указатель на массив structObj(j) был скопирован в объект каждой аксиальной ячейки. Разумеется, все потоки начинали работать с одним и тем же местом в памяти. Для того, чтобы это разрешить, необходимо полностью прописать то, что мы хотим скопировать, а именно содержимое структур, т.е.:

! копируем содержимое
axmeshes(i)%structObj(j)%var1 = structObj(j) )%var1
axmeshes(i)%structObj(j)%var2 = structObj(j) )%var1
axmeshes(i)%structObj(j)%var3 = structObj(j) )%var1
axmeshes(i)%structObj(j)%var4 = structObj(j) )%var1

Сколько бы переменных не было, их нужно перечислить все.

Сборка программы на Linux с компиляторами gcc и gfortran версии 11 не вызывала никаких проблем. Но вот на Windows с использованием Visual Studio 2017 и компилятором фортран Intel Parallel Studio XE 2019 Update5 Composer код собирался (с установкой всех нужных флагов подключения OpenMP в настройках), но не считал. Причина оказалась в конфликте версий OpenMP: в компиляторе C++ в Visual Studio 2017 стоит версия OpenMP 2.0, в упомянутом компиляторе фортрана - OpenMP 5.0. Однако, если компилятор Parallel Studio XE 2019 Update5 Composer заменить на компилятор фортрана на OneApi 2022, то проблема конфликта решается, - код считает параллельно.

Итак основные выводы:

• при параллельном запуске фортрановской функции следует в интерфейсе напрямую указывать объект структуры, а не указатель на этот объект

• приравнивание объекта структуры другому объекту структуры в фортране делает копирование указателей, а не содержимого структур

• программа, написанная на C++ и Fortran 90 не будет работать в Visual Studio 2017 с компилятором фортрана Parallel Studio XE 2019 - имеет место быть конфликт версий OpenMP в компиляторах разных языков. Следует ставить либо Visual Studio 2017 + компилятор фортран OneApi 2022, или Visual Studio 2019 и выше с компилятором фортрана OneApi 2024.