четверг, 2 марта 2017 г.

Подключение библиотек Intel MKL или OpenBLAS к GFortrаn

Раньше я был уверен, что, возможно, самая крутая в мире :) математическая библиотека Intel MKL – полностью коммерческая и с чистой совестью пользоваться ей можно только купив лицензию у Intel. Однако, разбираясь со свободными (совершенно бесплатными) системами для научных вычислений, я обнаружил, что Intel MKL входит в состав некоторых из них и используется там как основной решатель для СЛАУ, а также других задач. Как минимум, она входит:

  • в библиотеки NumPy и SciPy в составе дистрибутива Anaconda Python;

  • в SciLAB (известная MatLAB-подобная свободная система компьютерной математики).

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

Intel License for Intel® Math Kernel Library and Installer Program (version January 2017). Copyright (c) 2017 Intel Corporation.

Use and Redistribution. You may use and redistribute the software (the “Software”), without modification, provided the following conditions are met:

  • Redistributions must reproduce the above copyright notice and the following terms of use in the Software and in the documentation and/or other materials provided with the distribution.
  • Neither the name of Intel nor the names of its suppliers may be used to endorse or promote products derived from this Software without specific prior written permission.
  • No reverse engineering, decompilation, or disassembly of this Software is permitted.

. . .

Я не юрист, но из этого по-моему следует, что я имею полное право пользоваться MKL и даже ее распространять, если уведомлю в документации на свою программу, что копирайт MKL принадлежит Intel.

Видимо, разработчики Anaconda Python и SciLAB так и делают.

Конечно, если действительно на MKL придется завязать какие-то коммерческие проекты, этот вопрос потребует более глубокого юридического изучения. Для меня подобные юридические вопросы – полные дебри. Тут должны разбираться специальные люди. Но надеюсь, что простое экспериментирование с MKL, так же как и расчеты в NumPy и SciLAB, прав Intel не нарушают. :)

Мои «околонаучные» проекты в основном сделаны в GFortran. А Intel MKL для Windows, к сожалению, тесно интегрирована только с Intel Fortran, который уже точно коммерческий и весьма недешевый. (Под Linux MKL официально работает в том числе и с компиляторами GCC).

Подключить Intel MKL к GFortran в Windows на первый взгляд казалось нетривиальной задачей, но, после определенного разбирательства, это у меня получилось. Как – хочу зафиксировать ниже. Если потом все-таки окажется, что юридически это делать нельзя :), Intel MKL всегда можно заменить на OpenBLAS. Как он подключается – тоже запишу здесь для удобства. (Сказанное относится к библиотекам LAPACK и BLAS – больше мне из MKL пока не нужно).

Intel MKL в GFortrаn

Дистрибутив Intel MKL содержит скомпилированные библиотеки .dll, а также библиотеки .lib для компиляции и линковки в Intel Fortran и среде MS Visual Studio.

Мы не используем Visual Studio, а работаем с GFortrаn в составе компиляторов GCC. В принципе, GCC работает и с библиотеками .lib. Однако, корректно подключить все нужные для компиляции MKL .lib к GFortrаn – не получается. Судя по форумам, не у меня одного. Поэтому, чтобы вызывать из нашей программы MKL, попробуем подготовить библиотеки в родном для GCC формате .a и ограничимся динамическим связыванием. В принципе, .a содержит дистрибутив Intel MKL для Linux (я не поленился, установил MKL в Ubuntu), но .a в Ликунсе и в Виндовсе – несовместимы.

Нам надо получить из .dll и .lib библиотеки .a. Это возможно. Оказалось, что проще не использовать для этого .lib вообще. Библиотеки .a можно сделать просто из *.dll. Правда, это будут библиотеки только для динамической компоновки. Т.е. чтобы наша программа работала, в системе должно иметься несколько *.dll с Intel MKL. Кстати, если у нас инсталлирована, например, Anaconda Python, то все необходимые *.dll библиотеки с ней уже установлены.

Итак, у нас есть ряд .dll с Intel MKL. Их можно взять как из Anaconda Python, так и из официального дистрибутива Intel MKL, скачанного с их сайта. Первое, по-моему, даже проще т.к. на сайте Intel нужна регистрация.

Оказалось, что нужные нам функции LAPACK и BLAS в привычном формате вызова содержатся конкретно в файле mkl_rt.dll. Чтобы посмотреть, какие функции входят в *.dll есть разные утилиты. Проще всего воспользоваться утилитой dumpbin.exe из состава Visual Studio:

dumpbin /exports mkl_rt.dll

Она выдает следующее:


Ниже по списку упоминаются привычные функции LAPACK dgesv, dposv и другие.

Чтобы сделать библиотеку *.a из *.dll сначала надо сформировать текстовый файл .def со списком названий функций, т.е. из того, что обведено красным. В начале этого файла надо написать слово EXPORTS. Т.е. файл .def должен выглядеть так:

EXPORTS
CAXPBY
CAXPY
CAXPYI
CAXPY_DIRECT
CBBCSD
CBDSQR
CCOPY
CDOTC
CDOTCI
CDOTU
...

Файл .def можно сформировать из того, что вывел dumpbin.exe руками, а можно сделать скрипт, например, на питоне (привожу его ниже).

Получив .def, осталось сформировать библиотеку *.a (она должна иметь префикс lib и суффикс .dll) с помощью утилиты dlltool.exe из состава GCC:

dlltool -d mkl_rt.def -D mkl_rt.dll -l libmkl_rt.dll.a

Вот скрипт на Питоне, который делает все описанное. Чтобы он работал в системе должны быть доступны dumpbin.exe и dlltool.exe. Кстати, чтобы это проверить, в Windows можно набрать where dumpbin – система выведет путь, где он лежит, либо скажет, что не находит такого.

DllFile = "mkl_rt.dll"  # *.dll-файл, из которого будет сделан *.a-файл

import os

(FileBaseName, FileExt) = os.path.splitext(DllFile)
TmpFile = FileBaseName + ".tmp"
DefFile = FileBaseName + ".def"
AFile   = "lib" + FileBaseName + ".dll.a"

os.system("dumpbin /exports " + DllFile + " > " + TmpFile)

LineNo = 1
FirstLineNo = 0
LastLineNo = 0
TmpFile_ = open(TmpFile, 'r')
DefFile_ = open(DefFile, 'w')
DefFile_.write("EXPORTS\n")
for Line in TmpFile_:
    if Line.find('ordinal hint') >= 0:
        FirstLineNo = LineNo
    if Line.find('Summary') >= 0:
        LastLineNo = LineNo
    if FirstLineNo != 0 and LineNo > FirstLineNo and \
       (LastLineNo == 0 or  LineNo < LastLineNo):
        DefFile_.write(Line[26:]) # Названия функций – начиная с 26 символа
    LineNo += 1
TmpFile_.close()
DefFile_.close()

os.system("dlltool -d" + DefFile + " -D " + DllFile + " -l " + AFile)
os.system("DEL *.tmp")

Полученный файл *.a надо положить туда, где GCC хранит свои библиотеки, либо в папку с нашими собственными библиотеками, например, C:\gcc\mylibs.

В принципе, это все. Проверим как все работает.

Протестируем MKL на примере из прошлого поста с тестом решения СЛАУ.

Тестовая программа на Фортране (в скоростном варианте с \(LL^T\)-разложением) выглядит так (файл LinTest.f90):

program LinTest
    implicit none
    integer, parameter      :: n = 1000 ! Число уравнений
    integer, parameter      :: m = 100  ! Число шагов
    real(8), dimension(n,n) :: A, A0    ! Матрица коэффициентов
    real(8), dimension(n)   :: B        ! Вектор правых частей / результатов
    real(8)                 :: s        ! Контрольная сумма
    integer                 :: i, j, k  ! Вспомог. индексы
    integer                 :: info     ! Инфо об успешности решения
    integer                 :: t_0, t_1, t_rate, t_max ! Для измер. времени
    !------------------------------------------------------------------------ 
    forall(i=1:n, j=1:n) A0(i,j) = 1-(dble(i-j)/n)**2 
    forall(i=1:n, j=1:n, i==j) A0(i,j) = 10.
    s = 0.
    call system_clock(t_0, t_rate, t_max) 
    do k = 1, m
        A=A0
        B=k
        call dposv("L", n, 1, A, n, B, n, info)  ! Вызов процедуры LAPACK
        if (info/=0) stop ("Ошибка")
        s = s + sum(B)
    enddo
    call system_clock(t_1, t_rate, t_max)
    print *, "Время:", real(t_1-t_0)/t_rate
    print *, "Сумма:", s
endprogram

Условно статическая компиляция (пусть наша библиотека libmkl_rt.dll.a лежит в C:\gcc\mylibs):

gfortran LinTest.f90 -static -L C:\gcc\mylibs -lmkl_rt.dll -O2 -o LinTest.exe

Запуская программу, получаем:

LinTest.exe
 Время:  0.328000009
 Сумма:  22439.195088956156 

Время 0,33 с – точно такое же как и у лидера прошлого теста Intel Fortran + MKL.

Размер LinTest.exe получился 595 КБ (компилятор был из GCC 6.3.0 из дистрибутива MinGW64). В него статически прилинковались все библиотеки, кроме .dll-библиотек MKL.

Какие в принципе .dll нужны для работы нашей программы? Если мы хотим, чтобы наша программа работала на другой машине, не имеющей установленную Анаконду или другой вариант MKL, надо будет скопировать и обеспечить доступность 5-ти файлов:

Размер       Файл
----------------------------------
  1 328 896  libiomp5md.dll
 37 854 480  mkl_avx2.dll
 25 122 064  mkl_core.dll
 24 297 744  mkl_intel_thread.dll
 12 629 264  mkl_rt.dll
----------------------------------
101 232 448 байт

100 МБ – конечно, не мало. Это вариант с последним дистрибутивом MLK с сайта Intel, версия 2017.2.187. Из дистрибутива Анаконды набегает поменьше – «всего» 86 МБ. Современный софт, он такой. :)

Полностью динамическая компиляция (убираем ключ -static):

gfortran LinTest.f90 -L C:\gcc\mylibs -lmkl_rt.dll -O2 -o LinTest.exe

LinTest.exe получился крошечный – 57 КБ, но теперь для работы потребуются еще три .dll из GCC (данные для версии 6.3.0):

Размер      Файл
----------------------------------
    75 264  libgcc_s_seh-1.dll
 1 277 952  libgfortran-3.dll
   325 632  libquadmath-0.dll
----------------------------------
 1 678 848 байт

OpenBLAS в GFortrаn

Установить OpenBLAS – куда проще.

Скачиваем с официального сайта уже скомпилированные бинарные файлы в архиве OpenBLAS-v0.2.19-Win64-int32.zip (версия для 64-разрядных машин).

Переписываем файлы libopenblas.a и libopenblas.dll.a туда, где GCC хранит свои библиотеки, либо в папку с нашими собственными библиотеками, например, C:\gcc\mylibs.

Статическая компиляция (пусть наши библиотеки *.a лежат в C:\gcc\mylibs):

gfortran LinTest.f90 -static -L C:\gcc\mylibs -lopenblas -O2 -o LinTest.exe

Запускаем:

LinTest.exe
 Время:  0.483999997
 Сумма:   22439.195088956145 

Время на 48% больше, чем с Intel Fortran + MKL, но на самом деле тоже очень достойное.

Размер файла LinTest.exe , правда, громадный: 26 МБ.

Динамическая компиляция (убираем -static и добавляем к -lopenblas суффикс .dll):

gfortran LinTest.f90 -L C:\gcc\mylibs -lopenblas.dll -O2 -o LinTest.exe

Размер LinTest.exe – всего 57 КБ. Для работы потребуется libopenblas.dll из дистрибутива OpenBLAS и еще три DLL из GCC. Общий расклад (данные для версии 6.3.0):

Размер      Файл
----------------------------------
    75 264  libgcc_s_seh-1.dll
 1 277 952  libgfortran-3.dll
42 587 753  libopenblas.dll
   325 632  libquadmath-0.dll
----------------------------------
44 266 601 байт

«Всего» 44 МБ – не так и страшно. :)

понедельник, 27 февраля 2017 г.

Оценка скорости решения СЛАУ различными средствами

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


Введение

Значение решателей СЛАУ в современном мире невозможно переоценить. К СЛАУ сводится подавляющее большинство задач современной инженерии. При этом часто СЛАУ получаются очень большими. Зачастую, к тому же, и решаться они должны множество раз.

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

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

Для «простого» инженера вопрос о выборе средств решения СЛАУ на обычном персональном компьютере (в рамках решения насущных задач) – не менее актуален. Тема, на самом деле, сложная. Чтобы проанализировать вопрос корректно, нужно погрузиться достаточно глубоко. Надо изучить (или освежить в памяти) базовые методы решения СЛАУ, теорию и практику параллельных вычислений, архитектуру современных машин, многочисленные реализации современных решателей. Попробуем всего этого этого не делать и подойдем к вопросу «по-дилетантски». Опробуем разные доступные средства «по-бытовому»: на каждом проведем всего один тест со «взятой с потолка» СЛАУ, используя средства «из-коробки» (ну, или почти «из-коробки») и сравнивая время расчета. Сразу отдадим себе отчет в том, что такой подход мягко говоря не очень корректен:

  • эффективность расчета разными средствами может сильно зависеть от типа матрицы СЛАУ (и выбора соответствующего метода решения);

  • у разных средств может сильно различаться масштабируемость относительно размера матрицы;

  • «под капотом» различных средств могут лежать разные (или одни и те же) решатели СЛАУ, обычно это скомпилированные библиотеки. Иногда в рамках одного средства можно выбрать различные решатели, и каждый решатель тоже можно настроить по-разному (например, на работу в одно- или многопоточном режиме). Сейчас мы разбираться с этим не будем, используем все «из коробки»;

  • на разных машинах (например, с разным количеством ядер) эффективность различных средств может очень сильно разниться;

  • с решением СЛАУ в инженерных задачах часто связана основная доля потери времени работы программы, но далеко не единственная. Скорость работы остального кода на разных языках и в разной реализации тоже может различаться в разы или десятки, сотни, тысячи… :) раз.

Тем не менее, для первичной оценки такой тест провести интересно.

В качестве матрицы коэффициентов СЛАУ выберем симметрическую положительно определенную матрицу с вещественными коэффициентами. К таким матрицам сводится множество физических задач. При этом матрицу коэффициентов будем считать (и сделаем) плотной. Вообще говоря, в реальной практике значительно чаще получаются не плотные, а сильно разряженные матрицы. Разработка эффективных алгоритмов для таких матриц – отдельное направление современной науки. Хорошо известны решатели PARDISO, MUMPS, SPOOLES, TAUCS, SuperLU, UMFPACK, PSPASES, WSMP, PaStiX и другие. При работе с разряженными матрицами на любом языке обычно вызывается одна из этих программ. Их настройка «в реальной жизни» и тестирование – очень важная и интересная для нас тема, но сейчас пока оставим ее за скобками. Разряженные матрицы иногда можно свести к ленточным, что чрезвычайно эффективно при малой ширине ленты. Ленточные матрицы тоже пока не рассматриваем.

Те средства, где доступны специальные методы решения СЛАУ с симметрической положительно определенной матрицей (на основе \(LL^T\) или \(U^TU\)-разложения Холецкого), протестируем два раза: с решением обычным методом для матриц общего вида (на основе классического \(LU\)-разложения) и методом с \(LL^T\)-разложением.

Будем тестировать следующие программные средства:

  • язык Fortran с различными математическими библиотеками (реализациями LAPACK/BLAS и другими);

  • язык Julia «из коробки», в котором доступны как «стандартные средства» решения СЛАУ (скрывающие используемый метод), так и вызовы конкретных процедур библиотеки LAPACK;

  • язык Python с математическими библиотеками NumPy и SciPy из дистрибутива Anaconda;

  • MatLAB-подобные открытые системы компьютерной математики (СКМ), установленные «из коробки»: Octave, SсiLAB и FreeLAB;

  • некоторые «прочие» открытые СКМ. Тут пока остановимся на СКА Maxima.

Все тесты проводились на обычном современном персональном компьютере с 4-ядерным 8-поточным процессором Intel i7 с достаточным размером памяти. Конкретную марку, частоты и т.п., думаю, описывать нет нужды, т.к. нас интересует лишь относительное сравнение средств.

ОС: Windows 10, 64. Вест софт (где это было доступно), использовался в 64-разрядной реализации.

Тестовая задача

Решим \(m\) раз систему из \(n\) уравнений

\[ \mathbf{A}\cdot\mathbf{X}=\mathbf{B}, \]

где \(\mathbf{A}\) – матрица коэффициентов размером \(n \times n\), \(\mathbf{X}\) – вектор неизвестных и \(\mathbf{B}\) – вектор правых частей.

Пусть коэффициенты матрицы \(\mathbf{A}\) вычисляются следующим образом:

\[ a_{i,j} =
\begin{cases}
1- \left(\displaystyle\frac{i-j}{n}\right)^2, &\text{если } i \ne j,\\
10, &\text{если } i = j.
\end{cases}
\]

При этом получается симметрическая положительно определенная матрица.

Например, для \(n=5\) матрица \(\mathbf{A}\) будет такой:

  10.00   0.96   0.84   0.64   0.36
   0.96  10.00   0.96   0.84   0.64
   0.84   0.96  10.00   0.96   0.84
   0.64   0.84   0.96  10.00   0.96
   0.36   0.64   0.84   0.96  10.00

Пусть на каждом шаге решения \( k = 1 \dots m \) коэффициенты вектора правых частей равны номеру шага: \( b_{i} = k \).

Чтобы контролировать правильность решения каждым программным средством будем подсчитывать сумму коэффициентов вектора \(\mathbf{X}\) на каждом шаге решения и суммировать ее по шагам:

\[ S = \sum_{k=1}^{m} \sum_{i=1}^{n} x_i^{(k)}. \]

Все программы будем тестировать при размере матрицы \(n=1000\) и числе повторений \(m=100\) шагов. При этих параметрах должно получаться \( S \approx 22 439,195\). Во всех программах будем использовать векторы и матрицы из вещественных чисел двойной точности (занимающие в памяти 8 байт).

Договоримся, что разложение (факторизацию) матрицы \(\mathbf{A}\) будем делать на каждом шаге, т.е. каждый раз систему уравнений будем решать полностью.

Время решения будем измерять по циклу решения СЛАУ – без учета предварительной подготовки матрицы \(\mathbf{A}\).

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

Решение различными средствами

Fortran

GFortran + BLAS/LAPACK (LU)

Для начала опробуем решение задачи на GNU GFortran (версия 5.3) с классическими библиотеками BLAS и LAPACK (версия 3.7.0). Установка описана в прошлом посте.

Программа для решения тестовой задачи выглядит так:

program LinTest
    implicit none
    integer, parameter      :: n = 1000 ! Число уравнений
    integer, parameter      :: m = 100  ! Число шагов
    real(8), dimension(n,n) :: A, A0    ! Матрица коэффициентов
    real(8), dimension(n)   :: B        ! Вектор правых частей / результатов
    real(8)                 :: s        ! Контрольная сумма
    integer                 :: i, j, k  ! Вспомог. индексы
    integer, dimension(n)   :: ipiv     ! Вспомог. матрица для процед. LAPACK
    integer                 :: info     ! Инфо об успешности решения
    integer                 :: t_0, t_1, t_rate, t_max ! Для измер. времени
    !------------------------------------------------------------------------
    forall(i=1:n, j=1:n) A0(i,j) = 1-(dble(i-j)/n)**2  ! Формируем матрицу
    forall(i=1:n, j=1:n, i==j) A0(i,j) = 10.
    s = 0.
    call system_clock(t_0, t_rate, t_max) 
    do k = 1, m
        A = A0
        B = k   ! Формируем вектор правых частей
        call dgesv(n, 1, A, n, ipiv, B, n, info) ! Решаем СЛАУ 
        if (info/=0) stop ("Ошибка")
        s = s + sum(B)
    enddo
    call system_clock(t_1, t_rate, t_max)
    print *, "Время:", real(t_1-t_0)/t_rate
    print *, "Сумма:", s
endprogram

Здесь использована процедура LAPACK dgesv, реализующая решение СЛАУ с \(LU\)-разложением и подходящая для матриц общего вида.

Матрицу A приходится восстанавливать из эталонной A0 на каждом шаге цикла т.к. dgesv ее изменяет (как и вектор B, куда записывается решение).

Компиляция:

gfortran LinTest.f90 -static -llapack -lrefblas -O2 -o LinTest.exe

Время решения задачи в таком варианте составило 14,6 с. При этом задействовался только один процессор. Библиотеки LAPACK/BLAS в классическом варианте – непараллельные. Видно, что в данном варианте они не оптимизированы – время расчета вызывает разочарование.

Размер exe-файла получился 821 КБ.

Компиляция с более агрессивными ключами оптимизации ускорения почти не дала.

Если вынести факторизацию матрицы

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

Для решения СЛАУ с матрицей общего вида в два приема в LAPACK имеются функции dgetrf и dgetrs. Первая реализует \(LU\)-разложение, вторая на его основе решает СЛАУ. Фрагмент кода в зоне главного цикла станет таким:

call dgetrf(n, n, A0, n, ipiv, info)  ! LU-разложение
do k = 1, m
    B = k
    call dgetrs("N", n, 1, A0, n, ipiv, B, n, info)  ! Решение СЛАУ
    if (info/=0) stop ("Ошибка") 
    s = s + sum(B)
enddo

Время решения задачи стало всего лишь 0,21 с, т.е. в 70 раз меньше (!). Оно и понятно. Трудоемкая часть в решении СЛАУ – факторизация.

Результат решения, разумеется, не изменился т.к. матрица \(\mathbf{A}\) на разных шагах решения по условию задачи у нас одинаковая.

GFortran + BLAS/LAPACK (LLt)

Для решения СЛАУ с симметрическими положительно определенными матрицами в LAPACK имеются процедуры на основе \(LL^T\)-разложения Холецкого. В исходной программе достаточно заменить одну строчку:

call dposv("L", n, 1, A, n, B, n, info) ! Решаем СЛАУ 

Время решения составило 8,6 с.

Размер exe-файла – примерно тот же.

Если вынести факторизацию матрицы

Фрагмент кода в зоне главного цикла:

call dpotrf("L", n, A0, n, info)  ! LU-разложение 
do k = 1, m 
    B = k
    call dpotrs("L", n, 1, A0, n, B, n, info)  ! Решение СЛАУ
    if (info/=0) stop ("Ошибка") 
    s = s + sum(B)
enddo

Время решения стало 0,16 с – в 55 раз меньше.

GFortran + OpenBLAS (LU)

OpenBLAS (www.openblas.net/) – свободная библиотека с функциями BLAS+LAPACK, оптимизированная для современных процессоров.

Устанавливается в GFortran OpenBLAS очень просто. Достаточно скачать версию со скомпилированными бинарниками и положить библиотеку libopenblas.a туда, где GCC хранит библиотеки .a. Текущая версия 0.2.19.

Код программы идентичен предыдущим примерам.

Компиляция:

gfortran LinTest.f90 -static -lopenblas -O2 -o LinTest.exe

Время решения составило 0,79 с. При этом задействовались все процессоры.

Слинкованный статически exe-файл получился очень большим: 26 (!) мегабайт. Можно ли его уменьшить я не разобрался. В дистрибутиве OpenBLAS прилагается также dll-библиотека. Но динамическая линковка с ним у меня не заработала.

Вынос факторизации за тело цикла уменьшил время до 0,03 с (!).

GFortran + OpenBLAS (LLt)

Время решения составило 0,51 с.

Это – шикарно!

Если бы еще размер exe-файла получался поменьше…

Вынос факторизации за тело цикла уменьшил время до 0,07 с.

GFortran + LAIPE2 (LU)

Библиотека Laipe2 (www.equation.com) и способ ее установки были упомянуты в предыдущем посте.

Это оригинальная реализация параллельных решателей СЛАУ с отличной от LAPACK идеологией. Для матриц общего вида вызов совмещенного решателя осуществляется так:

call laipe$Solution_DAG_8(A, n, B, info) ! Вызов процедуры LAIPE2

(Доступны также процедуры с разложениями и последующим решением).

Предварительно нужно инициализировать используемое число ядер (потоков) процессора (в моей случае – максимум 8). Без этого работать не будет:

call laipe$use(8)

Компиляция:

gfortran LinTest.f90 -static -fdollar-ok -llaipe2 -lneuloop4 -O2 -o LinTest.exe

Время решения – 7,4 с. Не густо. При этом все ядра загружаются максимально.

Размер exe-файла – 910 КБ. Это – приятно.

GFortran + LAIPE2 (LLt)

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

Приведу полный код – он отличается от прошлых реализаций:

program LinTest
    implicit none
    integer, parameter      :: n = 1000 ! Число уравнений
    integer, parameter      :: m = 100  ! Число шагов
    real(8), dimension(((n+1)*n)/2) :: A, A0 ! Профиль матрицы коэффициентов
    real(8), dimension(n)   :: B        ! Вектор правых частей / результатов
    integer, dimension(n)   :: Label    ! Массив меток
    real(8)                 :: s        ! Контрольная сумма      
    integer                 :: i, j, k  ! Вспомог. индексы
    integer                 :: info     ! Инфо об успешности решения
    integer                 :: t_0, t_1, t_rate, t_max ! Для измер. времени
    !------------------------------------------------------------------------
    Label(1) = 1      ! Формируем массив меток
    do j = 2, n
        Label(j) = Label(j-1) + n-j+1
    enddo
    do j = 1, n       ! Формируем профиль матрицы коэффициентов
        do i = j, n
            A0(Label(j)+i-1) = 1-(dble(i-j)/n)**2  ! j-столбец, i-строка
            if (i==j) A0(Label(j)+i-1) = 10.
        enddo
    enddo
    s = 0.
    call laipe$use(8) ! Число используемых потоков
    call system_clock(t_0, t_rate, t_max)
    do k = 1, m       ! Решение задачи
        A = A0
        B = k
        call laipe$Solution_DSP_8(A, n, Label, B, info)  ! Вызов LAIPE2
        if (info==1) stop ("Ошибка") 
        s = s + sum(B)
    enddo
    call system_clock(t_1, t_rate, t_max)
    print *, "Время:", real(t_1-t_0)/t_rate
    print *, "Сумма:", s
endprogram

Время решения в данной реализации составило 4,9 с. Ядра загружаются максимально.

Размер exe-файла – 910 КБ.

Intel Fortran + MKL (LU)

В интернете оказался доступен триал коммерческого компилятора Intel Fortran с математической библиотекой MKL, содержащей все те же BLAS и LAPACK, но хорошо оптимизированные для современных процессоров. Он стоит немалых денег. Хотя настоящий блог посвящен свободному ПО, я не мог не попробовать. Очень интересно сравнить.

Код программы остается все тем же, что и для классических BLAS/LAPACK.

Компиляция:

ifort /Qmkl /static /O2 LinTest.f90

Время решения составило 0,58 с. Конечно, великолепно.

Размер exe-файла – 6,3 МБ. Тоже немало, но не 26 все же. :)

Intel Fortran + MKL (LLt)

Время решения составило всего 0,33 с. Круто, что скажешь. Забегая вперед, это абсолютный рекорд из всех протестированных решений. Программисты Intel едят хлеб не зря. Intel Fortran, определенно, своих денег стоит.

Julia

У меня установлена Julia 0.5.0. У нее «под капотом» – реализация LAPACK OpenBlas, очевидно, аналогичная той, которую мы уже «щупали» под Fortran (может быть версия постарее).

Julia, стандартное средство (LU)

Код решения задачи на Джулии – лаконичен и понятен без комментариев:

n = 1000; m = 100
s = 0.
A = [1-((i-j)/n)^2 for i in 1:n, j in 1:n] + 9*eye(n,n)
  @time for k=1:m
    B = fill(Float64(k), n)
    X = A \ B
    s += sum(X)
  end
println(s)

По сравнению с кодом на Fortran, что называется, «почувствуйте разницу».

Время решения в таком варианте составило 1,2 с. Задействовались все ядра.

Макрос @time выводит как время выполнения цикла, а также объем выделенной памяти. В данном случае 765 МБ.

Функция «\» скрывает фактический метод решения. На самом деле он выбирается в зависимости от типа матрицы А.

Например, можно явно указать, что матрица симметричная:

A = Symmetric(A)

После этого время решения почему-то несколько возрастает… :)

Скорее всего в исходном варианте использовалась LU-факторизация. Ее можно выполнить предварительно явно. Заменим строку с решением уравнения:

X = lufact(A) \ B

При этом функция «\» так же выберет необходимый метод с учетом типа аргументов, поступающих на вход. Первый аргумент будет иметь тип уже не матрицы, а Base.LinAlg.LU{Float64,Array{Float64,2}}, содержащий в себе факторизованную матрицу. Функция «\» это поймет.

Время решения в таком варианте не изменилось.

Если вынести факторизацию матрицы

Если lufact(A) вынести за тело цикла программа будет выглядеть так:

n = 1000; m = 100
s = 0.
A = [1-((i-j)/n)^2 for i in 1:n, j in 1:n] + 9*eye(n,n)
tic()
  luA = lufact(A)
  for k=1:m
    B = fill(Float64(k), n)
    X = luA \ B
    s += sum(X)
  end
toc()
println(s)

(Тут изменен также способ измерения времени на tic/toc).

Время решения уменьшилось до 0,04 с – в 3 раза.

Julia, стандартное средство (LLt)

Чтобы решить задачу с предварительным разложением Холецкого в исходном коде меняем строку с решением СЛАУ:

X = cholfact(A) \ B

Время получается 0,98 с.

Julia, LAPACK (LU)

Из Джулии также можно напрямую вызвать процедуры LAPACK. Вариант с LU-разложением и вызовом процедуры gesv! будет выглядеть так:

import Base.LinAlg.LAPACK: gesv!
n = 1000; m = 100
s = 0.
A0 = [1-((i-j)/n)^2 for i in 1:n, j in 1:n] + 9*eye(n,n)
  @time for i = k:m
    B = fill(Float64(k), n)
    A = copy(A0)
    gesv!(A, B)
    s += sum(B)
  end
println(s)

Здесь матрицу А приходится восстанавливать из эталонной A0, т.к. gesv! ее изменяет. Отметим особенности:

  • написать A = A0 нельзя, т.к. в этом случае будет присвоен только указатель на матрицу;

  • функции, меняющие значения аргументов, в Джулии принято называть с восклицательным знаком.

Время решения – 1,3 с.

Julia, LAPACK (LLt)

Вариант с LU-разложением и вызовом процедуры LAPACK posv!:

import Base.LinAlg.LAPACK: posv!
n = 1000; m = 100
s = 0.
A0 = [1-((i-j)/n)^2 for i in 1:n, j in 1:n] + 9*eye(n,n)
  @time for k = 1:m
    B = fill(Float64(k), n)
    A = copy(A0)
    posv!('L', A, B)
    s += sum(B)
  end
println(s)

Время решения – 0,78 с.

Python

У меня установлен Python 3.5.2 из дистрибутива Anaconda 4.2.0 (64-bit), в который «из коробки» входят библиотеки NumPy 1.11.1 и ScyPy 0.18.1.

«Под капотом» у NumPy/SciPy – Intel MKL (интересно, почему это остается бесплатным?). Ожидаем отличную производительность.

Python + NumPy (LU)

Код:

import numpy as np
import time
n = 1000; m = 100
s = 0.
A = np.fromfunction(lambda i,j: 1-((i-j)/n)**2, (n,n)) + 9*np.eye(n)
t = time.time()
for k in range(1, m+1):
    B = np.full((n), float(k))
    X = np.linalg.solve(A, B)
    s += np.sum(X)
print("Время:", time.time() - t)
print("Сумма:", s)

Время работы – 0,93 с. Отличный результат для матрицы общего вида.

Модуль линейной алгебры у NumPy скромный. Готовая функция для решения СЛАУ здесь только для матриц общего вида, вызывающая ту же функцию LAPACK gesv.

Библиотека SciPy – намного богаче. В частности, в нее входят функции решения СЛАУ с матрицами различного вида.

Python + SciPy (LU)

Изменения:

    ...
    from scipy import linalg
    ...
    X = linalg.solve(A, B)
    ...

Время – 1,2 с. Странно, почему, больше чем у NumPy? Видимо там немного разные версии MKL вызываются.

Python + SciPy (LLt)

Изменения:

    ...
    from scipy import linalg
    ...
    X = linalg.solve(A, B, sym_pos = True)
    ...

Время – 0,98 с.

MatLAB-подобные СКМ

Octave (LU)

У меня установлен Octave 4.2.1. Для решения СЛАУ он использует OpenBLAS.

Код:

n = 1000; m = 100;
tic
A = zeros(n,n);
for i = 1:n
    for j = 1:n
        A(i,j) = 1-((i-j)/n)**2;
        if (i == j)
            A(i,j) = 10.;
        endif
    endfor
endfor
toc
s = 0.;
tic
for k = 1:m
    B = k * ones(n, 1);
    X = A \ B;
    s = s + sum(X);
endfor
toc
s

Цикл с решением СЛАУ выполнился за достойное время 1,4 с (все-таки здесь работает OpenBLAS), а вот предшествующая часть с заполнением матрицы отняла неприличные 7 секунд (!). Я в ступоре и прострации – что я делаю не так??? Во всех остальных «приличных» языках эта часть занимает пренебрежимо малое время.

SciLab (LU)

У меня установлен SciLab 6.0.0. Для решения СЛАУ он использует Intel MKL.

Код:

n = 1000; m = 100
A = zeros(n,n)
for i = 1:n
    for j = 1:n
        A(i,j) = 1-((i-j)/n)**2
        if i == j then
            A(i,j) = 10.
        end
    end
end
s = 0.
tic()
for k = 1:m
    B = k * ones(n, 1)
    X = A \ B
    s = s + sum(X)
end
disp(toc(), s)

Время – 1,2 с.

FreeMat (LU)

FreeMat – MatLAB-подобная СКМ, свободная и относительно легковесная, за что многие коллеги ее очень любят. К сожалению, проект не развивается с 2013 года. Доступна только в 32-битной версии.

tic;
n = 1000; m = 100;
A=zeros(n,n);
for i = 1:n
    for j = 1:n 
        A(i,j) = 1-((i-j)/n)^2;
        if (i == j) 
            A(i,j) = 10.;
        end
    end
end
disp(toc)
s=0.;
tic;
for k = 1:m
    B = k * ones(n, 1);
    X = A \ B;
    s = s + sum(X);
end
disp(toc, s)

Время решения задачи составило 96 с.

Время заполнения матрицы А – 15 с.

Ясно, что FreeMat для больших вычислений не предназначен.

Maxima

Maxima – свободная, мощная и очень старая система компьютерной алгебры (СКА). Не смотря на преклонный возраст проект развивается, периодически появляются новые версии. Текущая – 5.39.0. Я привык работать с графической оболочной wxMaxima. Когда нужно что-то быстро посчитать, сделать символьные преобразования – отличная вещь. Грешен: старушку Максиму я по-прежнему нежно люблю.

Конечно, Maxima не предназначена для вычислений большого объема. Однако, система – универсальная. В частности, в ней есть модуль линейной алгебры с вызовом библиотек LAPACK (в какой реализации – непонятно, видимо несовременной).

Решим нашу задачу и в Maxima, хотя, это скорей «по приколу» (ну и чтобы не забывать ее несколько специфичный язык):

load(lapack)$

n: 1000$  m: 100$
t0: elapsed_run_time()$
A: zeromatrix (n, n)$
for i:1 thru n do
  for j:1 thru n do
    A[i,j]: if (i#j) then float(1-((i-j)/n)**2) 
                     else float(10)$
B: zeromatrix (n,1)$
s: 0.$
t1: elapsed_run_time()$
t1-t0;
for k:1 thru m do (
    B: B + 1.,
    X: dgesv (A, B),
    s: s + sum(X[i,1], i, 1, n)
)$
elapsed_run_time()-t1;
s;

Время заполнения матрицы А составило аж 70 с.

Время решения задачи… спустя несколько часов я не дождался. :) На малых объемах данных приведенное решение работает правильно, но при n=1000 и m=100 – все о-очень долго.

Сравнение и заключение

Мы решали систему уравнений с матрицей вещественных чисел двойной точности, размером 1000 х 1000, повторяя решение в цикле 100 раз.

Общие сравнительные результаты выглядят следующим образом:




Это иллюстрация очевидного факта: скорость решения СЛАУ разными средствами может колоссально разниться – в разы, на порядки…

Не предназначенная для решения подобных задач Maxima работала так долго, что я решения не дождался. Ей можно простить. Неоптимизированный FreeMat, позиционируемый как аналог MatLAB, тоже рекордно забуксовал. Разочаровал GFortran с классическими BLAS/LAPACK без распараллеливания и оптимизации.

Если откинуть перечисленных «аутсайдеров» получается следующая картина:




Абсолютный лидер – коммерческий Intel Fortran + MKL (LLt). Время решения задачи у него заняло всего 0,33 с. Инженеры Intel едят хлеб не зря. Если надо создать по-настоящему быструю реализацию на Fortran чего-либо серьезного, покупка лицензии Intel – правильное решение.

Догоняющий лидера свободный аналог – GFortran + OpenBLAS (LLt) с временем 0,51 с. Единственный недостаток этой реализации – огромный exe-файл (26 МБ). Это несколько напрягает. Надо будет попробовать с Фортраном другие оптимизированные реализации BLAS/LAPACK: например, GotoBLAS2 или ScaLAPACK.

Библиотека LAIPE2 от профессора Ляо как-то разочаровала.

В протестированных современных высокоуровневых средствах – языках Julia, Python и СКМ Octave и SciLAB в качестве решателей СЛАУ используются скомпилированные оптимизированные библиотеки Intel MKL или OpenBLAS. Поэтому результаты у них получились весьма достойными, а между собой – близкими. Выигрышней всего смотрится Julia + LAPACK (LLt) – 0,78 с. Единственное, неприятно удивил Octave, который собирал исходную матрицу неприличные 7 секунд. Может, я что-то сделал не так?

В целом же, если львиную долю времени работы какой-либо инженерной программы отнимает именно решение СЛАУ, а за рекордами скорости нет нужды гнаться, можно смело писать ее на одной из данных высокоуровневых систем. По сравнению с «низкоуровневым» Фортраном это сэкономит массу времени в разработке, а также нервные клетки.

Очевидная истина – решатель СЛАУ надо выбирать соответственно типа матрицы. В нашем примере была симметрическая положительно определенная матрица. Все решения с \(LL^T\)-разложением для нее были заметно быстрее аналогов с \(LU\)-разложением, подходящим для матриц общего вида. Для разряженных матриц следует использовать соответствующие решатели. Если можно сделать разряженную матрицу ленточной с малой шириной ленты – будет совсем хорошо. Аналогичный тест для разряженных матриц – надеюсь, еще проведем. Там свой «мир» решателей, изучить и освоить его – было бы очень полезно и интересно.

Если реализуется некоторый, скажем итерационный, алгоритм, в котором СЛАУ решается многократно, а при этом матрица коэффициентов не изменяется, обязательно надо выносить ее разложение за тело цикла. Это даст огромный выигрыш по времени – до десятков или даже сот раз.

В заключение еще раз отметим, что проведенный тест – «бытовой» и мягко говоря не объективный. Причины были перечислены во введении. Пожалуй, для меня самый интересный его результат – тренировка в работе с различными средствами.

четверг, 26 января 2017 г.

Механические расчеты конструкций в свободном ПО. Часть 1

Введение

Классическая «технология» расчета конструкций МКЭ включает 4 этапа:

  1. Твердотельное моделирование с применением 3D-CAD’ов.

  2. Предпроцессинг: разбиение на конечные элементы.

  3. Расчет.

  4. Постпроцессинг: визуализация и анализ результатов.

В мире существует огромная «армия» коммерческого ПО класса CAD/CAE, решающая эти задачи как вместе – в одной системе-комбайне, так и по отдельности.

В нашей отрасли наиболее популярны для 3D-моделирования SolidWorks, КОМПАС 3D и Autodesk Inventor, а для пунктов 2–4 – либо соответствующие модули в составе этих систем (SolidWorks Simulation, APM-FEM в КОМПАС-3D или Stress Analysis в Inventor Professional), либо – более серьезный ANSYS.

Расчеты – дело ответственное. Для сложных систем лучше использовать профессиональные CAD/CAE-системы.

Однако, например, в нашей практике в 90% случаев приходится иметь дело либо с «очень», либо с «довольно» простыми конструкциями. Для их расчетов отлично зарекомендовал себя рассмотренный ниже свободный софт, справляющиеся с возложенными на него функциями не хуже коммерческих монстров. Применение свободных программ оправдано, когда подобные расчеты приходится делать «время от времени», а сами конструкции (повторюсь) – не слишком сложны. Завязываться на свободный софт по его текущему уровню развития в технологическом цикле машиностроительного проектирования CAD/CAE/CAМ по моему мнению пока преждевременно.

Мы многократно тестировали описанную ниже цепочку программ на предмет адекватности результатов:

  • сравнением с модельными задачами, для которых можно получить аналитическое решение (балки, стержневые системы);

  • сравнением с результатами, полученными в «серьезных» CAD/CAE;

  • испытаниями.

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

Пожалуй, единственный значимый недостаток – это определенное неудобство использования. «Порог входа» довольно высок. На разных этапах приходится применять программы с разными идеологиями и разными не вполне интуитивными интерфейсами. Как в любом свободном ПО, у этих программ немало недоработок. Но, если к их особенностям чуть-чуть попривыкнуть, «рука набивается». Расчеты получается делать достаточно быстро и качественно.

Важное условие – инженер должен понимать, что он делает и иметь определенный багаж теоретических знаний в строительной механике / сопромате и МКЭ. Но без этого – легко «натворить дел» и в коммерческих CAE. Их обманчивая простота «для домохозяек» продуцирует ситуацию, когда за расчеты берутся мало понимающие в них люди. При этом мотивации к изучению теории – не возникает. В итоге мы имеем настоящую эпидемию – море неграмотно сконструированных объектов. Иногда все заканчивается очень печально. В этом смысле, свободный софт, наоборот, мотивирует знать теорию. Без нее, скорее всего, просто ничего не получится. :)

Отметим, что вариантов использования различных программ на разных этапах «технологии» – множество. Перечислю в таблице лишь некоторые:

Этап Возможный свободный (бесплатный) софт
1. 3D-моделирование OpenSCAD, ImplicitCad, FreeCAD, CadQuery, SolveSpace, OnShape, Salome, Calculix cgx, GMSH, BRL-CAD
2. Пред-процессинг GMSH, Salome, NetGen, TetGen
3. Расчет CalculiX cсx, Code_Aster, GetDP, NGSolve
4. Пост-процессинг CalculiX cgx, Salome, GMSH, NetGen, ParaView

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

Чтобы было не скучно, последовательность расчета будет проиллюстрирована не на классических балках или «пистонах». Поиздеваемся над символом этого блога – звездчатым многогранником «Соединение пяти тетраэдров». Поставим скромную задачу выполнить самый простой линейный расчет напряженно-деформированного состояния. Более интересные нелинейные расчеты, расчеты устойчивости и динамический анализ – также без проблем делаются. Надеюсь о них написать как-нибудь в будущем.

CAD’ы: подготовка 3D-геометрии

OpenSCAD: 3D-моделер для программистов

OpenSCAD – это легкая и простая, но совершенно крутейшая штука. 3D-модели здесь просто «пишутся» на удобном и очень простом языке программирования. Рядом – интерактивно строится результат.


Теоретически на нем можно делать весьма сильные вещи. Умельцы и делают. Когда нужно сформировать что-то «модельное» или «абстрактное» – я тоже с удовольствием его запускаю. Коллеги-конструкторы крутят у виска пальцем. :)

Вообще, идея описания конструкций на языке – на мой взгляд прекрасна. Она не получила развития в коммерческих CAD’ах скорее в силу традиций (не приучены конструкторы к текстам). А ведь достоинств у такого подхода немало:

  1. Текст – прозрачен и гибок. Модель параметризуется и «работает» в точности так, как вы задумали.

  2. Универсальность и выразительность языка – всегда выше GUI.

  3. Текст хорошо «накапливается» и обобщается в библиотеки, которые дальше повторно используются. Уровни абстракций – наращиваются.

  4. Текст не может «испортиться», как непонятный бинарный файл, в котором хранят модели обычные CAD’ы.

  5. С текстом можно использовать системы контроля версий, такие как Git. Одно это очень многого стоит.

  6. Языковая модель CAD наиболее легко интегрируется с узкоспециальными САПР вашей предметной области. Это особенно актуально, если части такой САПР приходится разрабатывать самостоятельно.

  7. С учетом впечатляющего прогресса в развитии языков программирования можно было бы сделать очень продвинутый и лаконичный язык описания конструкций. Можно было бы совместить описание на языке и визуальное проектирование элементов, которые в тексте задавать неудобно (например, как в конструкторах форм IDE).

Не все пункты это перечня реализуются в OpenSCAD. :) Вообще, он – скорее концепт. В реальной работе, когда нужно быстро посчитать что-то не элементарно простое, в нем работать, конечно, не очень удобно. Приходится выбирать средства визуального проектирования, о которых речь пойдет ниже.

Но когда надо сделать что-то «модельное», «математически абстрактное», параметрическое, OpenSCAD – крут. Например, вот описание нашего многогранника:

module tetrahedron(edgeLength){ // Правильный тетраэдр 
    points =  [ [ 1, 1, 1],
                [ 1,-1,-1],
                [-1, 1,-1],
                [-1,-1, 1] ];
    faces  =  [ [ 0, 2, 1],  
                [ 3, 0, 1],  
                [ 0, 3, 2],  
                [ 1, 2, 3] ];    
    scale([1,1,1]*edgeLength/sqrt(8))    
    polyhedron(points, faces);    
}

module fiveTetrahedra(edgeLength){ // Соединение пяти тетраэдров 
    for(i=[0:4]){
        rotate([0,-atan(sqrt(5)/2-0.5), i*360/5])
        tetrahedron(edgeLength);
    } 
}

// ---------------- Модель -------------------
rotate($t*360/5)
fiveTetrahedra(100);    

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

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

Предпоследняя строчка нашего текста, где выполняется преобразование rotate с параметром $t, позволяет сделать в окне OpenSCAD анимацию. Нельзя удержаться, чтобы не перегнать ее в .gif: :)


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

OpenSCAD очень популярен среди любителей 3D-печати. Кстати, они наплодили буквально тонны готовых текстов для разных моделей, открытых для повторного применения. OpenSCAD позволяет сохранять модели в сеточном формате STL, стандартном для 3D-печати.

OpenSCAD – самый известный, но далеко не единственный проект «языкового» CAD. Автор OpenSCAD Мариус Кинтель намерено сохраняет язык простым. Это по-своему здорово, но возможностей языка хватает далеко не всегда. В частности, огорчает отсутствие многих средств, к которым мы привыкли в языках программирования общего пользования. Имеется несколько альтернативных проектов. Основные из них представляют похожие средства описания геометрии вместе с продвинутыми возможностями различных базовых языков программирования.

Альтернативный проект (ссылка) Базовый язык Особенности
1 OpenJSCAD JavaScript Веб-интерфейс
2 CoffeeCAD CoffeeScript Веб-интерфейс
3 RapCAD OpenSCAD с дополнениями Заточен для RepRap-принтеров. GUI, похожий на OpenSCAD.
4 ImplicitCAD OpenSCAD / Haskell Написан целиком на Haskell. Добавляет к языку OpenSCAD много интересных возможностей. В частности, можно легко выполнить сглаживание граней и плавные переходы между частями моделей. Веб-интерфейс + работа из командной строки.
5 OpenSCAD.net OpenSCAD / OpenJSCAD Веб-интерфейс для OpenSCAD / OpenJSCAD
6 CadQuery Python Надстройка к FreeCAD либо веб-интерфейс. Позволяет работать с BREP-геометрией. Предоставляет высокоуровневые абстракции, позволяющие определять модель подобно тому, как описал бы ее человек на родном языке

Я с особым интересом слежу на проектом ImplicitCad (пункт 4), т.к. Haskell – это безумно круто. :) И вообще: проект создает интересный канадский чувак Кристофер Олах, мыслящий очень продвинуто. Но особый интерес для нас представляет пункт 6 – CadCuery (я планирую рассказать о нем немного подробней). Помимо того, что он открывает нам мощь Python, он интересен вот почему.

Огромный недостаток OpenSCAD (а также альтернатив пп. 1-5) для нас заключается в том, что он «мыслит» полигональными сетками, описывающими кривые поверхности приближенно. OpenSCAD базируется на геометрическом ядре CGAL и реализует технологию конструкционной блочной геометрии (CSG) с помощью сеток. По этой причине у него нет возможности сохранять модель в «точном» формате, описывающем геометрию на основе граничных поверхностей (BREP) в форматах .STEP или .BREP. А для последующих расчетов нам крайне желательно иметь геометрию в «точном» формате. Точно также OpenSCAD не может и импортировать геометрию, скажем, из .STEP. Очень жаль: иначе можно было бы использовать части модели, сделанные где-то еще, видоизменяя их на языковом CAD. Всех этих недостатков лишен CadQuery, и о нем я еще напишу.

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

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

// ---------------- Модель -------------------
rotate($t*360/5)
union() {
    fiveTetrahedra(100);    
    sphere(64*sqrt(6)/4, $fn=100);
}

Здесь радиус шара задан как \( 64\sqrt{6}/{4} \) т.к. для правильного тетраэдра радиус описанной сферы равен \( a \sqrt{6}/{4}\), где \(a\) – длина ребра. Можно проверить, что при \( 100\sqrt{6}/{4} \) сфера опишет наши тетраэдры.

Получаем:


Обратите внимание, что сферическая поверхность отрисована сеткой. $fn=100 – это ее «густота».


(Продолжение следует).

вторник, 24 января 2017 г.

GFortran: доступ к Windows API

Сам вызов функций WinAPI делается на Cи (на примере MessageBox):

// Файл: winapi.c
// *** Вызовы функций Windows API ***
#include <windows.h>

int c_messagebox_(HWND hWnd, LPCTSTR lpText, LPCTSTR lpCaption, UINT uType)
{ 
    return MessageBox(hWnd, lpText, lpCaption, uType); 
}

Названия функций, которые мы затем будем вызывать из Фортран (в примере – c_messagebox_), должны заканчиваться знаком подчеркивания. В Фортране его не будет – это связано с принципами линковки в GCC.

Для доступа к этим функциям из Фортран-программ целесообразно сделать модуль с обертками:

! Файл: windows.f08
! *** Модуль с интерфейсом к некоторым функциям Windows API ***
module Windows
    use iso_c_binding

    implicit none
    integer(kind=C_LONG), parameter, private :: HWND_DESKTOP       = '00000000'X
    integer(kind=C_LONG), parameter, private :: MB_OK              = '00000000'X
    integer(kind=C_LONG), parameter, private :: MB_ICONERROR       = '00000010'X
    integer(kind=C_LONG), parameter, private :: MB_ICONINFORMATION = '00000040'X

    interface
        function c_messagebox(hWnd, lpText, lpCaption, uType)
            use iso_c_binding
            integer(kind=C_LONG)          :: c_messagebox
            integer(kind=C_LONG), value   :: hWnd, uType
            character(kind=C_CHAR, len=*) :: lpText, lpCaption
        endfunction 
    endinterface

contains 

    subroutine MsgErr(str)   ! Выводит простой MessageBox с ошибкой
        character(*)    :: str
        integer         :: res
        res = c_messagebox(HWND_DESKTOP, str // C_NULL_CHAR, &
              "Ошибка" // C_NULL_CHAR, MB_ICONERROR + MB_OK)
    endsubroutine 

    subroutine MsgInf(str)   ! Выводит простой MessageBox c информацией
        character(*)    :: str
        integer         :: res
        res = c_messagebox(HWND_DESKTOP, str // C_NULL_CHAR, &
              "Информация" // C_NULL_CHAR, MB_ICONINFORMATION + MB_OK)
    endsubroutine

endmodule

Здесь вызов функций из модуля на Си осуществляется через interface с указанием специальных разновидностей типов (C_LONG, C_CHAR), определенных для совместимости с Си в модуле iso_c_binding, входящим в GCC.

NB! В процедурах-обертках к строкам приходится добавлять C_NULL_CHAR.

Для проверки:

! Файл: test.f08
program Test
    use Windows
    call MsgInf("Это информация.")
    call MsgErr("А это ошибка!")
endprogram

Компилируем и запускам:

gcc -c winapi.c
gfortran -c windows.f08
gfortran -c test.f08
gfortran *.o -o test.exe
test.exe

Если все хорошо, получаем:


понедельник, 23 января 2017 г.

Ламповый теплый Fortran

Introduction

Инженер на любом языке программирования программирует на Fortran.

Это шутка, но, как обычно, только отчасти. :)

Fortran – старейший язык высокого уровня, задуманный изначально как FORmula TRANslator для ученых и инженеров. Fortran создавался в IBM начиная с 1954 года суровым чуваком Джоном Бэкусом.


Формально первым языком высокого уровня считается Планкалкюль (1945), компилятор которого, правда, был создан только в 2000. Fortran в 1950-х был не просто изобретен (точнее – «интуитивно нащупан»), но сразу реализован. Он – первый «по-честноку».

Благодарные ученые и инженеры активно используют Fortran до сих пор, несмотря на очевидную его моральную устарелость. (Сейчас XXI век, четвертая пятилетка). Мы тоже не исключение. На Fortran меня подсадили в 2009 году седовласые профессора кафедры «Прикладная математика» петербургского Политеха, с которыми мы вместе создавали одну нетривиальную загогулину. С тех пор слезть с него не удается никак.

В рейтинге популярности языков программирования Fortran находится на 28 месте, между F# и Lua. Современные мейнстримные программисты заливаются хохотом, когда им говоришь про Fortran. (Боюсь, они представляют его весьма условно, в стандарте приблизительно 77). А зря. Среди куда более узкой аудитории – ученых и инженеров Fortran занимает если уже и не первое, то точно одно из ведущих мест.

Причин продолжения использования Fortran несколько:

  1. Огромная, накопленная за пол века база инженерно-научного кода.

  2. Скорость. Программы на Fortran – очень быстрые. Компиляторы Fortran чрезвычайно эффективны, а используемые математические библиотеки алгоритмически отточены как лезвия бритв. Кроме того, современный Fortran очень эффективно позволяет распараллеливать вычисления.

  3. Удобство работы с матрицами и вообще с математикой. В Си-подобных языках этого можно добиться только с помощью ряда библиотек, здесь – все «из коробки».

  4. Он простой как топор и императивный как танк. Инженеру/ученому не надо забивать голову программистскими премудростями (объектно-ориентированными или функциональными парадигмами, либо особенностями распределения памяти, указателями и т.д., и т.п.). Можно просто писать свои сладкие формулы.

  5. Он компилируемый. Программы очень просто кому-нибудь отдавать безо всяких ваших виртуальных машин.

  6. Как ни странно, он развивается и поддерживается. Стандарты регулярно обновляются (последний – 2015), компиляторы – совершенствуются.

На самом деле, Fortran действительно хорош для задач вычислительной математики, требующих очень высокой скорости исполнения кода с одной стороны, а с другой стороны – не слишком сложных в плане организации данных. Большинство околонаучных задач таковы. Однако, если данных все-таки много (не по объему, а по запутанности), а проект – сложен, то Fortran с удовольствием поспособствует его превращению в летающего макаронного монстра.

Если высокая скорость не требуется, лучше взять более современный язык «сверхвысокого» уровня: Python с библиотеками NumPy и SciPy или специализированный математический язык Mathematica или MatLAB (или свободные аналоги; для последнего – Octave или SciLAB). На них разработка «научного» проекта будет в разы быстрей и прозрачней, правда все будет работать чертовки по-черепашьи.

Часто проекты так и приходится создавать в два этапа:

  1. Прототипирование на «высоком и прозрачном» языке типа Python.

  2. Реализация на быстром Fortran (или С/С++).

В последние годы научным сообществом активно разрабатывается новый современный язык Julia, претендующий на замену Fortran. Он позволяет писать как «прозрачные» прототипы, так и быструю реализацию критичных участков. Julia, родившаяся в 2012, – многообещающий и очень интересный язык, который лично я изучать пока только начал. О нем еще будут посты.

Компиляторы Fortran

Живых компиляторов Fortran сегодня, как ни странно, не меньше десятка.

Лучший активно развивающийся коммерческий компилятор для Windows и Linux – вероятно, Intel Fortran. С ним «из коробки» идет мощнейшая математическая библиотека MKL, включающая известные пакеты линейной алгебры BLAS и LAPACK, а также супер-эффективный пакет для работы с разряженными матрицами PARDISO и прочие вкусности. (За дополнительные деньги можно приобрести дополнительные или альтернативные библиотеки, такие как знаменитая IMSL). С Intel Fortran поставляется очень мощный профайлер, позволяющий смотреть, на каких участках кода как тратится время (и другие параметры) и оптимизировать код. В Windows Intel Fortran ставится на MS Visual Studio и работает с линкером от Майкрософт.

Вторым лидером рынка вроде бы считается PGI Fortran (для Linux and Mac доступен во Free edition) – тоже монстр, на который, я правда особо и не смотрел. Еще есть Absoft Pro Fortran, IBM XL Fortran (для Linux и AIX), EKOPath (Linux), Lahey/Fujitsu Fortran, NAG Fortran, Oracle Solaris Studio (для Solaris и Linux, кажется недавно стала бесплатной) и другие.

Довольно интересный коммерческий компилятор – Silverfrost FTN95, бесплатный для персонального использования или «в целях оценки». Единственное неудобство – каждый раз при запуске скомпилированной программы появляется баннер. С FTN95 поставляется приятная легкая Plato IDE. Вообще, вся система производит впечатление очень легкой и быстрой.

Лучший свободный компилятор – GFortran из набора компиляторов GCC, разрабатываемых в рамках проекта GNU. GFortran активно совершенствуется «свободным сообществом».

Еще неплохой свободный компилятор был G95, но похоже загнулся в 2013.

Итого, наш выбор – GFortran, поэтому о нем – поподробнее.

Установка GFortran в Windows

На январь 2017 последняя стабильная версия GCC с GFortran – 6.3, заканчивается работа над седьмой версией.

GFortran – разумеется, кроссплатформенный. Для Windows его можно получить и установить разными способами:

  1. Из родного «набора свободных средств разработки под Windows для настоящих минималистов» :) – проекта MinGW или отпочковавшегося от него проекта MinGW-w64, реализующего преимущества 64-битной архитектуры. Или из проекта CygWin. Все три варианта – не самые простые пути. Что надо скачать и как это настроить – все довольно запутанно.

  2. Поставить GFortran вместе с хорошей свободной кроссплатформенной IDE Code::Blocs. Скачиваем с их сайта и устанавливаем файл codeblocks-16.01mingw_fortran-setup.exe (версия 16.01 на сегодня последняя). Возможно, это наиболее простой путь. Недостаток – GFortran будет несколько устаревшей (текущая версия в сборке 4.9.2) и только в разрядности 32. Зато поставил – и можно сразу переходить к сладким формулам в интуитивно понятном окружении IDE.

  3. От ребят, именующих себя Equation Solution, в лице китайца Джен-Чинга Ляо (www.equation.com). Это бывший профессор колумбийского университета, еженедельно собирающий самую актуальную версию дистрибутива GCC для разрядности как 32, так и 64. Сборка дополнена отладчиком и собственными библиотеками профессора Ляо, из которых особый интерес для нас представляют LAIPE2 (решение СЛАУ с распараллеливанием) и JCL (перенумерация строк и столбцов разряженных матриц для приведения их к ленточной форме). Это хорошие библиотеки, оказавшиеся для меня очень кстати. Из дистрибутива профессора Ляо GCC устанавливается очень легко через инсталлятор. Необходимые пути автоматом прописываются в PATH. У Ляо в принципе, все хорошо. С его дистрибутивами я жил какой-то время, единственное, там немножко все «собрано на коленке». Например, в комплекте нет dll (предполагается только статическая линковка библиотек gcc), получаемые exe в последних версиях стали странным образом разбухать, его собственные библиотеки – иногда глючат и т.п. :) Как обычно, Open Source – такой Open Source. :)

  4. Из специальной сборки компиляторов GCC для Windows с сайта TDM-GCC. Это очень хорошие и продуманные сборки, включающие полный набор нужных средств. К сожалению, здесь GCC тоже несколько устаревший (текущая версия 5.1), но сейчас этот вариант – мой фаворит. Скачиваем с их странички файл tdm64-gcc-5.1.0-2.exe для архитектуры 64 или tdm-gcc-5.1.0-3.exe для 32. На версии 64 можно (в том числе) компилировать свои программы и для 32-разрядных машин, т.е. две версии ставить не нужно. Устанавливается все элементарно через инсталлятор. Единственное, ему надо указать «Устанавливать все пакеты», т.к. по умолчанию компилятор фортран как «экзотика» пропускается. Инсталлятор прописывает путь к компиляторам в PATH.

Последний вариант, TDM-GCC, – очень хороший, за одним «но». В версиях 5.3 (и 4.9) допущен досадный баг, из-за которого не работает оператор фортрана Open. :) (Open Source – такой Open Source). Об этом не написал в баг-листе только ленивый. Надо скачать с их сайта старую версию 4.8 (там надо найти gcc-4.8.1-tdm64-2-fortran.zip) и заменить из нее файлы libgfortran_64-3.dll и libgfortran-3.dll, а также libgfortran.a, libgfortran.dll.a, libgfortran.spec, libgfortranbegin.a. Причем последние файлы имеются в 64 и 32-разрядных версиях, их надо аккуратно разложить по соответствующим папкам lib или lib32.

То же относится и к текущей версии 4.9 GFortran в CodeBlocs (он по сути там тот же самый из TDM-GCC).

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

В дистрибутивах от Ляо этой проблемы нет.

Компиляция в GFortran

Из IDE

IDE для того и нужны, чтобы не думать, а нажимать кнопочки. :) Все интуитивно понятно без дополнительных пояснений.

Из командной строки

Допустим, мы хотим запустить «чистый компилятор», не пользуясь IDE. Пишем программу в текстовом файле, например Hello.f08:

program Hello
    print *, "Здавствуй, мир! Hello, Space Boy!"
    call execute_command_line("Pause")
endprogram

Компилируем из командной строки:

gfortran Hello.f08 -o Hello.exe

Если пути в PATH «подхватились», после этого получается файл Hello.exe, запуск которого приведет к выводу приветствия на консоль (или «кракозябров» из-за несовпадения кодировки). Если пути в PATH не прописаны, перед gfortran надо указать к нему путь – все точно также будет работать.

Опция -o задает выходной файл.

Чтобы статически сгенерировать exe, не завязанный ни на какие внешние .dll-библиотеки, нужно добавить опцию -static:

gfortran Hello.f08 -static -o Hello.exe

О других опциях – см. ниже.

Далее. Допустим, проект состоит из нескольких модулей. Пусть, например, это три файла: Module1.f08, Module2.f08 и MainProgram.f08 с соответствующими текстами:

module Module1 
    contains 
    function Square(x)  ! Вычисление квадрата x
        real :: x, Square
        Square = x**2
    endfunction
endmodule
module Module2
    use Module1 
    contains 
    subroutine PrintSquare(y)  ! Печать квадрата y 
        real :: y
        print *, "Квадрат числа:", Square(y)
    endsubroutine
endmodule
program MainProgram ! Супер-программа «Нахождение квадрата числа»
    use Module1 ; use Module2
    real :: z = 3.   ! Наше число
    call PrintSquare(z)
    print *, "Проверка:", Square(z)
    call execute_command_line("Pause")
endprogram

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




Компилируется такой проект, чтобы не нарушался порядок зависимостей, следующим образом:

gfortran -c Module1.f08
gfortran -c Module2.f08
gfortran -c MainProgram.f08
gfortran Module1.o Module2.o MainProgram.o -static -o MyProgram.exe

Сначала исходные файлы компилируются в объектные .o (опция -c), а затем они линкуются в результирующий .exe.

Компиляция с GNU Make

Работать с командной строкой, когда фалов проекта много, не слишком удобно (даже если все прописывать в *.bat или *.cmd). Задачу автоматизации сборки берет на себя IDE и/или олдскульная технология составления Make-файлов, которые можно легко вызывать на исполнение, например, из любимого текстового редактора. Преимущество Make, в частности, в том, что каждый раз будут перекомпилироваться не все файлы проекта, а только те, которые изменились.

Если фортран был поставлен из дистрибутива TDM-GCC или CodeBlocs, утилита Make находится в папке bin под именем mingw32-make.exe. Для простоты скопируем (продублируем) ее под именем просто make.exe. В дистрибутиве от Ляо она сразу называется make.exe.

Если на машине инсталлировано несколько сред программирования, могут стоять и различные мэйки, которые также могут вызываться по команде make не совсем к месту. Убедимся, что так вызывается именно GNU Make, набрав make -v. Это должно выдавать приблизительно следующее:

GNU Make 3.82.90
Built for i686-pc-mingw32
Copyright (C) 1988-2012 Free Software Foundation, Inc.
License GPLv3+: GNU GPL version 3 or later <http://gnu.org/licenses/gpl.html>

Инструкции компиляции пишутся в текстовом файле, который обычно называется Makefile (без расширения) и располагается в папке с исходниками.

Make-файл состоит из правил и переменных. Правила имеют следующий синтаксис:

цель1 цель2 ...: реквизит1 реквизит2 ...
    команда1
    команда2
    ...

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

NB! Строки, в которых записаны команды, должны начинаться с символа табуляции. Пробелы здесь ставить нельзя!!!

Простейший Makefile для программы из одного модуля Hello.f08 (наш первый пример) может выглядеть так:

PRG_NAME = Hello

all: 
    gfortran $(PRG_NAME).f08 -o $(PRG_NAME).exe
    $(PRG_NAME).exe

clean:
    rm -rf *.o *.mod $(PRG_NAME).exe

Здесь вначале определена переменная PRG_NAME, значение которой затем «дословно» подставляется в команды – везде, где написано $(PRG_NAME).

Далее определены две цели: all и clean, без реквизитов, но зато с командами.

Если из рабочего каталога вызвать

make all

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

make

Можно также указать требуемый мейк-файл явно, тогда он может называться как нам угодно:

make -f MyMakeFile all

Запуск Make с другой целью, clean, сотрет все промежуточные файлы проекта:

make clean

Теперь составим чуть более сложный Makefile для нашего второго примера с тремя модулями. При этом предусмотрим возможность управления вариантами компиляции с разными опциями, путем их комментирования/раскомментирования. Пусть также exe помещается в отдельную папку bin, которая должна лежать рядом с папкой для исходного кода.

# 1. Флаги компилятора, меняемые оперативно (отладка, оптимизация):

# Быстрая компиляция без лишних вопросов:
# MY_FLAGS = -w -O0 

# Компиляция для отладки с проверками времени выполнения и предупреждениями:
MY_FLAGS = -Og -Wall -Wno-tabs -fcheck=all

# Компиляция для релиза с оптимизацией:
# MY_FLAGS = -w -Ofast -mfpmath=sse -ftree-vectorize -funroll-all-loops

# 2. Флаги компилятора - общие для всех модулей и вариантов компиляции:
FLAGS = -ffree-line-length-512 -freal-4-real-8 

# 3. Название программы и исходных файлов (без расш.) в порядке компиляции:
PRG_NAME = ..\\bin\\MyProgram
F1 = Module1
F2 = Module2
F3 = MainProgram

# 4. Компилятор фортран (здесь можно прописать путь):
CF = gfortran

# 5. Главная цель - сформировать программу и потом ее запустить:
all: $(PRG_NAME)
    $(PRG_NAME).exe

# 6. Для этого нужно достичь другой цели - слинковать объектные файлы в exe:
$(PRG_NAME): $(F1).o $(F2).o $(F3).o 
    $(CF) $^ $(FLAGS) $(MY_FLAGS) -o $@.exe 

# 7. А для этого нужно скомпилировать каждый файл из исходников:
$(F1).o: $(F1).f08 
    $(CF) -c $(FLAGS) $(MY_FLAGS) $^ 

$(F2).o: $(F2).f08 
    $(CF) -c $(FLAGS) $(MY_FLAGS) $^ 

$(F3).o: $(F3).f08 
    $(CF) -c $(FLAGS) $(MY_FLAGS) $^ 

# 8. Альтернативная цель - очистка проекта от промежуточных файлов:
clean:
    rm -rf *.o *.mod $(PRG_NAME).exe

Возможно, это не лучший стиль написания Make, но для меня он «родной» и наглядный.

В главной цели all перед командой запуска файла мы указали реквизит, одноименный с названием программы:

all: $(PRG_NAME)
    $(PRG_NAME).exe

Make – как хороший солдат, любой ценой (но желательно наименьшей), должен выполнить поставленную перед ним цель. А заключается она в исполнении команд, при этом обязательным условием перед их исполнением является наличие реквизита. В качестве реквизита основной цели мы указали необходимость достижения «дочерней» цели, правило для которой описали так:

$(PRG_NAME): $(F1).o $(F2).o $(F3).o 
    $(CF) $^ $(FLAGS) $(MY_FLAGS) -o $@.exe 

Оно требует наличия еще трех реквизитов – готовых объектных файлов $(F1).o $(F2).o $(F3).o. Если они имеются, должна быть выполнена команда по их линковке в exe.

Выражение $@ означает «подставить сюда текущую цель», а $^ – «подставить сюда текущий реквизит».

Цели, связанные с данными реквизитами, мы описали так:

$(F1).o: $(F1).f08 
    $(CF) -c $(FLAGS) $(MY_FLAGS) $^ 

Это значит, что для достижения очередной цели надо, чтобы присутствовал файл $(F1).f08, после чего его надо скомпилировать в объектный файл.

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

Здесь находится перевод официальной справки по GNU Make.

Здесь – хороший материал «Эффективное использование GNU Make» от Владимира Игнатова.

Некоторые опции компилятора

При компиляции и линковке GFortran’у можно указать огромную кучу опций – все они описаны в справках по GFortran и GСС. Перечислю здесь некоторые «типовые» (используемые мной в практике):

Опция Пояснение
-static Указывается при линковке. Собирает весь код зависимых библиотек в один exe-файл. Его можно будет отдавать и переносить без дополнительных dll.
-shared-libgcc Наоборот, dll от GCC – отлинковать динамически. exe-файл будет маленький, но без установленных в доступности PATH нескольких «родных» dll ничего работать не будет.
-m32 Если установлен 64-разрядный GFortran, по умолчанию проект компилируется как 64-разрядное приложение, этой опцией можно попросить сделать версию для 32-разрядных машин.
-O0, -O1, -O2 или -O3 или -Ofast Уровни оптимизации, кода. С -O0 сама компиляция будет быстрой (удобно при разработке), но код – медленный. -O3 или -Ofast – для релиза. -Og – оптимизация для отладки.
-mfpmath=sse -ftree-vectorize -funroll-all-loops Дополнительные оптимизации для релиза (но без фанатизма, для современных процессоров можно навернуть больше).
-fdollar-ok -llaipe2 -lneuloop4 -ljcl Опции для подключения библиотек от Ляо. Первая разрешает доллар в именах функций как у него принято, остальные просто подключают библиотеки laipe2, neuloop4, jcl.
-fopenmp Подключает библиотеку OpenMP для параллельных вычислений.
-fimplicit-none Запрещает неявное определение типов по первой букве у переменных (i – integer и т.д.). Эквивалентно указанию в каждой процедуре implicit none, где их можно забыть и поиметь проблемы с переменными, которые забыли определить явно.
-ffree-line-length-512 Позволяет писать длинные строки в коде (с современными мониторами трудно умещать код в 72 символа, считающиеся ранее максимальной шириной кода).
-freal-4-real-8 Считать все вещественные числа типа real(4), которые в фортране по-умолчанию можно задавать как «просто» real, числами с двойной точностью real(8) (8 байт). На мой взгляд, для большинства задач сейчас лучше использовать 8 байт, а с этой опцией будет меньше визуального мусора типа real(8) :: r = 1.0_8, будет просто real :: r = 1.0.
-w Подавлять все синтактические предупреждения
-Wall -Wno-tabs -fcheck=all Наоборот, выводить все предупреждения (кроме абсурдных), а также выполнять все проверки времени выполнения (выходы за пределы индексов в массивах и т.д., и т.п.)
-ggdb Генерировать информацию для отладчика gdb
-pg Генерировать информацию для профайлера gprof

О расширениях исходных файлов

Расширение файлов по традиции отражает стандарт языка: в наших примерах оно f08 – стандарт 2008. Это дань традиции: компилятор без доп. опций понимает f08, f03, f95, f90 одинаково, как файлы в свободном формате. Для старого фиксированного формата используются расширения f, for и др.

Если расширение начинается с большой буквы (F08, F95 и т.п.), а также если оно fpp, то файл будет сначала пропускаться через предпроцессор – в нем можно будет писать директивы компилятору с синтаксисом gcc-Си, например:

#define Compile_R16
! ...

#if defined Compile_R16
    real*16 :: MyReal
#else
    real*8  :: MyReal
#endif 

Подключение математических библиотек к GFortran

BLAS и LAPACK

BLAS и LAPACK – классические библиотеки линейной алгебры, которые не совсем очевидно подключаются к GFortran в Windows. Здесь описан процесс подключения, если вы работаете c Visual Studio с линкером от Майкрософт, который подключает библиотеки *.dll, используя *.lib. Родной линкер GFortran использует библиотеки в архивированном формате *.a.

Как подключить BLAS и LAPACK к GFortran в Windows:

  1. Скачиваем с официального сайта архив lapack-3.7.0.tgz. (На сегодня последняя версия 3.7.0). Распаковываем архив в папку lapack-3.7.0.

  2. В подпапке INSTALL находим файл make.inc.gfortran, копируем его в корневую папку (lapack-3.7.0) и переименовываем его в просто make.inc.

  3. Идем в каталог \BLAS\SRC (например, в файловом менеджере FAR) и набираем в нем команду make (GNU make должен быть подключен, см. выше). После компиляции ряда файлов в корневой папке (lapack-3.7.0) должен появиться файл librefblas.a.

  4. Аналогично, идем в каталог \SRC и набираем в нем команду make all. После компиляции ряда файлов (довольно долгой) в корне должен появиться файл liblapack.a.

  5. Переписываем файлы librefblas.a и liblapack.a туда, где GCC хранит свои библиотеки. В зависимости от установленного дистрибутива это папки вроде \i686-pc-mingw32\lib, \x86_64-w64-mingw32\lib или \MinGW\lib (в CodeBlocks).

  6. Если мы работаем с дистрибутивом разрядности 64 (например от TDM-GCC) и хотим по опции -m32 уметь компилировать наши программы для 32-разрядных машин, весь процесс придется повторить, скомпилировав librefblas.a и liblapack.a с опцией -m32 (придется дописать ее в make-файл), а потом положить 32-разрядные версии библиотек в соответствующую папку (у TDM-GCC это \x86_64-w64-mingw32\lib32).

  7. Все!

Протестируем. Решим систему линейных алгебраических уравнений (СЛАУ)

\[ \mathbf{A}\cdot\mathbf{X}=\mathbf{B}, \]

где \(\mathbf{A}\) – матрица коэффициентов, \(\mathbf{X}\) – вектор неизвестных и \(\mathbf{B}\) – вектор правых частей.

Пусть матрица \(\mathbf{A}\) – плотная, общего вида. Для СЛАУ с такими матрицами в подойдет LAPACK-процедура dgesv. Первая буква «d» означает вариант для вещественных чисел двойной точности, real(8). Пусть \(\mathbf{A}\) имеет размер \(10 \times 10\) с заполнением случайными числами от \(0\) до \(100\), а вектор \(\mathbf{B}\) заполнен числами \(100 \dots 1000\).

Наберем в файле test.f08:

program LAPACK_test
    implicit none 
    integer, parameter      ::  n = 10    ! Размерность матриц
    real(8), dimension(n,n) ::  A         ! Матрица коэффициентов
    real(8), dimension(n,1) ::  B, X      ! Векторы правых частей / результатов
    integer                 ::  i         ! Вспомог. индекс
    integer, dimension(n)   ::  ipiv      ! Вспомог. матрица для процедуры LAPACK
    integer                 ::  info      ! Инфо об успешности решения
    character(100)          ::  FRM       ! Формат для вывода матриц
    !----------------------------------------------------------------------------
    print *, "Проверка LAPACK. Решение СЛАУ А*X=B с плотной матрицей общего вида"
    call random_number(A);  A = 100*A     ! Заполняем матрицу А 
    forall(i = 1:n) B(i,1) = 100*i        ! Заполняем вектор B
    write(FRM, "(a,i10,a)") "(",n,"f8.2)" ! Формируем формат для вывода матриц
    print *, "Матрица коэффициентов A:"
    print FRM, A
    print *, "Вектор правых частей B:"
    print FRM, B
    X = B
    call dgesv(n, 1, A, n, ipiv, X, n, info)  ! Вызов LAPACK для решения СЛАУ
    print *, "Решение:"
    print FRM, X
    print *, "Успешность решения (флаг «info»):", info
endprogram

Компилируем, подключая библиотеки lapack / refblas:

gfortran test.f08 -llapack -lrefblas -static -o test.exe

Получаем:

Проверка LAPACK. Решение СЛАУ А*X=B с плотной матрицей общего вида
Матрица коэффициентов A:
  99.76   56.68   96.59   74.79   36.74   48.06    7.38    0.54   34.71   34.22
  21.80   13.32   90.05   38.68   44.55   66.19    1.61   65.09   64.64   32.30
  85.57   40.13   20.69   96.85   59.84   67.30   45.69   33.00   10.04   75.55
  60.57   71.90   89.73   65.82   15.07   61.23   97.87   99.91   25.68   55.09
  65.90   55.40   97.78   90.19   65.79   72.89   40.25   92.86   14.78   67.45
  76.96   33.93   11.58   61.44   82.06   94.71   73.11   49.76   37.48   42.15
  55.29   99.79   99.04   74.63   95.38    9.33   73.40   75.18   94.68   70.62
  81.38   55.86    6.17   48.04   59.77   13.75   58.74   52.00   88.59   30.38
  66.97   66.49   50.37   26.16    7.66   10.12   54.93   37.56    1.51   79.29
  62.09   77.36   95.36   11.42   31.85   59.68    4.82   11.42   21.60   10.06
Вектор правых частей B:
 100.00  200.00  300.00  400.00  500.00  600.00  700.00  800.00  900.00 1000.00
Решение:
  -8.79   16.14   12.52    3.00  -14.42   -0.18    4.31   -2.40    5.41   -1.10

Библиотеки дядюшки Ляо

В дистрибутиве GCC профессора Ляо (с сайта www.equation.com) имеются очень неплохие библиотеки:

  • NEOLOOP – для параллельных вычислений;

  • LAIPE2 – решение СЛАУ с матрицами коэффициентов различного вида, включая симметричные положительно определенные ленточные или просто разряженные, представляющие для нас особый интерес (почему – об этом потом). Решение выполняется с распараллеливанием процессов.

  • JCL – перенумерация строк и столбцов разряженных матриц для приведения их к ленточной форме.

Если нужно установить эти библиотеки в другие дистрибутивы, делаем следующее:

  1. Находим в дистрибутиве от Ляо файлы библиотек: libneuloop4.a, libneuloop6.a, liblaipe2.a, libjcl.a.

  2. Переписываем их в место, где хранит основные библиотеки ваш фортран.

Как обычно, если мы работаем с 64-битной версией GCC с желанием сохранить возможность компиляции программ под 32, это надо сделать для обеих версий библиотек, положив их в lib и lib32.

Компиляция программ с этими библиотеками выполняется с опциями:

-fdollar-ok -llaipe2 -lneuloop4 -ljcl

(Первая разрешает доллар в именах функций как у него принято, остальные просто подключают библиотеки laipe2, neuloop4, jcl).

Чтобы подпрограммы LAIPE2 заработали, надо предварительно вызвать процедуру laipe$use с указанием числа ядер процессора, которые вы хотите использовать (иначе не «заведется» – об этом Ляо забыл написать в справке).

Библиотеки NEOLOOP и LAIPE2 работают только на Windows 7 и более старших.

IDE и редакторы для GFortran

Для серьезной работы с языком хорошо иметь как IDE, со всем богатством возможностей, так и текстовый редактор (может быть, не один), позволяющий «что-нибудь сварганить» по-быстрому или подправить. Честно говоря, я IDE почему-то запускаю реже текстового редактора, с которым намного лучше постигается Дзен.

Какие есть варианты для GFortran в Windows?

IDE Code::Blocks

Это самый простой вариант начать работать в Fortran, пользуясь свободным ПО. Все ставится из коробки и заводится с пол-оборота. По сравнению с альтернативными IDE – относительно легкая. Тут есть все что нужно: автодополнение кода, подсказки, навигация по объектам программы, отладчик, профайлер, минимальные средства рефакторинга.


Когда-то я «щупал» ifort + Visual Studio, там таких удобных подсказок и умного автодополнения и близко не было (это не претензия к студии – в других языках все очень круто, но не в Intel Fortan). Т.е. в определенном смысле CodeBlocks – даже лучше. (Может конечно я отстал от жизни, последних версий ifort + VS не видел).

На 64-разрядных машинах желательно поставить 64-разрядный компилятор GFortran (см. выше) и подключить его в качестве основного (все делается мышкой в настройках элементарно).

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

Executable: ${TARGET_OUTPUT_BASENAME}
Working directory: ${TARGET_OUTPUT_DIR}
Launch tool hidden with standart output redirect

IDE Eclipse с расширением Photran

Eclipse – знаменитая мощнейшая IDE для множества языков программирования. Для Fortran есть расширение Photran, работающее в том числе c GCC. Доступно здесь.

Eclipse – очень тяжелая штука. Я побаловался, но дальше она как-то не прижилась.

Прочие IDE

Еще можно посмотреть в сторону легонькой Geany.

Не уверен, что своими силами можно толком настроить GFortran c Visual Studio. Есть коммерческий пакет, где это сделано – Lahey/Fujitsu Fortran. Это не наш путь.

Текстовые редакторы

Не могу отделаться от ассоциации программистских текстовых редакторов с людьми разных возрастных групп. :)

Emacs и Vim

Седовласые старики. Очень опытные, мудрые и крутые, но найти с ними общий язык – очень сложно. Я пока не сумел, хотя пару десятков попыток делал. :)

Sublime Text, Atom, VS Code

Энергичные мужики. Моложавые, хипстерской внешности. На этих редакторах получается сделать «почти IDE».

Atom – открытый, мощный, но жирный и неторопливый. Использую его для некоторых других языков, Фортран как-то на нем не прижился.

Sublime Text – мой выбор. Мощный, быстрый, правда не открытый и не бесплатный (70 долларов, хотя незарегистрированной версией можно пользоваться сколько угодно).


VS Code – удачная попытка Microsoft написать бесплатный Sublime. Очень приятный. Все намереваюсь на него перейти, но руки пока не доходят.

SciTE, NotePad++, AkelPad

Эдакие продвинутые юнцы. Меленькие, легкие, дерзкие, но далеко не такие мощные, как товарищи из первых двух групп. Не «почти IDE».

Раньше я очень любил SciTE, последнее время перешел в AkelPad. Он крошечный, работает со сверхзвуковой скоростью. Очень просто настраивается подсветка. Все, чего только только можно желать, реализуется на скриптах на простом человеческом JavaScript.


С GFortran нетрудно настроить любой редактор. Привожу некоторые скриншоты работы приведенной выше программы по тестированию правильности установки LAPACK.

Литература по Fortran

Лучшая книга на русском языке – О. В. Бартеньев «Современный Фортран». Она настолько полная и и исчерпывающая, что книги остальных авторов – не нужны. Хорош также его же трехтомник по библиотеке IMSL. Даже если вы не применяете IMSL, там много полезной теории. (На самом деле книги Бартеньева – по сути качественный перевод полной документации по Фортрану и IMSL и от минимума «отсебятины» они только лучше). Жаль только, что его книги несколько устарели, но это совсем не критично – Фортран развивается не так быстро.

Здесь – GFortran Wiki.

Здесь – PDF-справка по GFortran. Очень неплоха. Коротко и по делу. Для работы нужен второй PDF по GCC в целом.

Здесь – очень хорошее краткое описание языка в стандарте 95 в википедии. Если вы пишете на других языках, прочитав это, станете гуру в Фортране за каких-нибудь пару часов. :)

Здесь, в открытом доступе, – прекрасная справка по Intel Fortran, в основной части (в части самого языка) подходящая и для GFortran.

Здесь, аналогично, – прекрасная документация по Intel MKL, в основной части подходящая для BLAS и LAPACK.

Что прекрасно в Fortran

  1. Работа с матрицами.

С матрицами (например, A, B, C) можно работать как с переменными:

A = 2*B + sin(C)

К матрицам (поэлементно) могут быть применены любые чистые функции (как синус в примере).

Вырезки и сечения матриц:

B[1:5,1:5] = A[1:10:2,1:10:2]

В отличие от MatLAB-подобных языков, шаг - третий, а не второй элемент триплета.

Конструкция forall. Например, заполним матрицу \(n \times n\) единицами по диагонали (сделаем ее единичной):

A = 0 
forall(i = 1:n) A(i,i) = 1

Конструкция where. Например, так можно обнулить все отрицательные компоненты матрицы B:

where(B < 0) B = 0

С матрицами, вообще, множество удобных встроенных функций – практически на все случаи жизни.

Матрицы могут быть многомерными.

  1. Вызов процедур: поддерживается перегрузка параметров, необязательные и именованные параметры.

  2. Типы, включая комплексные числа. Удобные древовидные пользовательские типы (структуры данных).

  3. Как ни странно, форматный ввод-вывод, к которому привыкаешь, но который я также отмечу как одно из самых ужасных свойств языка. :)

Что ужасно в Fortran

Современный Fortran позволяет писать более-менее сносный код, если не пользоваться устаревшими конструкциями языка, поддерживаемыми ради совместимости, например:

  • типизацией переменных по умолчанию на основе первых букв их имен;

  • разными точками входа в процедуры;

  • переходами goto с метками;

  • указателями и прочими средствами работы с памятью «напрямую» и т.д.

Иногда чешутся руки, наоборот, это все использовать и писать программы аутентично «индийские». В идеале можно выйти на абсолютный Brainfuck. :)

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

  • Работа со строками.

В Fortran нет строк переменной длины, размер строк нужно задавать явно.

Специфичны преобразования других типов данных в строку и обратно. Чтобы преобразовать, например, вещественное число r в строку str или наоборот, приходится делать это через операторы ввода-вывода:

write(str, "(f8.3)") r  ! str <- r в формате f8.3
read(str,*) r           ! str -> r

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

str = trim(str1) // trim(str2)  ! Если строки выровнены влево
str = trim(AdjustL(str1)) // trim(AdjustL(str2)) ! Если не выровнены :)

Разумеется, в Fortran нет никаких регулярных выражений и других «вкусностей», к которым в современных языках мы успели уже «прикипеть».

  • Форматный ввод-вывод – очень мясной. Он был удобен для терминалов, но сейчас это – кошмарный анахронизм:
print "(a,i4,a,10f8.3)", "Шаг:", i, "Вектор:", V

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

  • Нет списков и вообще структур данных с изменяемой длиной «на лету». (Динамически размещаемые массивы, конечно же, есть).

  • Переусложненная система типов. Вещественные числа разной длины пишутся с указанием разновидности (KIND) весьма специфично:

real     :: r1 = 1.     ! По умолчанию real – это real(4)
real(8)  :: r2 = 1._8
real(16) :: r3 = 1._16

С другой стороны, KIND можно определить на уровне директивы компилятора и везде использовать просто real и числа без этих «_8»

  • Мало средств «декомпозиции и абстракции» для поддержания структуры больших проектов.

  • Избыточность языка, связанная с историей и совместимостью снизу вверх.

  • Наследуемый код, написанный за пол века «учеными и инженерами» – в подавляющей массе просто чудовищен. Чуть менее, чем весь созданный код. :) Почему так – загадка. Прямой вины языка в этом нет. Потемки – душа ученого.

Вместо Conclusions

В этой импровизированной статье пока не был раскрыт ключевой вопрос абстракта «как заставить себя в XXI веке писать на Fortran?». :)

Мейнстримовый C++ предлагает аналогичную скорость.

Лаконичный Python предоставляет неслыханную мощь и ясность/читаемость кода. Одна строчка идет за 5-10 Фортрана.

Julia совсем наступает на пятки, предлагая почти ту же мощь и прозрачность, что и Python, и вполне сравнимую с Fortran скорость.

На самом деле, у меня нет разумного ответа на сакраментальный вопрос, кроме лишь одного: Fortran надо любить. :) Как дедушку и стареющего отца.

Каждый раз, открывая текстовый редактор с Фортраном, чувствуешь, что «подключаешься» к чему-то ретрофутуристично-прекрасному. :) Будущее, как видели его в 60-х годах, космические полеты, покорение звезд. Элементарные частицы, атомные реакторы. Шуршащие пленкой вычислительные машины. Папиросы профессоров «Прикладной математики». Все это преломляется в вашем редакторе, хочется открывать его еще и еще. :) А на современных машинах все эти отблески проносятся с бешенной скоростью… Ну и конечно, по совокупности качеств, с чего «статья» начиналась…

Пускай в этом блоге Fortran будет отправной точкой в мир более современных решений.