C++ 方矩阵乘法 + Strassen矩阵

  这几天看算法导论，看到矩阵一章，就实现了一下。

下面是普通的矩阵乘法，复杂度为：n^3。

template<unsigned M,unsigned N, unsigned Q>
void Square_matrix_multiply(int(&A)[M][N], int(&B)[N][Q], int(&C)[M][Q]) {                 
	for (size_t i = 0;i != M;++i) {
		for (size_t j = 0;j != Q;++j) {
			C[i][j] = 0;
			for (size_t n = 0;n != N;++n) {
				C[i][j] += A[i][n] * B[n][j];
			}
		}
	}
}

函数接受三个二维数组，A * B得到的矩阵赋值给C。

下面是分治策略的算法。

template<typename T>
Matrix Square_max_matrix_multiply_recursive(const T &A, const T &B) {
	size_t n = A.rows();
	Matrix C(n, n);
	if (n == 1)
		return C = A.get()*B.get();
	else {
		MatrixRef A_11(A, 0, 0), A_12(A, 0, n / 2), A_21(A, n / 2, 0), A_22(A, n / 2, n / 2);    // 使用一个类MatirxRef
		MatrixRef B_11(B, 0, 0), B_12(B, 0, n / 2), B_21(B, n / 2, 0), B_22(B, n / 2, n / 2);    // 含有三个size_t类型。其中两个实现坐标，一个指明矩阵长度
		MatrixRef C_11(C, 0, 0), C_12(C, 0, n / 2), C_21(C, n / 2, 0), C_22(C, n / 2, n / 2);    // 进行分割
		C_11 = Square_max_matrix_multiply_recursive(A_11, B_11) + Square_max_matrix_multiply_recursive(A_12, B_21);  // Matrix::operator+;
		C_12 = Square_max_matrix_multiply_recursive(A_11, B_12) + Square_max_matrix_multiply_recursive(A_12, B_22);  // MatrixRef::operator=;
		C_21 = Square_max_matrix_multiply_recursive(A_21, B_11) + Square_max_matrix_multiply_recursive(A_22, B_21);
		C_22 = Square_max_matrix_multiply_recursive(A_21, B_12) + Square_max_matrix_multiply_recursive(A_22, B_22);
	}
	return C;
}

矩阵实现了一个Matrix类（具体实现在最下面），有一个构造函数：接受两个size_t值l、r，生成l*r大小值全为0的矩阵。

Matrix::Matrix(size_t l, size_t r) : hight(l), width(r), data(make_shared<vector<int>>()) {
	data->resize(l*r);
}

其中hight为矩阵行高，width为列宽，data为shared_ptr，矩阵用vector实现。

A.rows()返回A的width长度（即方矩阵的边长），Matrix（n，n）创建一个矩阵。

size_t rows() const {
		return width;
	}

如果n==1，通过Matrix的get函数返回第一个元素，也就是唯一的一个元素。

int Matrix::get() const {
	return (*data)[0];
}

为了不复制矩阵元素（如果可以复制矩阵元素的话，会简单很多），另实现了一个MatrixRef，其含有：两个size_t数据成员（实现坐标点）、一个size_t数据成员（实现矩阵长度）、一个weak_ptr（指向vector<int>)。

MatrixRef含有两个构造函数：一个接受Matrix加两个size_t；一个接受MatrixRef加两个size_t。都是为了指明引用范围。

MatrixRef::MatrixRef(const Matrix &m, size_t line, size_t row) : wptr(m.data), 
      hight_startptr(line), width_startptr(row), length(m.rows() / 2) { }
MatrixRef::MatrixRef(const MatrixRef &mref, size_t line, size_t row) : wptr(mref.wptr),                 
      hight_startptr(mref.hight_startptr + line),
      width_startptr(mref.width_startptr + row), length(mref.rows() / 2) { }

wptr用data或wptr初始化，避免拷贝。length为rows（）的返回值除以2，因为是分割为4个矩阵，行列各除以2。

要注意：接受MatrixRef的坐标要加上之前的坐标。

MatrixRef也有一个rows成员函数，为了递归调用。

size_t rows() const {
		return length;
	}

Square_max_matrix_multiply_recursive函数返回一个Matrix，Matrix实现了operator+，但是行列必须相等。

Matrix& Matrix::operator+=(const Matrix &rhs) {
	if (hight == rhs.hight && width == rhs.width) {
		for (size_t i = 0;i != size();++i)
			(*data)[i] += (*rhs.data)[i];
	}
	else
		throw std::logic_error("Not Matched");
	return *this;
}
Matrix operator+(const Matrix &lhs, const Matrix &rhs) {
	Matrix m(lhs);
	return m += rhs;
}

MatrixRef实现了一个operato