算法导论学习——分治矩阵乘法
BDplanDante 2016-04-09
算法导论学习——分治矩阵乘法
头文件 结构的定义
stdafx.h
// stdafx.h : 标准系统包含文件的包含文件，
// 或是经常使用但不常更改的
// 特定于项目的包含文件
//
#pragma once
#include "targetver.h"
#include <stdio.h>
#include <tchar.h>
// TODO: 在此处引用程序需要的其他头文件
#include <iostream>
using namespace std;
#define Maxelem 10
//求两个数的最大值
int inline max(int a, int b){
return a >= b ? a : b;
}
//求三个数的最大值
int inline max(int a, int b, int c){
return max(max(a, b), c);
}
//求min(2^x),ST. 2^x>=number。
int inline L2n(int number){
if ((number & number - 1) == 0)
{
return number;
}
else {
return pow(2, (int)log2(number) + 1);
}
}
//定义矩阵的结构
class Matrix{
public:
int data[Maxelem][Maxelem];
int M, N;
//以数组复制的方式构造矩阵对象，其中起始位置为0.
Matrix(int array[Maxelem][Maxelem], int m, int n){
M = m;
N = n;
for (int i = 0; i < m; i++){
for (int j = 0; j < n; j++){
this->data[i][j] = array[i][j];
}
}
}
//以数组复制的方式构造矩阵对象，以strat1,end1,strat2,end2为区间
Matrix(int array[Maxelem][Maxelem], int start1, int end1, int start2, int end2){
M = end1 - start1 + 1;
N = end2 - start2 + 1;
int p = 0, q = 0;
for (int i = 0; i < M; i++){
for (int j = 0; j < N; j++){
data[i][j] = array[start1 + i][start2 +j];
}
}
}
//仅构造矩阵，不填充数据
Matrix(int m, int n) :M(m), N(n){}
/*
打印输出
*/
void print() const {
for (int i = 0; i <M; i++){
for (int j = 0; j < N; j++){
cout << data[i][j] << " ";
}
cout << endl;
}
}
//为矩阵补充0，使其成为标准方阵
void fill(){
if (!(M == N && ((M & M - 1) == 0))){
int n = L2n(max(M, N));
for (int i = 0; i < n; i++){
for (int j = 0; j < n; j++){
data[i][j] = (i < M&&j < N ? data[i][j] : 0);
}
}
M = n;
N = n;
}
}
/*重载二维运算符[][]*/
int * const operator[](const int i)
{
return data[i];
}
Matrix friend operator +(Matrix m1, Matrix m2){
Matrix op = Matrix(m1.M, m1.N);
for (int i = 0; i < m1.M; i++)
{
for (int j = 0; j < m1.M; j++){
op[i][j] = m1[i][j] + m2[i][j];
}
}
return op;
}
Matrix friend operator -(Matrix m1, Matrix m2){
Matrix op = Matrix(m1.M, m1.N);
for (int i = 0; i < m1.M; i++)
{
for (int j = 0; j < m1.M; j++){
op[i][j] = m1[i][j] - m2[i][j];
}
}
return op;
}
/*
将计算过程中补充的0清除。计算完毕后才能用的方法，不加也能得到结果，不过行列数不对。
*/
void clean(int m, int n){
M = m;
N = n;
}
};
算法的实现：
// strassenAlgorithm.cpp : 定义控制台应用程序的入口点。
//
#include "stdafx.h"
//该方法只能相乘最简单的2*2矩阵。
Matrix mutilSimple(Matrix A, Matrix B){
int a = A[0][0],b=A[0][1],c=A[1][0],d=A[1][1];
int e = B[0][0], f = B[0][1], g = B[1][0], h = B[1][1];
int p1 = a*(f - h);
int p2 = (a + b)*h;
int p3 = (c + d)*e;
int p4 = d*(g - e);
int p5 = (a + d)*(e + h);
int p6 = (b - d)*(g + h);
int p7 = (a - c)*(e + f);
Matrix returnValue = Matrix(2, 2);
returnValue[0][0] = p5 + p4 - p2 + p6;
returnValue[0][1] = p1 + p2;
returnValue[1][0] = p3 + p4;
returnValue[1][1] = p1 + p5 - p3 - p7;
return returnValue;
}
//矩阵乘法 必须用了fill方法才能相乘
Matrix mutilMerge(Matrix A,Matrix B){
if (A.M == 2 && B.M == 2){
return mutilSimple(A, B);
}
int k = A.M;
Matrix
a = Matrix(A.data, 0, k / 2 - 1, 0, k / 2 - 1),
b = Matrix(A.data, 0, k / 2 - 1, k / 2, k - 1),
c = Matrix(A.data, k / 2, k - 1, 0, k / 2 - 1),
d = Matrix(A.data, k / 2, k - 1, k / 2, k - 1),
e = Matrix(B.data, 0, k / 2 - 1, 0, k / 2 - 1),
f = Matrix(B.data, 0, k / 2 - 1, k / 2, k - 1),
g = Matrix(B.data, k / 2, k - 1, 0, k / 2 - 1),
h = Matrix(B.data, k / 2, k - 1, k / 2, k - 1),
op = Matrix(k, k),
p1 = mutilMerge(a, f - h),
p2 = mutilMerge(a + b, h),
p3 = mutilMerge(c + d, e),
p4 = mutilMerge(d, g - e),
p5 = mutilMerge(a + d, e + h),
p6 = mutilMerge(b - d, g + h),
p7 = mutilMerge(a - c, e + f),
op1=p5+p4-p2+p6,
op2=p1+p2,
op3=p3+p4,
op4 = p1 + p5 - p3 - p7;
int x1 = 0, y1 = 0, x2 = 0, y2 = 0, x3 = 0,y3=0,x4=0, y4 = 0; //4个变量的游标
int u = 0, v = 0;
for (int i = 0; i < k; i++)
{
for (int j = 0; j < k; j++){
if (i >= 0 && i <= k / 2 - 1 && j >= 0 && j <= k / 2 - 1){
op[i][j] = op1[x1][y1];
y1++;
if (y1 == op1.M) { y1 = 0; x1++; }
}
if (i >= 0 && i <= k / 2 - 1 && j >= k / 2 && j <= k - 1){
op[i][j] = op2[x2][y2];
y2++;
if (y2 == op2.M) { y2 = 0; x2++; }
}
if (i >= k/2 && i <= k - 1 && j >= 0 && j <= k / 2 - 1){
op[i][j] = op3[x3][y3];
y3++;
if (y3 == op3.M) { y3 = 0; x3++; }
}
if (i >= k / 2 && i <= k - 1 && j >= k/2 && j <= k - 1){
op[i][j] = op4[x4][y4];
y4++;
if (y4 == op4.M) { y4 = 0; x4++; }
}
}
}
return op;
}
Matrix zeroclear( Matrix result,int M,int N){
result.M = M;
result.N = N;
return result;
}
int _tmain(int argc, _TCHAR* argv[])
{
int matrixA[Maxelem][Maxelem] = {
{ 10, 3, 3 ,7,4},
{ 5, 3, 8,2,1 },
{ -2, 3, 7, 5, 2 },
{1,10,-2,1,8},
{3,3,3,3,3}
};
int matrixB[Maxelem][Maxelem] = {
{ -4, 6, 1,2,1 },
{ 9, 10, 8,0,3 },
{ 2, 3, -7 ,-1,-1},
{ 1, -6, 2, 1, 7 },
{1,2,3,4,5}
};
Matrix ma = Matrix(matrixA,5,5);
Matrix mb = Matrix(matrixB,5,5);
cout << "A="<<endl;
ma.print();
cout << "B="<<endl;
mb.print();
ma.fill();
mb.fill();
Matrix result = mutilMerge(ma, mb);
cout << "A×B=" << endl;
result.clean(5, 5);
result.print();
cout << "B×A=" << endl;
result = mutilMerge(mb, ma);
result.clean(5, 5);
result.print();
system("pause");
return 0;
}
