Fortran или C. Что быстрее?

Введение

Не так давно я нашел исходные коды программ по вычислительной гидродинамике, написанных на языке Fortran 77. Я решил эти коды изучить на предмет полезности и перевести их на C, потому что на фортране ранее не писал и привык уже к сишнику. К тому же можно, наконец, разобраться в правдивости слухов, бытующих в научной среде, о том что программы, написанные на фортране, быстрее и точнее.

Перевод самой простой программы, решающей одномерную стационарную задачу конвекции-диффузии, особого труда не составил и занял не более часа. Самое интересно началось дальше, когда я собрал обе программы компиляторами gcc и gfortran с опцией оптимизации -O3 и сравнил результаты.

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

В качестве основных объектов исследования я выбрал:

1. Компиляторы GNU: gcc, gfortran. Этот выбор обусловлен не только моим неравнодушием в Linux, но и тем фактом, что большинстве вычислительных систем управляются Unix-подобными ОС, для которых эти компиляторы являются стандартными.
2. Intel-овские компиляторы: icc, ifort. Поскольку на моем ноутбуке стоит процессор Intel Core i7 3630QM, то любопытно будет узнать на что они способны.

Основные тесты я проводил на своем ноутбуке:
Процессор Intel Core i7 3630QM
Память: 8 ГБ (DDR3 1600 МГц / 4ГБ*2)
ОС: Debian 7.4

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

Опции оптимизации

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

При запуске в консоли gcc или gfrotran с параметром --help=optimizers увидим следующее

-O<number>  Set optimization level to <number> (Установить уровень оптимизации равным <number>)

Intel-овские компиляторы дают более развернутый ответ на запрос о помощи. При запуске в консоли icc или ifort с параметром --help увидим следующее

-O1  optimize for maximum speed, but disable some optimizations which increase code size for a small speed benefit (оптимизировать для достижения максимальной скорости, но с отключением некоторых оптимизаций, которые увеличивают размер кода для получения небольшого прироста производительности)

-O2  optimize for maximum speed (DEFAULT) (оптимизировать для достижения максимальной скорости)

-O3  optimize for maximum speed and enable more aggressive optimizations that may not improve performance on some programs (оптимизировать для достижения максимальной скорости и разрешить более агрессивную оптимизацию, которая может не увеличить производительность некоторых программ)

Таким образом чем выше цифра в параметре, тем более высокого быстродействия можно ожидать от программы.

Решение обыкновенного дифференциального уравнения 2 порядка

Рассмотрим задачу Коши вида

 

Точное решение имеет вид:

При численном решении подобные задачи приводят к систему ОДУ 1 порядка, в нашем случае такая система примет вид:

Как видно из математической постановки, необходимо определить две неизвестных функции x и y.

Теперь обратимся к программной реализации. Здесь можно предложить 2 варианта:

1. Нас интересуют значения неизвестных функций в конкретный момент времени, до которого следует проводить вычисления. В этом случае достаточно иметь две переменных для хранения значений функций в текущий момент времени и несколько буферных переменных для хранения промежуточных значений.
2. Нас интересуют значения неизвестных функций на протяжении всего исследуемого промежутка. При этом необходимы массивы для хранения сеточных функций x и y и несколько буферных переменных для хранения промежуточных значений.

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

Для замера времени я использовал функцию omp_get_wtime(). Интегрирование проводилось на сетке в 100000000 отрезков.

На следующих графиках представлены результаты измеремения времени, затраченного на вычисления.

Из диаграм видно:

1. Вариант кода без использования массивов быстрее варианта с массивами.
2. Приложения, собранные intel-компиляторами, быстрее своих GNU собратьев. При опциях -O2 и -O3 Intel-овские компиляторы обеспечили выигрыш во времени исполнения в 2 раза.
3. При максимальном уровне оптимизации O3 производительность C/C++ и Fortran приложений практически одинаковая.

Следует отметить, что результаты расчетов интеловскими приложениями отличаются  от GNU на величину не превосходящую 0.001 %. Этот факт может стать темой для отдельного исследования.

Исходные коды доступны по ссылке

Решение одномерной краевой задачи

Рассмотрим краевую задачу вида

Точное решение задачи находится простым интегрированием

Конечно-разностный аналог дифференциального уравнения имеет вид

а полученная система решается методом прогонки

Подробную информацию о прогонке можно найти в книгах по численным методам, например:

1. Вержбицкий В.М. Основы численных методов
2. Флетчер К. Вычислительные методы в динамике жидкостей

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

Для тестовых вычислений использовалась сетка с 40000001 узлом. На следующем графике представлены результаты измеремения времени, затраченного на вычисления.

Из диаграмы видно:

1. Как и в предыдущем примере приложения, собранные intel-компиляторами, быстрее своих GNU собратьев. Обеспечивается выигрыш во времени исполнения в 1.5 раза.
2. При опции оптимизации -O3 Fortran-приложения работают быстрее своих сишных аналогов.

Исходные коды доступны по ссылке

Решение уравнения Пуассона в прямоугольной области

Рассмотрим краевую задачу вида

Точное решение известно

Дифференциальный оператор аппроксимируется конечной разностью

Из конечно-разностной аппроксимации видно, что каждому элементу сеточной функции и разностному коэффициенты поставлена в соответствие пара индексов. Следовательно, для представления сеточных функций и хранения коэффициентов удобно применять двумерные массивы. Здесь я выбрал динамические массивы.

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

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

Описанную выше систему я решал методом SIP (Strongly Implicit Procedure), предложенным Stone'ом в статье Stone H.L. Iterative solution of implicit approximations of multidimensional partial differential equations // SIAM Journal on Numerical Analysis. 1968. 5(3), 530 - 558. Согласно этому методу на каждой итерации решение уточняется по формуле

а поправки определяются из выражения

коэффициенты которого находятся из формул

 

 

Расчетная сетка состояла из 201 узла в обоих координатных направлениях. Итерационный процесс продолжался до достижение равномерной нормой невязки величины 10^-8. На следующем графике представлены результаты измеремения времени, затраченного на вычисления.

Из диаграмы видно:

1. Независимо от опции оптимизации приложения, собранные gfortran, работают очень медленно. Можно предположить, что это связано с двумерными динамическими массивами.
2. Наиболее эффективно работает приложения, собранные ifort с опциями O2 и O3. Быстродействие программы, собранной icc с опцией O3, имеет соизмеримую производительность.
3. Приложение, собранные gcc с опциями O2 и O3, также обладают неплохой производительностью.

Исходные коды и скрипты сборки доступны по ссылке

Послесловие

Из проведенных тестов видно, что при должном выборе языка программирования и компилятора можно добиться максимальной производительности вычислений. Причем именно выбор компилятора может оказать наибольшее влияние на скорость вычислений.

Открытым остается вопрос о подобном тестировании рассмотренных компиляторов, на других процессорах.

В отдельном исследовании имеет смысл сравнить приведенные выше компиляторы с компиляторами Microsoft. Небольшие эксперименты с методом Эйлера-Коши и Visual C++ оказались достаточно любопытными. На разных машинах получены различные результаты, где-то занимал второе место после intel, а где-то уступал gcc.

Стоит также упомянуть про .NET. Эта платформа создавалась для эффективной разработки бизнес-приложений, поэтому не стоит ждать высокой производительности вычислений, например, от C# программ. Это предположение подтверждается простым экспериментом с методом Эйлера-Коши.