← Все примеры

Умножение матрицы на вектор

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

Задача

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

$y_i = \sum_{j=0}^{n-1} A_{i,j} \cdot x_j, \quad i = 0, \ldots, n-1$

В примере используется нижнетреугольная матрица из единиц:

$A_{i,j} = \begin{cases} 1, & j \le i \ 0, & j > i \end{cases}$

и единичный вектор $x_j = 1$. Ожидаемый результат: $y_i = i + 1$.

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

1. Заполнить матрицу $A$ (нижний треугольник -- единицы, остальное -- нули) и вектор $x$ (все элементы равны $1$).
2. Для каждой строки $i$ вычислить скалярное произведение строки $A_i$ на вектор $x$, записать результат в $y_i$.
3. Проверить результат: вычислить ошибку $e_i = y_i - (i + 1)$ и вывести нормы $|e|2$ и $|e|\infty$.

Разбор

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

Задаём размер $n = 256$ и объявляем двумерный массив A для матрицы, одномерные массивы x и y для входного и выходного векторов:

  цел n,i,j
  n := 256
  вещ таб A[0:n-1,0:n-1]
  вещ таб x[0:n-1]
  вещ таб y[0:n-1]
  вещ sum, absv
  вещ res_l2_sq, res_l2, res_inf, err

Инициализация матрицы и вектора

Заполняем $A$ как нижнетреугольную матрицу из единиц ($A_{i,j} = 1$ при $j \le i$, иначе $0$), а вектор $x$ -- единицами:

  | Инициализация A и x
  нц для i от 0 до n-1
    нц для j от 0 до n-1
      если j <= i то
        A[i,j] := 1.0
      иначе
        A[i,j] := 0.0
      все
    кц
    x[i] := 1.0
  кц

Умножение матрицы на вектор

Основной вычислительный блок -- два вложенных цикла. Для каждой строки $i$ считаем сумму $y_i = \sum_j A_{i,j} \cdot x_j$:

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

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

Для каждого $i$ вычисляем ошибку $e_i = y_i - (i+1)$ и накапливаем нормы -- $|e|2$ (евклидова норма, корень из суммы квадратов) и $|e|\infty$ (максимальное отклонение):

  | Проверка: ожидаем y[i] = i+1
  res_l2_sq := 0.0
  res_inf := 0.0
  нц для i от 0 до n-1
    err := y[i] - (i + 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)

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

Сложность умножения матрицы на вектор -- $O(n^2)$. Это фундаментальная операция, используемая как строительный блок во всех итерационных методах решения систем линейных алгебраических уравнений (СЛАУ) -- набора уравнений вида $Ax = b$.

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

алг main
нач
  | Матвект: y = A * x, где A — нижнетреугольная матрица из единиц,
  | x — вектор из единиц. Ожидается y[i] = i+1.

  цел n,i,j
  n := 256
  вещ таб A[0:n-1,0:n-1]
  вещ таб x[0:n-1]
  вещ таб y[0:n-1]
  вещ sum, absv
  вещ res_l2_sq, res_l2, res_inf, err

  | Инициализация A и x
  нц для i от 0 до n-1
    нц для j от 0 до n-1
      если j <= i то
        A[i,j] := 1.0
      иначе
        A[i,j] := 0.0
      все
    кц
    x[i] := 1.0
  кц

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

  | Проверка: ожидаем y[i] = i+1
  res_l2_sq := 0.0
  res_inf := 0.0
  нц для i от 0 до n-1
    err := y[i] - (i + 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)

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

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