linux多线程编写矩阵乘法,算法导论学习——分治矩阵乘法

算法导论学习——分治矩阵乘法

头文件 结构的定义

stdafx.h

// stdafx.h : 标准系统包含文件的包含文件，

// 或是经常使用但不常更改的

// 特定于项目的包含文件

//

#pragma once

#include "targetver.h"

#include

#include

// TODO:  在此处引用程序需要的其他头文件

#include

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

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="<

ma.print();

cout << "B="<

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;

}