矩阵乘法的分治算法实现思路

Question

为了提高编程能力和算法，一边看算法导论，一边实现。

到矩阵乘法分治求解的时候，思路很简单，我也理解了，可是怎么实现伪代码时出现了问题，搞了一天也没搞出来。

大致意思就是把矩阵切田字块，求解子矩阵，再合并。

导论提示用下标。我就试着按它来，自己写了个矩阵结构体二位数组，并且为了整洁地定位一个矩阵的二维下标，又写了个位置结构体。

可是递归函数的调用和传参就是弄不明白
我现在想到的递归函数参数有 a,b,c三个矩阵的引用(因为算完值要放进c中)和三个矩阵要操作的部分(子矩阵)位置，一开始觉得没什么问题，可是写到两个递归函数调用再相加赋值给C矩阵，我就懵了，我不能把c传进去，因为求和才是c那个元素的值，可是递归出口(子矩阵row为1时)又的确赋值给c了。。。

智商不够的我彻底懵了，思维乱掉了，想了一天，画了一天图也没搞出来，求前辈指点。

算法的实现太繁琐，不用麻烦，只要前辈，大神们给个思路就好，先谢谢了

Answer 1

转载自这里，稍有改动。与原书中的伪代码思路是相同的。

// C++
#include <iostream>
#include <vector>

template<typename T>
struct Matrix {
  Matrix(size_t r, size_t c) {
    Data.resize(c, std::vector<T>(r, 0));
  }

  void SetSubMatrix(const int r, const int c, const int rn, const int cn,
                    const Matrix<T>& A, const Matrix<T>& B) {
    for (int cl = c; cl < cn; ++cl)
      for (int rl = r; rl < rn; ++rl)
        Data[cl][rl] = A.Data[cl - c][rl - r] + B.Data[cl - c][rl - r];
  }

  static Matrix<T> SquareMultiplyRecursive(Matrix<T>& A, Matrix<T>& B,
                                           int ar, int ac, int br, int bc, int n) {
    Matrix<T> C(n, n);

    if (n == 1) {
      C.Data[0][0] = A.Data[ac][ar] * B.Data[bc][br];
    } else {
      C.SetSubMatrix(0, 0, n / 2, n / 2,
        SquareMultiplyRecursive(A, B, ar, ac, br, bc, n / 2),
        SquareMultiplyRecursive(A, B, ar, ac + (n / 2), br + (n / 2), bc, n / 2));

      C.SetSubMatrix(0, n / 2, n / 2, n,
        SquareMultiplyRecursive(A, B, ar, ac, br, bc + (n / 2), n / 2),
        SquareMultiplyRecursive(A, B, ar, ac + (n / 2), br + (n / 2), bc + (n / 2), n / 2));

      C.SetSubMatrix(n / 2, 0, n, n / 2,
        SquareMultiplyRecursive(A, B, ar + (n / 2), ac, br, bc, n / 2),
        SquareMultiplyRecursive(A, B, ar + (n / 2), ac + (n / 2), br + (n / 2), bc, n / 2));

      C.SetSubMatrix(n / 2, n / 2, n, n,
        SquareMultiplyRecursive(A, B, ar + (n / 2), ac, br, bc + (n / 2), n / 2),
        SquareMultiplyRecursive(A, B, ar + (n / 2), ac + (n / 2), br + (n / 2), bc + (n / 2), n / 2));
    }

    return C;
  }

  void Print() {
    for (size_t c = 0; c < Data.size(); ++c) {
      for (size_t r = 0; r < Data[0].size(); ++r)
        std::cout << Data[c][r] << " ";
      std::cout << "\n";
    }
    std::cout << "\n";
  }

  std::vector<std::vector<T> > Data;
};

int main() {
  Matrix<int> A(2, 2);
  Matrix<int> B(2, 2);

  A.Data[0][0] = 2; A.Data[0][1] = 1;
  A.Data[1][0] = 1; A.Data[1][1] = 2;

  B.Data[0][0] = 2; B.Data[0][1] = 1;
  B.Data[1][0] = 1; B.Data[1][1] = 2;

  A.Print();
  B.Print();

  Matrix<int> C(Matrix<int>::SquareMultiplyRecursive(A, B, 0, 0, 0, 0, 2));

  C.Print();
}