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# программ. Это предположение подтверждается простым экспериментом с методом Эйлера-Коши.

Точность машинных вычислений. Часть 2.

Введение

Возвращаемся к теме точности машинных вычислений. В первой статье

http://moderncomputingtechnology.blogspot.ru/2014/02/ComputingPrecision.html.

я рассматривал три типа данных, описывающих вещественные числа: float, double, __float128.

На примере вычисления числа e было показано, что при использовании типа __float128 затраты процессорного времени на несколько порядков выше нежели для типов float и double. Продемонстрировано, что тип float не пригоден для вычисления конечных разностей. Тип double более чувствителен к ошибке округления по сравнению с __float128.

В данной статье я рассматриваю методы интегрирования систем обыкновенных дифференциальных уравнений (ОДУ), реализованные с помощью типов float, double, long double и __float128. В механике очень часто изучаются колебательные системы, поэтому все тесты я проводил на системах с периодическим решением.

Вид уравнений

Общий вид системы ОДУ первого порядка, разрешенной относительно производных, в векторной форме имеет вид

где x - n-мерная вектор-функция одной переменной, F - заданная n-мерная вектор-функция n+1 переменной. Система, правая часть которой явно зависит от переменной t, называется неавтономной.

Если функция правой части независит от t, то систему называют автономной, и она может быть записана в виде

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

Стоит отметить, что для исследования таких систем разрабатываются различные программные пакеты. Например, в Ижевском институте компьютерных исследований разрабатывается пакет "Компьютерная динамика: хаос" для исследования динамических систем, демо-версию которого можно найти по адресу http://lab.ics.org.ru/lab/page/kompyuternaya-dinamika/.

Методы решения

Для решения ОДУ разработано множество методов. Здесь я рассмотрю четыре явных одношаговых метода: метод Эйлера, метод Эйлера-Коши, метод Рунге-Кутты 4 порядка, метод Вернера.

Метод Эйлера является самым простым и исторически самым первым методом решения ОДУ, был предложен бессмертным механиком, физиком, математиком и астрономом в одном лице Леонардом Эйлером. Метод Эйлера обладает первым порядком точности и записывается следующим образом:

Метод Эйлера-Коши является усовершенствованием предыдущего, обладает вторым порядком точности и состоит из двух этапов:

Метод Рунге-Кутты 4 порядка - наиболее популярный явный метод интегрирования ОДУ и систем ОДУ, записывается в виде

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

Метод Вернера (DVERK) относится к семейству методов Рунге-Кутты и обладает 6-м порядком точности. Уравнения метода здесь не приводятся из-за громоздкости. Описание метода и его таблицу коэффициентов можно найти в книге Hairer, Ernst, Syvert Paul Nørsett, and Gerhard Wanner. Solving ordinary differential equations. Vol. 1. Berlin: springer-Verlag, 1991. или попросить коэффициенты у меня.

При описании методов я несколько раз использовал словосочетание "порядок точности". Если точное решение является полиномом степени p, то метод p-го порядка точности теоретически даст точное решение, а решение по методу порядка p-1 будет неточным. К тому же, при уменьшении шага интегрирования, например, в 10 раз, ошибка решения, опять же теоретически, будет уменьшаться в 10^p раз.

Тестирование на автономной системе

Автономная система вида

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

Начальные условия были выбраны следующими:

Значениях параметров:

Интегрирование проводилось на отрезке

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

выбор типа данных не влияет результат при шагах интегрирования до 10^-5. Однако, при шаге 10^-6 реализация с типом float начинает терять точность.

Метода Эйлера-Коши (второй порядок точности) ведет себя аналогично методу Эйлера.

Однако float-реализация теряет при шаге 10^-4, использование типов double, long double, __float128 приводит к одинаковым результатам.

Существенные различия в представлении вещественных чисел проявляются при интегрировании системы ОДУ методом Рунге-Кутты 4 порядка

Наименее точный тип float не дает заявленной методом сходимости. Тип double обеспечивает 4 порядок до шага 10^-3 включительно, тип long double - до шага 10^-4 включительно, тип __float128 - для всех представленных шагов интегрирования.

Метода Вернера 6 порядка дает результаты на несколько порядков более точные, чем метод Рунге-Кутты 4 порядка, но и теряет порядок точности быстрее.

float-реализация вообще не дает сеточной сходимости. double- и long double- реализации начинают терять порядок точности с шага 10^-3. __float128-реализация теряет точность только при шаге 10^-6.

Таким образом, видно две закономерности:

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

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

Видно, что реализации __float128 требуют на 2 порядка больше времени. Таким образом, тип __float128 нецелесообразно использовать для проведения всевозможных параметрических исследований или расчетов длительных процессов. long double требует до трех раз больше процессорного времени нежели тип double. Затраты на обработку типов float и double примерно одинаковы.

Тестирование на неавтономной системе

Теперь рассмотрим неавтономную систему

решение которой задается следующими выражениями:

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

или рекуррентно

Аналитическое определение момента времени использует операцию умножения, рекуррентное - операцию сложения. С точки зрения математики разницы нет, но в программной реализации выявляются некоторые особенности, которые демонстрируют следующие 9 графиков.

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

Полученный результат позволяют сделать два вывода:

1. При фиксированном методе. Чем точнее тип данных, тем меньше влияет дополнительная погрешность.
2. При фиксированном типе данных. Чем выше порядок точности метода, тем сильнее он чувствует погрешность.

Послесловие

Продемонстрированные в данной статье примеры численного решения показывают, что, используя численные методы (не обязательно методы интегрирования ОДУ), нужно понимать их свойства и ограничения, накладываемые вычислительными системами. Даже простая формула определения момента времени может существенно повлиять на точность решения.

О числах с плавающей точкой можно найти краткую информацию по ссылке.

Если Вас заинтересовали реализации методов, пишите: veugene186@gmail.com

Множество кодов по методам решения ОДУ можно найти по ссылке http://www.unige.ch/~hairer/software.html

Автор выражает благодарность к. ф.-м. н. А.О. Казаков за полезные обсуждения и ценные замечания.