Матрицы, matmul, dot_product и BLAS/LAPACK

Fortran делает линейную алгебру естественной: операции над массивами, встроенные matmul и dot_product, а для серьёзных задач — оптимизированные BLAS и LAPACK.

BLAS (Basic Linear Algebra Subprograms) — стандартный набор низкоуровневых операций линейной алгебры; LAPACK — построенная над ним библиотека решения линейных систем, собственных значений и разложений матриц. Обе — эталон производительности численных вычислений.

Почему линейная алгебра — родная стихия Fortran

Огромная доля научных вычислений сводится к линейной алгебре: решение систем уравнений, метод наименьших квадратов, собственные значения, разложения матриц. Fortran проектировался с массивами как первоклассными объектами, и это делает запись линейной алгебры предельно близкой к математике. Двумерный массив — это матрица, и операции над ним пишутся естественно. Более того, исторически именно на Fortran написаны эталонные численные библиотеки (BLAS, LAPACK), которыми сегодня пользуется весь мир, включая NumPy и MATLAB под капотом. Понимать линейную алгебру в Fortran — значит понимать фундамент современных научных вычислений.

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

Поэлементные операции (+, -, *, /) применяются к массивам покомпонентно. Важно: a * b для двух матриц — это поэлементное умножение (произведение Адамара), а не матричное! Матричное произведение — это встроенная функция matmul. Эту разницу новички путают постоянно.

program matrix_basics
  implicit none
  real :: a(2,2), b(2,2), c(2,2)
  a = reshape([1.0, 2.0, 3.0, 4.0], [2,2])   ! по столбцам!
  b = reshape([5.0, 6.0, 7.0, 8.0], [2,2])

  c = a * b           ! ПОЭЛЕМЕНТНО (Адамар)
  print *, "a*b (поэлементно), c(1,1)=", c(1,1)

  c = matmul(a, b)    ! настоящее матричное умножение
  print *, "matmul(a,b), c(1,1)=", c(1,1)

  print *, "транспонированная a, элемент (1,2)=", transpose(a)
end program matrix_basics

Вывод:

 a*b (поэлементно), c(1,1)=   5.00000000    
 matmul(a,b), c(1,1)=   23.0000000    
 транспонированная a, элемент (1,2)=   1.00000000       3.00000000       2.00000000       4.00000000

Обратите внимание на reshape([1,2,3,4], [2,2]): Fortran заполняет матрицу по столбцам (column-major), поэтому элементы ложатся как a(1,1)=1, a(2,1)=2, a(1,2)=3, a(2,2)=4 — то есть первый столбец это (1,2), второй (3,4). Это фундаментальное свойство Fortran, отличающее его от C (row-major), и оно влияет и на reshape, и на производительность (см. урок про кэш).

Встроенные операции линейной алгебры

Fortran предоставляет ключевые операции встроенными функциями, что и компактно, и быстро:

Скалярное произведение через dot_product и норму через norm2 предпочитают ручным циклам: они и читаются лучше, и оптимизированы. Для матрицы на вектор тоже работает matmul: y = matmul(A, x).

program linalg_ops
  implicit none
  real :: u(3) = [1.0, 2.0, 2.0]
  real :: a(2,3) = reshape([1.,0., 0.,1., 1.,1.], [2,3])
  real :: x(3) = [1.0, 1.0, 1.0]

  print *, "dot_product(u,u) =", dot_product(u, u)   ! 1+4+4 = 9
  print *, "norm2(u)         =", norm2(u)            ! sqrt(9) = 3
  print *, "matmul(a,x)      =", matmul(a, x)        ! [2, 2]
end program linalg_ops

Вывод:

 dot_product(u,u) =   9.00000000    
 norm2(u)         =   3.00000000    
 matmul(a,x)      =   2.00000000       2.00000000

Когда встроенного matmul недостаточно: BLAS

Встроенный matmul хорош для умеренных размеров, но для больших матриц (сотни-тысячи строк) он, как правило, проигрывает специализированным библиотекам. BLAS — это стандартизованный набор процедур линейной алгебры, реализации которого (OpenBLAS, Intel MKL, ATLAS) доведены до предела: они учитывают иерархию кэшей, используют векторные инструкции (AVX), блочные алгоритмы и многопоточность. Разница в скорости на больших матрицах — десятки раз. BLAS делится на три уровня: уровень 1 — операции вектор-вектор (ddot, daxpy); уровень 2 — матрица-вектор (dgemv); уровень 3 — матрица-матрица (dgemm). Именно уровень 3, и прежде всего dgemm (умножение матриц), даёт максимальный выигрыш, потому что переиспользует данные в кэше.

! Вызов BLAS dgemm: C := alpha*A*B + beta*C (double precision)
! Интерфейс обычно подключают через модуль или внешний interface.
external :: dgemm
real(8) :: A(m,k), B(k,n), C(m,n)
real(8) :: alpha = 1.0d0, beta = 0.0d0
call dgemm('N', 'N', m, n, k, alpha, A, m, B, k, beta, C, m)
! 'N','N' — матрицы не транспонированы; m,k — ведущие размерности

Интерфейс BLAS архаичен (много позиционных аргументов, ведущие размерности lda), но это плата за десятилетия совместимости и оптимизации. В современном коде поверх BLAS часто кладут удобные обёртки, но под капотом работает именно dgemm.

Решение систем и разложения: LAPACK

Когда нужно решить систему Ax=b, найти собственные значения или сделать LU/QR/SVD-разложение, используют LAPACK — библиотеку поверх BLAS. Решать систему «в лоб» обращением матрицы (x = A⁻¹b) численно плохо: дорого и неустойчиво. LAPACK предлагает правильные алгоритмы. Решение линейной системы — это dgesv (LU-разложение с частичным выбором главного элемента).

! Решить A*x = b методом LU (LAPACK dgesv), double precision
external :: dgesv
integer :: n, nrhs, lda, ldb, info
integer :: ipiv(n)
real(8) :: A(n,n), b(n,1)
n = 3; nrhs = 1; lda = n; ldb = n
! ... заполнить A и b ...
call dgesv(n, nrhs, A, lda, ipiv, b, ldb, info)
if (info == 0) then
  print *, "Решение в b:", b(:,1)   ! dgesv пишет решение в b
else
  print *, "Матрица вырождена или ошибка, info =", info
end if

Здесь dgesv на месте превращает A в LU-разложение, а b — в решение x. Выходной параметр info сообщает успех (0) или проблему (вырожденность, неверный аргумент). Этот шаблон — «вызвать LAPACK и проверить info» — основа надёжного численного кода. LAPACK покрывает практически всю плотную линейную алгебру: dgels (наименьшие квадраты), dsyev (собственные значения симметричной матрицы), dgesvd (сингулярное разложение) и десятки других.

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

Производительность BLAS уровня 3 держится на повторном использовании данных. Наивное умножение матриц n×n выполняет O(n³) операций над O(n²) данными — значит каждый элемент в среднем участвует в n операциях. Если организовать вычисление блоками, помещающимися в кэш, эти n операций над элементом происходят, пока он горяч в кэше, а не подгружается из медленной памяти заново. Оптимизированный dgemm именно так и устроен: разбивает матрицы на блоки под размер кэша L1/L2, применяет векторные инструкции к блокам и распараллеливает по ядрам. Встроенный matmul компилятора обычно проще и для больших размеров отстаёт. LAPACK же строит сложные алгоритмы (LU, QR) из блочных вызовов BLAS уровня 3, наследуя его эффективность. Вот почему практическое правило численного программиста: не пишите собственное умножение больших матриц и не решайте системы обращением — вызывайте BLAS/LAPACK, за которыми стоят сотни человеко-лет оптимизации.

Почему мир считает на Fortran-библиотеках

Стоит осознать историческую и практическую глубину связи Fortran с численной линейной алгеброй, потому что она объясняет, почему эти библиотеки до сих пор эталон. BLAS появился в конце 1970-х, LAPACK — в начале 1990-х, оба изначально на Fortran. За десятилетия в их реализации вложены колоссальные усилия лучших специалистов по численным методам и оптимизации. Результат настолько хорош, что весь остальной мир научных вычислений построен поверх него. Когда вы пишете в Python numpy.dot(A, B) или numpy.linalg.solve, под капотом вызывается BLAS/LAPACK (через OpenBLAS или Intel MKL). MATLAB, R, Julia — все они опираются на те же библиотеки. То есть даже программист, никогда не видевший Fortran, ежедневно пользуется плодами Fortran-экосистемы. Это уникальное положение: язык, который многие считают «устаревшим», на деле остаётся фундаментом, на котором стоит современная вычислительная наука и значительная часть машинного обучения.

Практический вывод для вас как для Fortran-программиста двоякий. С одной стороны, вы имеете прямой доступ к этим библиотекам, без обёрток и накладных расходов чужого языка — ваш численный код может быть настолько же быстр, как самый быстрый код в мире, потому что он вызывает те же dgemm и dgesv напрямую. С другой стороны, это налагает ответственность: раз эталонные реализации существуют и доступны, писать своё умножение матриц или свой решатель систем для production-кода — почти всегда ошибка. Ваша задача — грамотно скомпоновать вызовы оптимизированных библиотек, а не переизобретать их. Исключение — учебные цели и очень специфические структуры данных, где общие библиотеки неприменимы.

Структурированные матрицы и выбор алгоритма

LAPACK содержит не одну, а семейство процедур для каждой задачи, и выбор правильной даёт огромный выигрыш. Ключ — структура матрицы. Матрица общего вида, симметричная, положительно определённая, ленточная, треугольная, разреженная — для каждой есть свой оптимальный алгоритм, и LAPACK кодирует это в именах процедур. Префикс задаёт тип данных (s/d/c/z — single/double/complex/double-complex), средние буквы — тип матрицы (ge general, sy symmetric, po positive-definite, gb general band), окончание — операцию (sv solve, trf factorize, ev eigenvalues).

Например, для симметричной положительно определённой матрицы (типичной в задачах метода конечных элементов и наименьших квадратов) разложение Холецкого dposv вдвое дешевле общего LU и численно устойчивее — грех не воспользоваться. Для ленточной матрицы (возникает при дискретизации дифференциальных уравнений) специальный решатель dgbsv хранит и обрабатывает только ненулевую ленту, экономя и память, и время на порядки по сравнению с трактовкой её как плотной. Поэтому квалифицированный вычислитель сначала анализирует структуру своей задачи, а затем выбирает соответствующую процедуру LAPACK. Это и есть настоящее искусство численного программирования на Fortran: не писать алгоритмы (они уже написаны и отточены), а правильно подобрать и скомпоновать готовые, понимая математическую структуру задачи. Тот, кто умеет сопоставить структуру матрицы с нужной процедурой LAPACK, получает и максимальную скорость, и максимальную точность — и именно так выглядит компетентный современный научный код на Fortran.

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

• Путать * и matmul. a * b — поэлементное умножение, matmul(a,b) — матричное. Это разные операции с разным результатом.
• Решать систему обращением матрицы. x = matmul(inverse(A), b) дороже и менее устойчиво, чем dgesv. Используйте LU-решатель LAPACK.
• Писать своё умножение больших матриц. Самописный тройной цикл проиграет dgemm в разы. Для серьёзных размеров — BLAS.
• Игнорировать info у LAPACK. Ненулевой info означает вырожденность или ошибку аргумента; без проверки получите мусор как «решение».
• Забыть про column-major при reshape и интеропе. Fortran заполняет и хранит по столбцам; данные из C (по строкам) нужно транспонировать.

Итоги

• Поэлементное a*b ≠ матричное matmul(a,b) — не путать; для скаляра векторов есть dot_product.
• Fortran хранит и заполняет массивы по столбцам (column-major) — важно для reshape и интеропа.
• Встроенные функции хороши для умеренных размеров; для больших матриц — BLAS (особенно dgemm, уровень 3).
• Систему Ax=b решают LU-методом через LAPACK dgesv, а не обращением матрицы; всегда проверяют info.
• Скорость BLAS/LAPACK — это блочные кэш-оптимизированные алгоритмы; своё умножение/решатель писать не нужно.