What is the source code for Strassen's matrix multiplication program?

Answer 1

#include
#include
#include
#include "2DArray.h"

#define GRAIN 1024 /* product size below which matmultleaf is used */

void seqMatMult(int m, int n, int p, double** A, double** B, double** C)
{
for (int i = 0; i < m; i++)
for (int j = 0; j < n; j++)
{
C[i][j] = 0.0;
for (int k = 0; k < p; k++)
C[i][j] += A[i][k]*B[k][j];
}
}


void matmultleaf(int mf, int ml, int nf, int nl, int pf, int pl, double **A, double **B, double **C)
/*
subroutine that uses the simple triple loop to multiply
a submatrix from A with a submatrix from B and store the
result in a submatrix of C.
*/
// mf, ml; /* first and last+1 i index */
// nf, nl; /* first and last+1 j index */
// pf, pl; /* first and last+1 k index */
{
for (int i = mf; i < ml; i++)
for (int j = nf; j < nl; j++)
for (int k = pf; k < pl; k++)
C[i][j] += A[i][k]*B[k][j];
}


void copyQtrMatrix(double **X, int m, double **Y, int mf, int nf)
{
for (int i = 0; i < m; i++)
X[i] = &Y[mf+i][nf];
}

void AddMatBlocks(double **T, int m, int n, double **X, double **Y)
{
for (int i = 0; i < m; i++)
for (int j = 0; j < n; j++)
T[i][j] = X[i][j] + Y[i][j];
}

void SubMatBlocks(double **T, int m, int n, double **X, double **Y)
{
for (int i = 0; i < m; i++)
for (int j = 0; j < n; j++)
T[i][j] = X[i][j] - Y[i][j];
}


void strassenMMult(int mf, int ml, int nf, int nl, int pf, int pl, double **A, double **B, double **C)
{
if ((ml-mf)*(nl-nf)*(pl-pf) < GRAIN)
matmultleaf(mf, ml, nf, nl, pf, pl, A, B, C);

else {
int m2 = (ml-mf)/2;
int n2 = (nl-nf)/2;
int p2 = (pl-pf)/2;

double **M1 = Allocate2DArray< double >(m2, n2);
double **M2 = Allocate2DArray< double >(m2, n2);
double **M3 = Allocate2DArray< double >(m2, n2);
double **M4 = Allocate2DArray< double >(m2, n2);
double **M5 = Allocate2DArray< double >(m2, n2);
double **M6 = Allocate2DArray< double >(m2, n2);
double **M7 = Allocate2DArray< double >(m2, n2);

double **A11 = new double*[m2];
double **A12 = new double*[m2];
double **A21 = new double*[m2];
double **A22 = new double*[m2];

double **B11 = new double*[p2];
double **B12 = new double*[p2];
double **B21 = new double*[p2];
double **B22 = new double*[p2];

double **C11 = new double*[m2];
double **C12 = new double*[m2];
double **C21 = new double*[m2];
double **C22 = new double*[m2];

double **tAM1 = Allocate2DArray< double >(m2, p2);
double **tBM1 = Allocate2DArray< double >(p2, n2);
double **tAM2 = Allocate2DArray< double >(m2, p2);
double **tBM3 = Allocate2DArray< double >(p2, n2);
double **tBM4 = Allocate2DArray< double >(p2, n2);
double **tAM5 = Allocate2DArray< double >(m2, p2);
double **tAM6 = Allocate2DArray< double >(m2, p2);
double **tBM6 = Allocate2DArray< double >(p2, n2);
double **tAM7 = Allocate2DArray< double >(m2, p2);
double **tBM7 = Allocate2DArray< double >(p2, n2);

copyQtrMatrix(A11, m2, A, mf, pf);
copyQtrMatrix(A12, m2, A, mf, p2);
copyQtrMatrix(A21, m2, A, m2, pf);
copyQtrMatrix(A22, m2, A, m2, p2);

copyQtrMatrix(B11, p2, B, pf, nf);
copyQtrMatrix(B12, p2, B, pf, n2);
copyQtrMatrix(B21, p2, B, p2, nf);
copyQtrMatrix(B22, p2, B, p2, n2);

copyQtrMatrix(C11, m2, C, mf, nf);
copyQtrMatrix(C12, m2, C, mf, n2);
copyQtrMatrix(C21, m2, C, m2, nf);
copyQtrMatrix(C22, m2, C, m2, n2);

// M1 = (A11 + A22)*(B11 + B22)
AddMatBlocks(tAM1, m2, p2, A11, A22);
AddMatBlocks(tBM1, p2, n2, B11, B22);
strassenMMult(0, m2, 0, n2, 0, p2, tAM1, tBM1, M1);

//M2 = (A21 + A22)*B11
AddMatBlocks(tAM2, m2, p2, A21, A22);
strassenMMult(0, m2, 0, n2, 0, p2, tAM2, B11, M2);

//M3 = A11*(B12 - B22)
SubMatBlocks(tBM3, p2, n2, B12, B22);
strassenMMult(0, m2, 0, n2, 0, p2, A11, tBM3, M3);

//M4 = A22*(B21 - B11)
SubMatBlocks(tBM4, p2, n2, B21, B11);
strassenMMult(0, m2, 0, n2, 0, p2, A22, tBM4, M4);

//M5 = (A11 + A12)*B22
AddMatBlocks(tAM5, m2, p2, A11, A12);
strassenMMult(0, m2, 0, n2, 0, p2, tAM5, B22, M5);

//M6 = (A21 - A11)*(B11 + B12)
SubMatBlocks(tAM6, m2, p2, A21, A11);
AddMatBlocks(tBM6, p2, n2, B11, B12);
strassenMMult(0, m2, 0, n2, 0, p2, tAM6, tBM6, M6);

//M7 = (A12 - A22)*(B21 + B22)
SubMatBlocks(tAM7, m2, p2, A12, A22);
AddMatBlocks(tBM7, p2, n2, B21, B22);
strassenMMult(0, m2, 0, n2, 0, p2, tAM7, tBM7, M7);

for (int i = 0; i < m2; i++)
for (int j = 0; j < n2; j++) {
C11[i][j] = M1[i][j] + M4[i][j] - M5[i][j] + M7[i][j];
C12[i][j] = M3[i][j] + M5[i][j];
C21[i][j] = M2[i][j] + M4[i][j];
C22[i][j] = M1[i][j] - M2[i][j] + M3[i][j] + M6[i][j];
}

Free2DArray< double >(M1);
Free2DArray< double >(M2);
Free2DArray< double >(M3);
Free2DArray< double >(M4);
Free2DArray< double >(M5);
Free2DArray< double >(M6);
Free2DArray< double >(M7);

delete[] A11; delete[] A12; delete[] A21; delete[] A22;
delete[] B11; delete[] B12; delete[] B21; delete[] B22;
delete[] C11; delete[] C12; delete[] C21; delete[] C22;

Free2DArray< double >(tAM1);
Free2DArray< double >(tBM1);
Free2DArray< double >(tAM2);
Free2DArray< double >(tBM3);
Free2DArray< double >(tBM4);
Free2DArray< double >(tAM5);
Free2DArray< double >(tAM6);
Free2DArray< double >(tBM6);
Free2DArray< double >(tAM7);
Free2DArray< double >(tBM7);
}
}

void matmultS(int m, int n, int p, double **A, double **B, double **C)
{
int i,j;
for (i=0; i < m; i++)
for (j=0; j < n; j++)
C[i][j] = 0;
strassenMMult(0, m, 0, n, 0, p, A, B, C);
}


int CheckResults(int m, int n, double **C, double **C1)
{
#define THRESHOLD 0.001
//
// May need to take into consideration the floating point roundoff error
// due to parallel execution
//
for (int i = 0; i < m; i++) {
for (int j = 0; j < n; j++) {
if (abs(C[i][j] - C1[i][j]) > THRESHOLD ) {
printf("%f %f\n", C[i][j], C1[i][j]);
return 1;
}
}
}
return 0;
}



int main(int argc, char* argv[])
{
clock_t before, after;

int M = atoi(argv[1]);
int N = atoi(argv[2]);
int P = atoi(argv[3]);

double **A = Allocate2DArray< double >(M, P);
double **B = Allocate2DArray< double >(P, N);
double **C = Allocate2DArray< double >(M, N);
double **C4 = Allocate2DArray< double >(M, N);

int i, j;

for (i = 0; i < M; i++) {
for (j = 0; j < P; j++) {
A[i][j] = 5.0 - ((double)(rand()%100) / 10.0);
}
}

for (i = 0; i < P; i++) {
for (j = 0; j < N; j++) {
B[i][j] = 5.0 - ((double)(rand()%100) / 10.0);
}
}

printf("Execute Standard matmult\n\n");
before = clock();
seqMatMult(M, N, P, A, B, C);
after = clock();
printf("Standard matrix function done in %7.2f secs\n\n\n",(float)(after - before)/ CLOCKS_PER_SEC);

before = clock();
matmultS(M, N, P, A, B, C4);
after = clock();
printf("Strassen matrix function done in %7.2f secs\n\n\n",(float)(after - before)/ CLOCKS_PER_SEC);

if (CheckResults(M, N, C, C4))
printf("Error in matmultS\n\n");
else
printf("OKAY\n\n");

Free2DArray(A);
Free2DArray(B);
Free2DArray(C);
Free2DArray(C4);

return 0;
}

Answer 2

C code of two 2 by 2 matrix multiplication using Strassen algorithm:

#include

int main(){

int a[2][2],b[2][2],c[2][2],i,j;

int m1,m2,m3,m4,m5,m6,m7;

printf("Enter the 4 elements of first matrix: ");

for(i=0;i<2;i++)

for(j=0;j<2;j++)

scanf("%d",&a[i][j]);

printf("Enter the 4 elements of second matrix: ");

for(i=0;i<2;i++)

for(j=0;j<2;j++)

scanf("%d",&b[i][j]);

printf("\nThe first matrix is\n");

for(i=0;i<2;i++){

printf("\n");

for(j=0;j<2;j++)

printf("%d\t",a[i][j]);

}

printf("\nThe second matrix is\n");

for(i=0;i<2;i++){

printf("\n");

for(j=0;j<2;j++)

printf("%d\t",b[i][j]);

}

m1= (a[0][0] + a[1][1])*(b[0][0]+b[1][1]);

m2= (a[1][0]+a[1][1])*b[0][0];

m3= a[0][0]*(b[0][1]-b[1][1]);

m4= a[1][1]*(b[1][0]-b[0][0]);

m5= (a[0][0]+a[0][1])*b[1][1];

m6= (a[1][0]-a[0][0])*(b[0][0]+b[0][1]);

m7= (a[0][1]-a[1][1])*(b[1][0]+b[1][1]);

c[0][0]=m1+m4-m5+m7;

c[0][1]=m3+m5;

c[1][0]=m2+m4;

c[1][1]=m1-m2+m3+m6;

printf("\nAfter multiplication using \n");

for(i=0;i<2;i++){

printf("\n");

for(j=0;j<2;j++)

printf("%d\t",c[i][j]);

}

return 0;

}

Answer 3

# include <iostream.h> # include <conio.h> int matrix_A[4][4]={0}; int matrix_B[4][4]={0}; int matrix_C[4][4]={0}; void get_matrix_A( ); void get_matrix_B( ); void multiply_matrices( ); void show_matrix_C( ); void add_2x2_matrices(constint [2][2],constint [2][2],int [2][2]); void subtract_2x2_matrices(constint [2][2],constint [2][2],int [2][2]); void multiply_2x2_matrices(constint [2][2],constint [2][2],int [2][2]); int main( ) { clrscr( ); textmode(C4350); get_matrix_A( ); get_matrix_B( ); multiply_matrices( ); show_matrix_C( ); cout<<"\n\n\n\n\n Press any key to exit..."; getch( ); return 0; } /***********************************************************************///-------------------------- get_matrix_A( ) ------------------------///***********************************************************************/void get_matrix_A( ) { gotoxy(1,2); cout<<" Enter the values of Matrix-A row by row :\n "<<endl; cout<<"\t\t\t Ú ¿"<<endl; cout<<"\t\t\t ³ ³"<<endl; cout<<"\t\t\t ³ ³"<<endl; cout<<"\t\t\t ³ ³"<<endl; cout<<"\t\t\t ³ ³"<<endl; cout<<"\t\t\t À Ù"<<endl; gotoxy(18,6); cout<<" A = "<<endl; int x=28; int y=5; for(int i=0;i<4;i++) { for(int j=0;j<4;j++) { gotoxy(x,y); cin>>matrix_A[i][j]; x+=5; } x=28; y++; } } /***********************************************************************///------------------------- get_matrix_B( ) -------------------------///***********************************************************************/void get_matrix_B( ) { gotoxy(1,15); cout<<" Enter the values of Matrix-B row by row :\n "<<endl; cout<<"\t\t\t Ú ¿"<<endl; cout<<"\t\t\t ³ ³"<<endl; cout<<"\t\t\t ³ ³"<<endl; cout<<"\t\t\t ³ ³"<<endl; cout<<"\t\t\t ³ ³"<<endl; cout<<"\t\t\t À Ù"<<endl; gotoxy(18,19); cout<<" B = "<<endl; int x=28; int y=18; for(int i=0;i<4;i++) { for(int j=0;j<4;j++) { gotoxy(x,y); cin>>matrix_B[i][j]; x+=5; } x=28; y++; } } /*************************************************************************///------------------------ add_2x2_matrices( ) ------------------------///*************************************************************************/void add_2x2_matrices(constint a[2][2],constint b[2][2],int c[2][2]) { for(int i=0;i<2;i++) { for(int j=0;j<2;j++) c[i][j]=(a[i][j]+b[i][j]); } } /*************************************************************************///----------------------- subtract_2x2_matrices( ) --------------------///*************************************************************************/void subtract_2x2_matrices(constint a[2][2],constint b[2][2],int c[2][2]) { for(int i=0;i<2;i++) { for(int j=0;j<2;j++) c[i][j]=(a[i][j]-b[i][j]); } } /*************************************************************************///---------------------- multiply_2x2_matrices( ) ---------------------///*************************************************************************/void multiply_2x2_matrices(constint a[2][2],constint b[2][2],int c[2][2]) { for(int i=0;i<2;i++) { for(int j=0;j<2;j++) { c[i][j]=0; for(int k=0;k<2;k++) c[i][j]+=(a[j][k]*b[k][j]); } } } /*************************************************************************///----------------------- multiply_matrices( ) ------------------------///*************************************************************************/void multiply_matrices( ) { int A11[2][2]={0}; int A12[2][2]={0}; int A21[2][2]={0}; int A22[2][2]={0}; int B11[2][2]={0}; int B12[2][2]={0}; int B21[2][2]={0}; int B22[2][2]={0}; int C11[2][2]={0}; int C12[2][2]={0}; int C21[2][2]={0}; int C22[2][2]={0}; int i; int j; for(i=0;i<2;i++) { for(j=0;j<2;j++) { A11[i][j]=matrix_A[i][j]; B11[i][j]=matrix_B[i][j]; } } for(i=0;i<2;i++) { for(j=2;j<4;j++) { A12[i][(j-2)]=matrix_A[i][j]; B12[i][(j-2)]=matrix_B[i][j]; } } for(i=2;i<4;i++) { for(j=0;j<2;j++) { A21[(i-2)][j]=matrix_A[i][j]; B21[(i-2)][j]=matrix_B[i][j]; } } for(i=2;i<4;i++) { for(j=2;j<4;j++) { A22[(i-2)][(j-2)]=matrix_A[i][j]; B22[(i-2)][(j-2)]=matrix_B[i][j]; } } int M1[2][2]={0}; int M2[2][2]={0}; int M3[2][2]={0}; int M4[2][2]={0}; int M5[2][2]={0}; int M6[2][2]={0}; int M7[2][2]={0}; int Temp1[2][2]={0}; int Temp2[2][2]={0}; subtract_2x2_matrices(A12,A22,Temp1); add_2x2_matrices(B21,B22,Temp2); multiply_2x2_matrices(Temp1,Temp2,M1); add_2x2_matrices(A11,A22,Temp1); add_2x2_matrices(B11,B22,Temp2); multiply_2x2_matrices(Temp1,Temp2,M2); subtract_2x2_matrices(A11,A21,Temp1); add_2x2_matrices(B11,B12,Temp2); multiply_2x2_matrices(Temp1,Temp2,M3); add_2x2_matrices(A11,A12,Temp1); multiply_2x2_matrices(Temp1,B22,M4); subtract_2x2_matrices(B12,B22,Temp1); multiply_2x2_matrices(Temp1,A11,M5); subtract_2x2_matrices(B21,B11,Temp1); multiply_2x2_matrices(Temp1,A22,M6); add_2x2_matrices(A21,A22,Temp1); multiply_2x2_matrices(Temp1,B11,M7); add_2x2_matrices(M1,M6,Temp1); subtract_2x2_matrices(M2,M4,Temp2); add_2x2_matrices(Temp1,Temp2,C11); add_2x2_matrices(M4,M5,C12); add_2x2_matrices(M6,M7,C21); subtract_2x2_matrices(M2,M3,Temp1); subtract_2x2_matrices(M5,M7,Temp2); add_2x2_matrices(Temp1,Temp2,C22); for(i=0;i<2;i++) { for(j=0;j<2;j++) matrix_C[i][j]=C11[i][j]; } for(i=0;i<2;i++) { for(j=2;j<4;j++) matrix_C[i][j]=C12[i][(j-2)]; } for(i=2;i<4;i++) { for(j=0;j<2;j++) matrix_C[i][j]=C21[(i-2)][j]; } for(i=2;i<4;i++) { for(j=2;j<4;j++) matrix_C[i][j]=C22[(i-2)][(j-2)]; } } /***********************************************************************///------------------------ show_matrix_C( ) -------------------------///***********************************************************************/void show_matrix_C( ) { gotoxy(1,28); cout<<" Values of Matrix-C row by row :\n "<<endl; cout<<"\t\t\t Ú ¿"<<endl; cout<<"\t\t\t ³ ³"<<endl; cout<<"\t\t\t ³ ³"<<endl; cout<<"\t\t\t ³ ³"<<endl; cout<<"\t\t\t ³ ³"<<endl; cout<<"\t\t\t À Ù"<<endl; gotoxy(18,32); cout<<" C = "<<endl; int x=28; int y=31; for(int i=0;i<4;i++) { for(int j=0;j<4;j++) { gotoxy(x,y); cout<<matrix_C[i][j]; x+=5; } x=28; y++; } }