Умножение матриц

Умножение двух матриц -- фундаментальная операция линейной алгебры, используемая в компьютерной графике, машинном обучении, решении систем уравнений и множестве других задач. Программа вычисляет произведение $C = A \cdot B$ для нижнетреугольной и верхнетреугольной матриц и проверяет результат по аналитической формуле.

Задача

Вычислить произведение двух матриц $A$ и $B$ размера $n \times n$:

$C_{i,j} = \sum_{k=0}^{n-1} A_{i,k} \cdot B_{k,j}, \quad i,j = 0, \ldots, n-1$

В примере используются:

$A$ -- нижнетреугольная матрица: $A_{i,k} = 1$ при $k \le i$, иначе $0$
$B$ -- верхнетреугольная матрица: $B_{k,j} = 1$ при $k \le j$, иначе $0$

Ожидаемый результат: $C_{i,j} = \min(i, j) + 1$ -- симметричная матрица.

Идея алгоритма

Заполнить матрицу $A$ (нижний треугольник -- единицы) и матрицу $B$ (верхний треугольник -- единицы).
Вычислить произведение $C = A \cdot B$ тройным вложенным циклом по $i$, $j$, $k$.
Проверить результат: для каждого элемента вычислить ошибку $e_{i,j} = C_{i,j} - C_{i,j}^{\text{exact}}$ и вывести нормы $|e|2$ и $|e|\infty$.

Разбор

Объявление данных

Задаём размер $n = 128$ и объявляем три двумерных массива для матриц $A$, $B$ и результата $C$:

  цел n,i,j,k
  n := 128
  вещ таб A[0:n-1,0:n-1]
  вещ таб B[0:n-1,0:n-1]
  вещ таб C[0:n-1,0:n-1]
  вещ sum, err, absv, res_l2_sq, res_l2, res_inf

Инициализация матриц

Заполняем $A$ и $B$ в одном двойном цикле. $A_{i,j} = 1$ при $j \le i$ (нижний треугольник), $B_{i,j} = 1$ при $j \ge i$ (верхний треугольник):

  | Инициализация A (нижний треугольник) и B (верхний треугольник)
  нц для i от 0 до n-1
    нц для j от 0 до n-1
      если j <= i то A[i,j] := 1.0 иначе A[i,j] := 0.0 все
      если j >= i то B[i,j] := 1.0 иначе B[i,j] := 0.0 все
    кц
  кц

Тройной цикл умножения

Классический алгоритм: для каждой пары $(i, j)$ вычисляем скалярное произведение $i$-й строки $A$ на $j$-й столбец $B$:

  | Вычисляем C = A * B
  нц для i от 0 до n-1
    нц для j от 0 до n-1
      sum := 0.0
      нц для k от 0 до n-1
        sum := sum + A[i,k] * B[k,j]
      кц
      C[i,j] := sum
    кц
  кц

Проверка результата

Аналитически $C_{i,j} = \min(i, j) + 1$. В коде это выражено условием: если $i < j$, ожидаемое значение равно $i + 1$, иначе $j + 1$. Вычисляем ошибку и её нормы:

  | Проверка: ожидаем C[i,j] = max(0, i - j + 1)
  res_l2_sq := 0.0
  res_inf := 0.0
  нц для i от 0 до n-1
    нц для j от 0 до n-1
      если i < j то
        err := C[i,j] - (i + 1.0)
      иначе
        err := C[i,j] - (j + 1.0)
      все
      absv := err; если absv < 0.0 то absv := -absv все
      res_l2_sq := res_l2_sq + err * err
      если absv > res_inf то res_inf := absv все
    кц
  кц
  res_l2 := sqrt(res_l2_sq)

  вывод "matmat: ||e||_2^2 = ", res_l2_sq, нс
  вывод "matmat: ||e||_2    = ", res_l2, нс
  вывод "matmat: ||e||_inf  = ", res_inf, нс

Сложность наивного умножения матриц -- $O(n^3)$ (три вложенных цикла). Существуют асимптотически более быстрые алгоритмы (например, алгоритм Штрассена со сложностью $O(n^{2.81})$), но на практике для небольших матриц тройной цикл работает хорошо.

Полная программа

алг main
нач
  | Матматр: C = A * B. A — нижнетреугольная матрица из единиц,
  | B — верхнетреугольная матрица из единиц. Тогда
  | C[i,j] = число k, таких что j <= k <= i = max(0, i - j + 1).

  цел n,i,j,k
  n := 128
  вещ таб A[0:n-1,0:n-1]
  вещ таб B[0:n-1,0:n-1]
  вещ таб C[0:n-1,0:n-1]
  вещ sum, err, absv, res_l2_sq, res_l2, res_inf

  | Инициализация A (нижний треугольник) и B (верхний треугольник)
  нц для i от 0 до n-1
    нц для j от 0 до n-1
      если j <= i то A[i,j] := 1.0 иначе A[i,j] := 0.0 все
      если j >= i то B[i,j] := 1.0 иначе B[i,j] := 0.0 все
    кц
  кц

  | Вычисляем C = A * B
  нц для i от 0 до n-1
    нц для j от 0 до n-1
      sum := 0.0
      нц для k от 0 до n-1
        sum := sum + A[i,k] * B[k,j]
      кц
      C[i,j] := sum
    кц
  кц

  | Проверка: ожидаем C[i,j] = max(0, i - j + 1)
  res_l2_sq := 0.0
  res_inf := 0.0
  нц для i от 0 до n-1
    нц для j от 0 до n-1
      если i < j то
        err := C[i,j] - (i + 1.0)
      иначе
        err := C[i,j] - (j + 1.0)
      все
      absv := err; если absv < 0.0 то absv := -absv все
      res_l2_sq := res_l2_sq + err * err
      если absv > res_inf то res_inf := absv все
    кц
  кц
  res_l2 := sqrt(res_l2_sq)

  вывод "matmat: ||e||_2^2 = ", res_l2_sq, нс
  вывод "matmat: ||e||_2    = ", res_l2, нс
  вывод "matmat: ||e||_inf  = ", res_inf, нс
кон

▶ Запустить пример