Операции над массивами целиком

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

Операции над массивами целиком (whole-array operations) — применение арифметических операторов и функций сразу ко всем элементам массива поэлементно, без явного цикла, что выражает математику напрямую и позволяет компилятору векторизовать код.

Вот ради чего инженеры выбирают Fortran. Там, где в C нужно написать цикл, объявить индекс и следить за границами, в Fortran вы пишете c = a + b для целых векторов — и язык сам выполняет сложение поэлементно. Это не просто короче: код выражает математическую идею напрямую, а компилятор получает свободу превратить операцию в эффективные векторные инструкции процессора. Этот урок — о массивной арифметике, поэлементных функциях, broadcasting и встроенных операциях вроде суммы и скалярного произведения.

Поэлементная арифметика

Арифметические операторы +, -, *, /, ** применяются к массивам одинаковой формы поэлементно. Результат — массив той же формы.

program elementwise
  implicit none
  real :: a(4), b(4), c(4)
  a = [1.0, 2.0, 3.0, 4.0]
  b = [10.0, 20.0, 30.0, 40.0]
  c = a + b                 ! поэлементное сложение
  print *, "a + b =", c
  c = a * b                 ! ВНИМАНИЕ: поэлементное, не матричное!
  print *, "a * b =", c
  c = b / a
  print *, "b / a =", c
end program elementwise

Вывод:

 a + b =   11.0000000  22.0000000  33.0000000  44.0000000
 a * b =   10.0000000  40.0000000  90.0000000  160.000000
 b / a =   10.0000000  10.0000000  10.0000000  10.0000000

Критически важно: a * b — это поэлементное умножение (произведение Адамара), а не матричное произведение и не скалярное. Каждый элемент результата — произведение соответствующих элементов. Для настоящего матричного умножения есть функция matmul, для скалярного — dot_product (ниже). Путаница здесь — частый источник ошибок у тех, кто ждёт «как в линейной алгебре».

Broadcasting: массив и скаляр

Когда один операнд — скаляр, а другой — массив, скаляр «распространяется» (broadcast) на все элементы. Это позволяет масштабировать и сдвигать массивы предельно лаконично.

program broadcasting
  implicit none
  real :: celsius(4), fahrenheit(4)
  celsius = [0.0, 25.0, 50.0, 100.0]
  fahrenheit = celsius * 1.8 + 32.0    ! скаляр применяется ко всем
  print *, "По Цельсию:  ", celsius
  print *, "По Фаренгейту:", fahrenheit
end program broadcasting

Вывод:

 По Цельсию:     0.00000000  25.0000000  50.0000000  100.000000
 По Фаренгейту:   32.0000000  77.0000000  122.000000  212.000000

Выражение celsius * 1.8 + 32.0 умножает каждый элемент на 1.8 и прибавляет 32 — целая формула перевода температур в одну строку для всего массива. Скаляр в операции с массивом всегда трактуется так, будто он повторён до нужной формы.

Элементные (elemental) функции

Встроенные математические функции — sqrt, sin, exp, abs и десятки других — являются элементными: применённые к массиву, они действуют на каждый элемент и возвращают массив. Это устраняет циклы при векторных вычислениях.

program elemental_funcs
  implicit none
  real :: x(5), y(5)
  integer :: i
  x = [(real(i), i = 1, 5)]
  y = sqrt(x)              ! корень от каждого элемента
  print *, "x   =", x
  print *, "sqrt=", y
  print *, "sin =", sin(x)
  print *, "Сумма квадратов:", sum(x**2)
end program elemental_funcs

Вывод:

 x   =   1.00000000  2.00000000  3.00000000  4.00000000  5.00000000
 sqrt=   1.00000000  1.41421354  1.73205078  2.00000000  2.23606801
 sin =  0.841470957 0.909297407  0.141120002 -0.756802499 -0.958924294
 Сумма квадратов:   55.0000000

Запись sqrt(x) возвращает массив корней, sin(x) — массив синусов. Это позволяет переносить математические формулы в код почти буквально: где в учебнике стоит √x для вектора, там в Fortran пишется sqrt(x).

Встроенные операции редукции

Fortran предоставляет функции, «сворачивающие» массив в одно значение или меньший массив — редукции. Это рабочий инструмент численного анализа.

program reductions
  implicit none
  real :: v(5)
  v = [3.0, 1.0, 4.0, 1.0, 5.0]
  print *, "Сумма:", sum(v)
  print *, "Максимум:", maxval(v), " на позиции", maxloc(v)
  print *, "Среднее:", sum(v) / size(v)
  print *, "Скалярное произведение v·v:", dot_product(v, v)
end program reductions

Вывод:

 Сумма:   14.0000000
 Максимум:   5.00000000  на позиции           5
 Среднее:   2.79999995
 Скалярное произведение v·v:   52.0000000

Редукции — это мост между «сырым» массивом и числом, которое нужно человеку или следующему этапу расчёта: норма вектора, энергия системы, невязка, среднее по выборке. Почти всё это выражается через sum, dot_product и элементные функции в одну-две строки. Например, евклидова норма вектора — это sqrt(sum(v**2)) или, элегантнее, sqrt(dot_product(v, v)); среднеквадратичное отклонение прогноза от факта — sqrt(sum((pred - obs)**2) / size(obs)). Эти выражения читаются почти как формулы из учебника, и именно в этом ценность: между математической записью и кодом почти нет «переводческого» слоя, где обычно и заводятся ошибки.

Почему математика ложится в код один в один

Главная мысль этого урока шире, чем синтаксис: массивная нотация устраняет разрыв между тем, как задачу формулирует математик, и тем, как её записывает программист. В обычном языке формула вроде «новое поле = старое поле плюс шаг по времени, умноженный на производную» превращается в цикл с индексами, проверкой границ и временными переменными — и читателю кода приходится мысленно «собирать» формулу обратно из этих деталей. В Fortran та же формула пишется как u_new = u + dt * du для целых массивов, и она остаётся формулой. Это не косметика: код, который выглядит как математика, легче проверять на соответствие выводу на бумаге, легче сопровождать и труднее испортить опечаткой в индексе, потому что индексов попросту нет.

Отсюда вытекает практическое правило стиля в численном Fortran: предпочитайте массивную запись явным циклам везде, где это выражает намерение. Цикл уместен, когда логика по-настоящему поэлементная и нерегулярная (сложная ветвистая обработка каждого элемента), но для арифметики над полями целиком массивная форма и короче, и быстрее, и ближе к предметной области. Опытные расчётчики читают c = a + b мгновенно как «сложить векторы», тогда как трёхстрочный цикл заставляет их останавливаться и проверять, не закралась ли в границы ошибка.

Адамар, скаляр и матрица: три разных умножения

Стоит ещё раз заострить различие, потому что оно — источник самой частой ошибки новичка. Для двух векторов одинаковой длины в Fortran существуют три принципиально разные операции «умножения». Поэлементное a * b даёт вектор, где перемножены соответствующие элементы, — это произведение Адамара, полезное, например, при применении маски или весов. Скалярное dot_product(a, b) даёт одно число — сумму попарных произведений, то есть привычное из линейной алгебры скалярное произведение. А матричное matmul работает уже с матрицами (и матрица-вектор), реализуя правило «строка на столбец». Знак * в Fortran — это всегда поэлементная операция, и язык не делает исключения для «похожих на матрицы» массивов; ожидание иного — наследие математической нотации, где AB по умолчанию означает матричное произведение.

matmul, BLAS и почему важны библиотеки

Функция matmul заслуживает отдельного слова, потому что за её скромным вызовом стоит целая индустрия оптимизации. Наивное матричное умножение — это тройной вложенный цикл со сложностью порядка n³, и написанный «в лоб» он чудовищно недогружает процессор из-за промахов кэша. Поэтому компиляторы Fortran для больших матриц нередко перенаправляют matmul в специализированную библиотеку семейства BLAS, где тот же алгоритм реализован с блочной разбивкой под размер кэша, векторизацией и иногда многопоточностью. Разница в скорости между наивной и библиотечной реализацией — порядки, и именно поэтому в серьёзном расчётном коде матричные операции стараются выражать через стандартные вызовы, а не переписывать руками. Для совсем тяжёлых случаев расчётчики идут дальше matmul и вызывают BLAS/LAPACK напрямую — но даже встроенная функция уже даёт доступ к этому миру оптимизаций бесплатно.

Полезно сравнить и с тем, как устроена скорость в Python: там numpy.dot и оператор @ тоже в конечном счёте вызывают BLAS, написанную на C и Fortran. То есть когда питонист пишет быстрый матричный код, он опосредованно пользуется фортрановской машинерией. Разница в том, что в Fortran этот слой — родной и не отделён от остального кода границей «интерпретатор — нативная библиотека», через которую данные приходится упаковывать и распаковывать. Для алгоритма, где много мелких разнородных операций над массивами, отсутствие этой границы и даёт Fortran преимущество: нет постоянных переходов туда-обратно, всё компилируется в единый поток инструкций.

Как работает под капотом

Почему массивная запись не только короче, но и быстрее? Когда вы пишете c = a + b, компилятор знает, что это поэлементная операция над непрерывными блоками одинаковой формы, и итерации заведомо независимы. Это идеальный кандидат для векторизации: современный процессор имеет SIMD-инструкции, обрабатывающие сразу 4, 8 или 16 чисел за один такт. Ручной цикл на C компилятор тоже иногда векторизует, но только доказав отсутствие алиасинга; массивная нотация Fortran объявляет независимость явно, упрощая оптимизатору жизнь. Кроме того, операции вроде sum и matmul часто реализованы внутри компилятора оптимально — matmul может вызывать высокопроизводительную BLAS-подобную реализацию. Важная тонкость: в выражении c = a + b * d компилятор не обязан создавать видимый временный массив для b * d — он вправе слить операции в один проход (fusion), экономя память и проходы по данным. Именно эта свобода и делает Fortran быстрым на числодробилках.

Частые ошибки

• a * b как матричное умножение. Это поэлементное произведение! Для матричного — matmul, для скалярного — dot_product.
• Несовпадение форм. Поэлементные операции требуют одинаковой формы массивов; иначе ошибка (выявляется компилятором или в рантайме с -fcheck).
• Ручные циклы вместо массивной нотации. Работают, но длиннее и хуже оптимизируются; предпочитайте c = a + b.
• Среднее как sum(v)/size(v) с целыми. Если массив целый, деление целочисленное; приводите к real.
• Ожидание порядка суммирования. sum может суммировать в любом порядке (для скорости), что влияет на последний бит при сложении вещественных.

Итоги

• Операторы + - * / ** применяются к массивам одной формы поэлементно.
• a * b — поэлементное (Адамарово) произведение, НЕ матричное; для матричного — matmul.
• Скаляр в операции с массивом распространяется (broadcast) на все элементы.
• Встроенные функции (sqrt, sin, exp) элементны: к массиву применяются поэлементно.
• Редукции (sum, maxval, dot_product) сворачивают массив в значение.
• Массивная нотация позволяет компилятору векторизовать код и часто быстрее ручных циклов.