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

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

Проблема

Требуется эффективно выполнить умножение двух матриц.

Решение

Пример 11.32 показывает, как можно выполнить умножение матриц, причем эта реализация подходит как для динамических, так и для статических матриц. Формально этот алгоритм реализует соотношение A=A+B*C, которое (возможно, неожиданно) вычисляется более эффективно, чем A=B*C.

Пример 11.32. Умножение матриц

#include "matrix.hpp" // рецепт 11.13

#include "kmatrix.hpp" // рецепт 11.14

#include <iostream>

#include <cassert>

using namespace std;

template<class M1, class M2, class M3>

void matrixMultiply(const M1& m1, const M2& m2, M3& m3) {

 assert(m1.cols() == m2.rows());

 assert(m1.rows() == m3.rows());

 assert(m2.cols() == m3.cols());

 for (int i=m1.rows()-1; i >= 0; --i) {

  for (int j=m2.cols()-1; j >= 0; --j) {

   for (int k = m1.cols()-1; k >= 0; --k) {

    m3[i][j] += m1[i][k] * m2[k][j];

   }

  }

 }

}

int main() {

 matrix<int> m1(2, 1);

 matrix<int> m2(1, 2);

 kmatrix<int, 2, 2> m3;

 m3 = 0;

 m1[0][0] = 1;

 m1[1][0] = 2;

 m2[0][0] = 3;

 m2[0][1] = 4;

 matrixMultlply(m1, m2, m3);

 cout << "(" << m3[0][0] << ", " << m3[0][1] << ")" << endl;

 cout << "(" << m3[1][0] << ", " << m3[1][1 ] << ")" << endl;

}

Программа примера 11.32 выдает следующий результат.

(3, 4)

(6, 8)

Обсуждение

При умножении двух матриц число столбцов первой матрицы должно равняться числу строк второй матрицы. Число строк полученной матрицы равно числу строк первой матрицы, а число столбцов равно числу столбцов второй матрицы. Я обеспечиваю эти условия в отладочной версии с помощью макроса assert, определенного в заголовочном файле <cassert>.

Решающее значение для эффективной реализации умножения имеет отсутствие избыточных операций по созданию и копированию временных объектов. Так, представленная в примере 11.32 функция умножения матриц передает результат по ссылке. Если бы алгоритм умножения я реализовал впрямую путем перегрузки оператора operator*, это привело бы к лишним операциям распределения, копирования и освобождения памяти, занимаемой временной матрицей. Потенциально такой подход может оказаться очень затратным при работе с большими матрицами.

В примере 11.32 реализуется равенство A=A+B*C, а не A=B*C, для того чтобы избежать лишней инициализации значений матрицы A.

Смотри также

Рецепт 11.17.