算法导论 矩阵相乘(Strassen方法)

#include <stdio.h>
#include <stdlib.h>
#define LEN 4

typedef struct
{
int** S;
int n;
}Square;

void printSquare(Square A)
{
printf("length is %d\n",A.n);
for(int i=0;i<A.n;i++)
{
for(int j=0;j<A.n;j++)
{
printf("%4d",A.S[i][j]);
}
printf("\n");
}
}

void freeSquare(Square A)
{
for(int i=0;i<A.n;i++)
{
free(*(A.S+i));
}
free(A.S);
}

Square makeSquare(int n)
{
int** S=(int **)malloc(sizeof(int *)*n);
for(int i=0;i<n;i++)
{
*(S+i)=(int*)malloc(sizeof(int)*n);
}
for(int i=0;i<n;i++)
{
for(int j=0;j<n;j++)
{
S[i][j]=0;
}
}
Square sq={S,n};
return sq;
}

Square divideSquare(Square A,int rowstart,int rowend,int colstart,int colend)
{
int n=rowend-rowstart+1;
Square sq=makeSquare(n);
for(int i=rowstart,ii=0;i<=rowend;i++,ii++)
{
for(int j=colstart,jj=0;j<=colend;j++,jj++)
{
sq.S[ii][jj]=A.S[i][j];
}
}
return sq;
}

Square squareCalc(Square A,Square B,char op)
{
int n=A.n;
Square sq=makeSquare(n);
for(int i=0;i<n;i++)
{
for(int j=0;j<n;j++)
{
if(op=='+')
sq.S[i][j]=A.S[i][j]+B.S[i][j];
else if(op=='-')
sq.S[i][j]=A.S[i][j]-B.S[i][j];
}
}
return sq;
}

void StrassenMul(Square C,Square A,Square B)
{
if(A.n==1)
{
C.S[0][0]=A.S[0][0]*B.S[0][0];
return;
}

Square A11=divideSquare(A,0,A.n/2-1,0,A.n/2-1);
Square A12=divideSquare(A,0,A.n/2-1,A.n/2,A.n-1);
Square A21=divideSquare(A,A.n/2,A.n-1,0,A.n/2-1);
Square A22=divideSquare(A,A.n/2,A.n-1,A.n/2,A.n-1);
Square B11=divideSquare(B,0,B.n/2-1,0,B.n/2-1);
Square B12=divideSquare(B,0,B.n/2-1,B.n/2,B.n-1);
Square B21=divideSquare(B,B.n/2,B.n-1,0,B.n/2-1);
Square B22=divideSquare(B,B.n/2,B.n-1,B.n/2,B.n-1);

Square S1=squareCalc(B12,B22,'-');
Square S2=squareCalc(A11,A12,'+');
Square S3=squareCalc(A21,A22,'+');
Square S4=squareCalc(B21,B11,'-');
Square S5=squareCalc(A11,A22,'+');
Square S6=squareCalc(B11,B22,'+');
Square S7=squareCalc(A12,A22,'-');
Square S8=squareCalc(B21,B22,'+');
Square S9=squareCalc(A11,A21,'-');
Square S10=squareCalc(B11,B12,'+');

int n=A11.n;
Square P1=makeSquare(n);
Square P2=makeSquare(n);
Square P3=makeSquare(n);
Square P4=makeSquare(n);
Square P5=makeSquare(n);
Square P6=makeSquare(n);
Square P7=makeSquare(n);

StrassenMul(P1,A11,S1);
StrassenMul(P2,S2,B22);
StrassenMul(P3,S3,B11);
StrassenMul(P4,A22,S4);
StrassenMul(P5,S5,S6);
StrassenMul(P6,S7,S8);
StrassenMul(P7,S9,S10);

Square C111=squareCalc(P5,P4,'+');
Square C112=squareCalc(C111,P2,'-');
Square C11=squareCalc(C112,P6,'+');
Square C12=squareCalc(P1,P2,'+');
Square C21=squareCalc(P3,P4,'+');
Square C221=squareCalc(P5,P1,'+');
Square C222=squareCalc(C221,P3,'-');
Square C22=squareCalc(C222,P7,'-');

int i,j,ii,jj;
for(i=0,ii=0;i<=C.n/2-1;i++,ii++)
{
for(j=0,jj=0;j<=C.n/2-1;j++,jj++)
{
C.S[i][j]=C11.S[ii][jj];
}
}
for(i=0,ii=0;i<=C.n/2-1;i++,ii++)
{
for(j=C.n/2,jj=0;j<=C.n-1;j++,jj++)
{
C.S[i][j]=C12.S[ii][jj];
}
}
for(i=C.n/2,ii=0;i<=C.n-1;i++,ii++)
{
for(j=0,jj=0;j<=C.n/2-1;j++,jj++)
{
C.S[i][j]=C21.S[ii][jj];
}
}
for(i=C.n/2,ii=0;i<=C.n-1;i++,ii++)
{
for(j=C.n/2,jj=0;j<=C.n-1;j++,jj++)
{
C.S[i][j]=C22.S[ii][jj];
}
}
freeSquare(A11);
freeSquare(A12);
freeSquare(A21);
freeSquare(A22);
freeSquare(B11);
freeSquare(B12);
freeSquare(B21);
freeSquare(B22);

freeSquare(S1);
freeSquare(S2);
freeSquare(S3);
freeSquare(S4);
freeSquare(S5);
freeSquare(S6);
freeSquare(S7);
freeSquare(S8);
freeSquare(S9);
freeSquare(S10);

freeSquare(P1);
freeSquare(P2);
freeSquare(P3);
freeSquare(P4);
freeSquare(P5);
freeSquare(P6);
freeSquare(P7);

freeSquare(C111);
freeSquare(C112);
freeSquare(C11);
freeSquare(C12);
freeSquare(C21);
freeSquare(C221);
freeSquare(C222);
freeSquare(C22);
}

int main()
{
int AS[LEN][LEN]={{1,2,3,4},{3,4,5,6},{5,6,7,8},{7,8,9,10}};
int BS[LEN][LEN]={{5,6,7,8},{7,8,9,10},{11,12,13,14},{15,16,17,18}};
Square A=makeSquare(LEN);
for(int i=0;i<LEN;i++)
{
for(int j=0;j<LEN;j++)
{
A.S[i][j]=AS[i][j];
}
}
Square B=makeSquare(LEN);
for(int i=0;i<LEN;i++)
{
for(int j=0;j<LEN;j++)
{
B.S[i][j]=BS[i][j];
}
}
Square C=makeSquare(LEN);
StrassenMul(C,A,B);
printSquare(C);
getchar();
}