The Program of Algorithms ------- Diveide and Conquer ---- Matrix Multiplication

Q:

Input:     A=[a_ij] B=[b_ij]    i,d=1,2,3...n

Output:  C=[c_ij]=A*B

 

Naive Algorithm:  Θ(n^3)

Realization --- CPP:

long int** NaiveMatrixMulitiple(long int** lhs,long int** rhs,long int size){
	long int ** c=new long int*[size];
	for(int i=0;i<size;i++){
		c[i]=new long int[size];
	}
	
	for(int i=0;i<size;i++){
		for(int j=0;j<size;j++){
			c[i][j]=0;
			for(int k=0;k<size;k++){
				c[i][j]+=lhs[i][k]*rhs[k][j];
			}
		}
	}
	return c;
}

Using Divide-and-Conquer Algorithm, the time cost is also Θ(n^3).

Recursively equation as follow:

   T(n)=8T(n/2)+Θ(n^2)    ----a=8 b=2 n^logb(a)=n^3  So, the T(n)=Θ(n^3). Using Master Method Case1.

 

Improve: Strassen`s algorithm.

   Idea: reduce the Multiplication from 8 to 7.

         P1=.....P7=....

         r=P5+P4-P2+P6...s=....t=.....u=....

   Divide and Conquer Process:

         1. Divide:       Divide A and B, compute terms for P(compute the a,b,c...f).

         2. Conquer:   recursively compute P1,....P7

         3. Combine:   combine the P1...P7 to form the r,s,t,u(Result Matrix)

   Time Costs: T(n)=7T(n/2)+Θ(n^2)=Θ(n^(lg7))=Θ(n^2.81)

   Realization --- CPP:

long int** MatrixAdd(long int** lhs,long int** rhs,long int size){
	long int **Res=new long int*[size];
	for(int i=0;i<size;i++){
		Res[i]=new long int[size];
	}
	for(int i=0;i<size;i++){
		for(int j=0;j<size;j++){
			Res[i][j]=lhs[i][j]+rhs[i][j];
		}
	}
	return Res;
}
long int** MatrixSub(long int** lhs,long int** rhs,int size){
	long int **Res=new long int*[size];
	for(int i=0;i<size;i++){
		Res[i]=new long int[size];
	}
	for(int i=0;i<size;i++){
		for(int j=0;j<size;j++){
			Res[i][j]=lhs[i][j]-rhs[i][j];
		}
	}
	return Res;
}

long int** MatrixMulitiple(long int** lhs,long int** rhs,int size){

	long int **Res=new long int*[size];
	for(int i=0;i<size;i++){
		Res[i]=new long int[size];
	}
	if(size==2){
		Res[0][0]=lhs[0][0]*rhs[0][0]+lhs[0][1]*rhs[1][0];
		Res[0][1]=lhs[0][0]*rhs[0][1]+lhs[0][1]*rhs[1][1];
		Res[1][0]=lhs[1][0]*rhs[0][0]+lhs[1][1]*rhs[1][0];
		Res[1][1]=lhs[1][0]*rhs[0][1]+lhs[1][1]*rhs[1][1];
		return Res;
	}

	int smallSize=size/2;
	long int **a=lhs;
	long int **b=new long int*[smallSize];
	for(int i=0;i<smallSize;i++){
		b[i]=new long int[smallSize];
	}
	for(int i=0;i<smallSize;i++){
		for(int j=0;j<smallSize;j++)
			b[i][j]=lhs[i][smallSize+j];
	}
	long int **c=&(lhs[smallSize]);
	long int **d=new long int*[smallSize];
	for(int i=0;i<smallSize;i++){
		d[i]=new long int[smallSize];
	}
	for(int i=0;i<smallSize;i++){
		for(int j=0;j<smallSize;j++)
			d[i][j]=lhs[smallSize+i][smallSize+j];
	}


	long int **e=rhs;
	long int **f=new long int*[smallSize];
	for(int i=0;i<smallSize;i++){
		f[i]=new long int[smallSize];
	}
	for(int i=0;i<smallSize;i++){
		for(int j=0;j<smallSize;j++)
			f[i][j]=rhs[i][smallSize+j];
	}
	long int **g=&(rhs[smallSize]);
	long int **h=new long int*[smallSize];
	for(int i=0;i<smallSize;i++){
		h[i]=new long int[smallSize];
	}
	for(int i=0;i<smallSize;i++){
		for(int j=0;j<smallSize;j++)
			h[i][j]=rhs[smallSize+i][smallSize+j];
	}

	long int **P1=MatrixMulitiple(a,MatrixSub(f,h,smallSize),smallSize);
	long int **P2=MatrixMulitiple(MatrixAdd(a,b,smallSize),h,smallSize);
	long int **P3=MatrixMulitiple(MatrixAdd(c,d,smallSize),e,smallSize);
	long int **P4=MatrixMulitiple(d,MatrixSub(g,e,smallSize),smallSize);
	long int **P5=MatrixMulitiple(MatrixAdd(a,d,smallSize),MatrixAdd(e,h,smallSize),smallSize);
	long int **P6=MatrixMulitiple(MatrixSub(b,d,smallSize),MatrixAdd(g,h,smallSize),smallSize);
	long int **P7=MatrixMulitiple(MatrixSub(a,c,smallSize),MatrixAdd(e,f,smallSize),smallSize);

	long int **r=MatrixAdd(MatrixSub(MatrixAdd(P5,P4,smallSize),P2,smallSize),P6,smallSize);
	long int **s=MatrixAdd(P1,P2,smallSize);
	long int **t=MatrixAdd(P3,P4,smallSize);
	long int **u=MatrixSub(MatrixSub(MatrixAdd(P5,P1,smallSize),P3,smallSize),P7,smallSize);
	
	for(int i=0;i<size;i++){
		for(int j=0;j<size;j++){
			if(i<smallSize&&j<smallSize){
				Res[i][j]=r[i][j];
			}else if(i<smallSize&&j>=smallSize){
				Res[i][j]=s[i][j-smallSize];
			}else if(i>=smallSize&&j<smallSize){
				Res[i][j]=t[i-smallSize][j];
			}else if(i>=smallSize&&j>=smallSize){
				Res[i][j]=u[i-smallSize][j-smallSize];
			}
		}
	}
	return Res;
}