Source Code of The Strassen Algorithm in C Language
/*------------------------------------------------------------------------------*/
/* Compile this without linking since there is no main method. */
/* Assuming that the file name is Strassen.c this can be done using gcc: */
/* gcc -c Strassen.c */
#include
#include
void strassen(double **a, double **b, double **c, int tam);
void sum(double **a, double **b, double **result, int tam);
void subtract(double **a, double **b, double **result, int tam);
double **allocate_real_matrix(int tam, int random);
double **free_real_matrix(double **v, int tam);
void strassen(double **a, double **b, double **c, int tam) { // trivial case: when the matrix is 1 X 1: if (tam == 1) { c = a * b; return; } // other cases are treated here: int newTam = tam/2; double **a11, **a12, **a21, **a22; double **b11, **b12, **b21, **b22; double **c11, **c12, **c21, **c22; double **p1, **p2, **p3, **p4, **p5, **p6, **p7; // memory allocation: a11 = allocate_real_matrix(newTam, -1); a12 = allocate_real_matrix(newTam, -1); a21 = allocate_real_matrix(newTam, -1); a22 = allocate_real_matrix(newTam, -1); b11 = allocate_real_matrix(newTam, -1); b12 = allocate_real_matrix(newTam, -1); b21 = allocate_real_matrix(newTam, -1); b22 = allocate_real_matrix(newTam, -1); c11 = allocate_real_matrix(newTam, -1); c12 = allocate_real_matrix(newTam, -1); c21 = allocate_real_matrix(newTam, -1); c22 = allocate_real_matrix(newTam, -1); p1 = allocate_real_matrix(newTam, -1); p2 = allocate_real_matrix(newTam, -1); p3 = allocate_real_matrix(newTam, -1); p4 = allocate_real_matrix(newTam, -1); p5 = allocate_real_matrix(newTam, -1); p6 = allocate_real_matrix(newTam, -1); p7 = allocate_real_matrix(newTam, -1); double **aResult = allocate_real_matrix(newTam, -1); double **bResult = allocate_real_matrix(newTam, -1); int i, j; //dividing the matrices in 4 sub-matrices: for (i = 0; i < newTam; i++) { for (j = 0; j < newTam; j++) { a11 = a; a12 = a; a21 = a; a22 = a; b11 = b; b12 = b; b21 = b; b22 = b; } } // Calculating p1 to p7: sum(a11, a22, aResult, newTam); // a11 + a22 sum(b11, b22, bResult, newTam); // b11 + b22 strassen(aResult, bResult, p1, newTam); // p1 = (a11+a22) * (b11+b22) sum(a21, a22, aResult, newTam); // a21 + a22 strassen(aResult, b11, p2, newTam); // p2 = (a21+a22) * (b11) subtract(b12, b22, bResult, newTam); // b12 - b22 strassen(a11, bResult, p3, newTam); // p3 = (a11) * (b12 - b22) subtract(b21, b11, bResult, newTam); // b21 - b11 strassen(a22, bResult, p4, newTam); // p4 = (a22) * (b21 - b11) sum(a11, a12, aResult, newTam); // a11 + a12 strassen(aResult, b22, p5, newTam); // p5 = (a11+a12) * (b22) subtract(a21, a11, aResult, newTam); // a21 - a11 sum(b11, b12, bResult, newTam); // b11 + b12 strassen(aResult, bResult, p6, newTam); // p6 = (a21-a11) * (b11+b12) subtract(a12, a22, aResult, newTam); // a12 - a22 sum(b21, b22, bResult, newTam); // b21 + b22 strassen(aResult, bResult, p7, newTam); // p7 = (a12-a22) * (b21+b22) // calculating c21, c21, c11 e c22: sum(p3, p5, c12, newTam); // c12 = p3 + p5 sum(p2, p4, c21, newTam); // c21 = p2 + p4 sum(p1, p4, aResult, newTam); // p1 + p4 sum(aResult, p7, bResult, newTam); // p1 + p4 + p7 subtract(bResult, p5, c11, newTam); // c11 = p1 + p4 - p5 + p7 sum(p1, p3, aResult, newTam); // p1 + p3 sum(aResult, p6, bResult, newTam); // p1 + p3 + p6 subtract(bResult, p2, c22, newTam); // c22 = p1 + p3 - p2 + p6 // Grouping the results obtained in a single matrix: for (i = 0; i < newTam ; i++) { for (j = 0 ; j < newTam ; j++) { c = c11; c = c12; c = c21; c = c22; } } // deallocating memory (free): a11 = free_real_matrix(a11, newTam); a12 = free_real_matrix(a12, newTam); a21 = free_real_matrix(a21, newTam); a22 = free_real_matrix(a22, newTam); b11 = free_real_matrix(b11, newTam); b12 = free_real_matrix(b12, newTam); b21 = free_real_matrix(b21, newTam); b22 = free_real_matrix(b22, newTam); c11 = free_real_matrix(c11, newTam); c12 = free_real_matrix(c12, newTam); c21 = free_real_matrix(c21, newTam); c22 = free_real_matrix(c22, newTam); p1 = free_real_matrix(p1, newTam); p2 = free_real_matrix(p2, newTam); p3 = free_real_matrix(p3, newTam); p4 = free_real_matrix(p4, newTam); p5 = free_real_matrix(p5, newTam); p6 = free_real_matrix(p6, newTam); p7 = free_real_matrix(p7, newTam); aResult = free_real_matrix(aResult, newTam); bResult = free_real_matrix(bResult, newTam);
} // end of Strassen function
/*------------------------------------------------------------------------------*/
// function to sum two matrices
void sum(double **a, double **b, double **result, int tam) { int i, j; for (i = 0; i < tam; i++) { for (j = 0; j < tam; j++) { result = a + b; } }
}
/*------------------------------------------------------------------------------*/
// function to subtract two matrices
void subtract(double **a, double **b, double **result, int tam) { int i, j; for (i = 0; i < tam; i++) { for (j = 0; j < tam; j++) { result = a - b; } }
}
/*------------------------------------------------------------------------------*/
// This function allocates the matrix using malloc, and initializes it. If the variable random is passed
// as zero, it initializes the matrix with zero, if it's passed as 1, it initializes the matrix with random
// values. If it is passed with any other int value (like -1 for example) the matrix is initialized with no
// values in it. The variable tam defines the length of the matrix.
double **allocate_real_matrix(int tam, int random) { int i, j, n = tam, m = tam; double **v, a; // pointer to the vector // allocates one vector of vectors (matrix) v = (double**) malloc(n * sizeof(double*)); if (v == NULL) { printf ("** Error in matrix allocation: insufficient memory **"); return (NULL); } // allocates each row of the matrix for (i = 0; i < n; i++) { v = (double*) malloc(m * sizeof(double)); if (v == NULL) { printf ("** Error: Insufficient memory **"); free_real_matrix(v, n); return (NULL); } // initializes the matrix with zeros if (random == 0) { for (j = 0; j < m; j++) v = 0.0; } // initializes the matrix with random values between 0 and 10 else { if (random == 1) { for (j = 0; j < m; j++) { a = rand; v = (a - (int)a) * 10; } } } } return (v); // returns the pointer to the vector.
}
/*------------------------------------------------------------------------------*/
// This function unallocated the matrix (frees memory)
double **free_real_matrix(double **v, int tam) { int i; if (v == NULL) { return (NULL); } for (i = 0; i < tam; i++) { if (v) { free(v); // frees a row of the matrix v = NULL; } } free(v); // frees the pointer / v = NULL; return (NULL); //returns a null pointer /
}
/*------------------------------------------------------------------------------*/