From 4be0bc64807215636602f9f742c9255169b253d9 Mon Sep 17 00:00:00 2001 From: dbjdbj Date: Thu, 22 Jul 2021 00:04:06 +0200 Subject: [PATCH 01/16] all changed nothing works --- .gitignore | 6 ++ Matrix.c | 73 +++++++------- Matrix.h | 25 ----- Strassen_Algorithm.c | 220 +++++++++++++++++-------------------------- 4 files changed, 130 insertions(+), 194 deletions(-) diff --git a/.gitignore b/.gitignore index bbf313b..5b10f1f 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,9 @@ +# DBJ +*.ilk +*.pdb +.vs +.vscode + # Object files *.o *.ko diff --git a/Matrix.c b/Matrix.c index 5cb41a5..3c798ed 100644 --- a/Matrix.c +++ b/Matrix.c @@ -1,7 +1,8 @@ +#include #include #include "Matrix.h" -static int **Matrix_Allocate(uint32_t row, uint32_t column) +static inline int **matrix_allocate_value(uint32_t row, uint32_t column) { int **matrix = (int **) calloc(row + row * column, sizeof(**matrix)); int *matrixTemp = (int *) (matrix + row); @@ -14,62 +15,60 @@ static int **Matrix_Allocate(uint32_t row, uint32_t column) return matrix; } -__attribute__((always_inline)) -inline void Matrix_Free(void *ptr) +static inline void Matrix_Free(Matrix **ptr) { - Matrix *matrix = (Matrix *) ptr; - return free(matrix->values); + assert(*ptr && (*ptr)->values); + free((*ptr)->values); + free(*ptr); + *ptr = 0 ; } -static Matrix Matrix_Addition(Matrix res, const Matrix a, const Matrix b) +#define matrix_autofree __attribute__((cleanup(Matrix_Free))) + +static inline Matrix * Matrix_Addition(Matrix * res, const Matrix * a, const Matrix * b) { - for (uint32_t i = 0; i < res.row; i++) - for (uint32_t j = 0; j < res.column; j++) - res.values[i][j] = a.values[i][j] + b.values[i][j]; + assert( res && a && b); + for (uint32_t i = 0; i < res->row; i++) + for (uint32_t j = 0; j < res->column; j++) + res->values[i][j] = a->values[i][j] + b->values[i][j]; return res; } -static Matrix Matrix_Subtract(Matrix res, const Matrix a, const Matrix b) +static inline Matrix * Matrix_Subtract(Matrix * res, const Matrix * a, const Matrix * b) { - for (uint32_t i = 0; i < res.row; i++) - for (uint32_t j = 0; j < res.column; j++) - res.values[i][j] = a.values[i][j] - b.values[i][j]; + assert( res && a && b); + for (uint32_t i = 0; i < res->row; i++) + for (uint32_t j = 0; j < res->column; j++) + res->values[i][j] = a->values[i][j] - b->values[i][j]; return res; } -static Matrix Matrix_Multiply(Matrix res, const Matrix a, const Matrix b) +static inline Matrix * Matrix_Multiply(Matrix * res, const Matrix * a, const Matrix * b) { - int m1 = (a.values[0][0] + a.values[1][1]) * (b.values[0][0] + b.values[1][1]); - int m2 = (a.values[1][0] + a.values[1][1]) * b.values[0][0]; - int m3 = (b.values[0][1] - b.values[1][1]) * a.values[0][0]; - int m4 = (b.values[1][0] - b.values[0][0]) * a.values[1][1]; - int m5 = (a.values[0][0] + a.values[0][1]) * b.values[1][1]; - int m6 = (a.values[1][0] - a.values[0][0]) * (b.values[0][0] + b.values[0][1]); - int m7 = (a.values[0][1] - a.values[1][1]) * (b.values[1][0] + b.values[1][1]); + assert( res && a && b); + int m1 = (a->values[0][0] + a->values[1][1]) * (b->values[0][0] + b->values[1][1]); + int m2 = (a->values[1][0] + a->values[1][1]) * b->values[0][0]; + int m3 = (b->values[0][1] - b->values[1][1]) * a->values[0][0]; + int m4 = (b->values[1][0] - b->values[0][0]) * a->values[1][1]; + int m5 = (a->values[0][0] + a->values[0][1]) * b->values[1][1]; + int m6 = (a->values[1][0] - a->values[0][0]) * (b->values[0][0] + b->values[0][1]); + int m7 = (a->values[0][1] - a->values[1][1]) * (b->values[1][0] + b->values[1][1]); - res.values[0][0] = m1 + m4 - m5 + m7; - res.values[0][1] = m3 + m5; - res.values[1][0] = m2 + m4; - res.values[1][1] = m1 - m2 + m3 + m6; + res->values[0][0] = m1 + m4 - m5 + m7; + res->values[0][1] = m3 + m5; + res->values[1][0] = m2 + m4; + res->values[1][1] = m1 - m2 + m3 + m6; return res; } -void Matrix_Initializer(Matrix *matrix, uint32_t row, uint32_t column) +static inline Matrix * Matrix_Initializer( uint32_t row, uint32_t column) { + Matrix *matrix = calloc(1, sizeof(Matrix)); matrix->row = row; matrix->column = column; - matrix->values = NULL; - matrix->New = Matrix_Allocate; - matrix->Delete = free; - - matrix->values = matrix->New(row, column); + matrix->values = matrix_allocate_value(row, column); + return matrix ; } - -Matrix_Arith Naive_Matrix_Arith __attribute__((section("Matrix_Arith"))) = { - .Addition = Matrix_Addition, - .Subtract = Matrix_Subtract, - .Multiply = Matrix_Multiply, -}; diff --git a/Matrix.h b/Matrix.h index 6521458..1932de2 100644 --- a/Matrix.h +++ b/Matrix.h @@ -3,38 +3,13 @@ #include -#define MATRIX_ARITH_BEGIN __start_Matrix_Arith -#define MATRIX_ARITH_END __stop_Matrix_Arith - -#define MATRIX_INITIALIZER(X, ROW, COLUMN) \ - Matrix_Initializer(&(X) ,(ROW), (COLUMN)) - -#define autofree \ - __attribute__((cleanup(Matrix_Free))) - enum { NAIVE_ARITHMETIC, }; typedef struct _Matrix { uint32_t row; uint32_t column; int **values; - - int **(*New)(uint32_t row, uint32_t column); - void (*Delete)(void *); } Matrix; -typedef struct _Matrix_Arith { - Matrix (*Addition)(Matrix, const Matrix, const Matrix); - Matrix (*Subtract)(Matrix, const Matrix, const Matrix); - Matrix (*Multiply)(Matrix, const Matrix, const Matrix); -} Matrix_Arith; - -void Matrix_Initializer(Matrix *, uint32_t, uint32_t); -void Matrix_Free(void *); - -extern Matrix_Arith Naive_Matrix_Arith; - -extern Matrix_Arith __start_Matrix_Arith[]; -extern Matrix_Arith __stop_Matrix_Arith[]; #endif /* MATRIX_H_ */ diff --git a/Strassen_Algorithm.c b/Strassen_Algorithm.c index f9b7b7b..3227c60 100644 --- a/Strassen_Algorithm.c +++ b/Strassen_Algorithm.c @@ -1,164 +1,120 @@ -#ifdef _cplusplus -extern "C" { - #include - #include - #include - #include "Matrix.h" -} -#else + #include #include - #include - #include "Matrix.h" -#endif + #include + #include "Matrix.c" -#define LOG2(LENGTH) \ - (8 * sizeof(unsigned int) - Count_Leading_Zero((LENGTH)) - 1) +#define N 64 -int Count_Leading_Zero(unsigned int number); - -Matrix Strassen(Matrix, const Matrix, const Matrix, int); - -int main(int argc, char *argv[]) +static inline Matrix * strassen(Matrix * dest, const Matrix * srcA, const Matrix * srcB, int length) { - int dimensions, matrix_length; + + if (length == 2) return Matrix_Multiply(dest, srcA, srcB); - printf("Please enter the dimensions of the matrix: "); - scanf("%d", &dimensions); + int len = length / 2; - if (dimensions <= 0) { printf("The number you entered is invalid."); exit(1); } + matrix_autofree Matrix * a11= Matrix_Initializer( len, len), * a12= Matrix_Initializer( len, len), * a21= Matrix_Initializer( len, len), * a22= Matrix_Initializer( len, len), + * b11= Matrix_Initializer( len, len), * b12= Matrix_Initializer( len, len), * b21= Matrix_Initializer( len, len), * b22= Matrix_Initializer( len, len), * c11= Matrix_Initializer( len, len), * c12= Matrix_Initializer( len, len), * c21= Matrix_Initializer( len, len), * c22= Matrix_Initializer( len, len), * + m1= Matrix_Initializer( len, len), * m2= Matrix_Initializer( len, len), * m3= Matrix_Initializer( len, len), * m4= Matrix_Initializer( len, len), * m5= Matrix_Initializer( len, len), * m6= Matrix_Initializer( len, len), * m7= Matrix_Initializer( len, len), * temp1= Matrix_Initializer( len, len), * temp2 = Matrix_Initializer( len, len) ; - /* Check if dimensions of matrix is the power of two or not */ - if ((dimensions & (dimensions - 1)) || (dimensions == 1)) - matrix_length = 2 << LOG2((unsigned int) dimensions); - else - matrix_length = dimensions; + /* Divide matrix into four parts */ + for (int i = 0; i < len; i++) { + for (int j = 0; j < len; j++) { + a11->values[i][j] = srcA->values[i][j]; + a12->values[i][j] = srcA->values[i][j + len]; + a21->values[i][j] = srcA->values[i + len][j]; + a22->values[i][j] = srcA->values[i + len][j + len]; + + b11->values[i][j] = srcB->values[i][j]; + b12->values[i][j] = srcB->values[i][j + len]; + b21->values[i][j] = srcB->values[i + len][j]; + b22->values[i][j] = srcB->values[i + len][j + len]; + } + } - Matrix matrixA, matrixB, matrixC; + /* Calculate seven formulas of strassen Algorithm */ + strassen(m1, Matrix_Addition(temp1, a11, a22), Matrix_Addition(temp2, b11, b22), len); + strassen(m2, Matrix_Addition(temp1, a21, a22), b11, len); + strassen(m3, a11, Matrix_Subtract(temp1, b12, b22), len); + strassen(m4, a22, Matrix_Subtract(temp1, b21, b11), len); + strassen(m5, Matrix_Addition(temp1, a11, a12), b22, len); + strassen(m6, Matrix_Subtract(temp1, a21, a11), Matrix_Addition(temp2, b11, b12), len); + strassen(m7, Matrix_Subtract(temp1, a12, a22), Matrix_Addition(temp2, b21, b22), len); - MATRIX_INITIALIZER(matrixA, matrix_length, matrix_length); - MATRIX_INITIALIZER(matrixB, matrix_length, matrix_length); - MATRIX_INITIALIZER(matrixC, matrix_length, matrix_length); + /* Merge the answer of matrix dest */ + /* c11 = m1 + m4 - m5 + m7 = m1 + m4 - (m5 - m7) */ + Matrix_Subtract(c11, Matrix_Addition(temp1, m1, m4), Matrix_Subtract(temp2, m5, m7)); + Matrix_Addition(c12, m3, m5); + Matrix_Addition(c21, m2, m4); + Matrix_Addition(c22, Matrix_Subtract(temp1, m1, m2), Matrix_Addition(temp2, m3, m6)); - srand((unsigned int) time(NULL)); - for (int i = 0; i < dimensions; i++) { - for (int j = 0; j < dimensions; j++) { - matrixA.values[i][j] = rand() % 1000; - matrixB.values[i][j] = rand() % 1000; + /* Store the answer of matrix multiplication */ + for (int i = 0; i < len; i++) { + for (int j = 0; j < len; j++) { + dest->values[i][j] = c11->values[i][j]; + dest->values[i][j + len] = c12->values[i][j]; + dest->values[i + len][j] = c21->values[i][j]; + dest->values[i + len][j + len] = c22->values[i][j]; } } - /* Matrix multiplication */ - Strassen(matrixC, matrixA, matrixB, matrixC.row); + return dest; +} - /* Print the answer of matrix multiplication */ - printf("Matrix A:\n"); - for(int i = 0; i < dimensions; i++) { - for (int j = 0; j < dimensions; j++) { - printf("%5d ", matrixA.values[i][j]); +static void matrix_print ( const char * prompt, Matrix * mx ) +{ + assert( prompt && mx ); + printf("%s\n", prompt ); + for(int i = 0; i < mx->row ; i++) { + for (int j = 0; j < mx->column; j++) { + printf(" %5d ", mx->values[i][j]); } printf("\n"); } +} - printf("\nMatrix B:\n"); - for(int i = 0; i < dimensions; i++) { - for (int j = 0; j < dimensions; j++) { - printf("%5d ", matrixB.values[i][j]); +static void matrix_foreach ( Matrix * mx, Matrix * (*callback)(Matrix *, unsigned, unsigned) ) +{ + for(int i = 0; i < mx->row ; i++) { + for (int j = 0; j < mx->column; j++) { + callback(mx,i,j); } - printf("\n"); } - - printf("\nMatrix C:\n"); - for(int i = 0; i < dimensions; i++) { - for (int j = 0; j < dimensions; j++) { - printf("%10d ", matrixC.values[i][j]); - } - printf("\n"); - } - - matrixA.Delete(matrixA.values); - matrixB.Delete(matrixB.values); - matrixC.Delete(matrixC.values); - - return 0; } -int Count_Leading_Zero(unsigned int number) +Matrix * ordinal_as_val (Matrix * mx, unsigned i, unsigned j) { - if (number == 0) return 32; - - int count = 0; - if (number <= 0x0000FFFF) { count += 16; number <<= 16; } - if (number <= 0x00FFFFFF) { count += 8; number <<= 8; } - if (number <= 0x0FFFFFFF) { count += 4; number <<= 4; } - if (number <= 0x3FFFFFFF) { count += 2; number <<= 2; } - if (number <= 0x7FFFFFFF) { count += 1; number <<= 1; } - - return count; + mx->values[i][j] = (i * mx->column + j) ; + return mx ; } -Matrix Strassen(Matrix dest, const Matrix srcA, const Matrix srcB, int length) +/* +-------------------------------------------------------------------------------------------------------------------- +*/ +int main(int argc, char *argv[]) { - Matrix_Arith *arith = &__start_Matrix_Arith[NAIVE_ARITHMETIC]; - - if (length == 2) return arith->Multiply(dest, srcA, srcB); - - int len = length / 2; - - autofree Matrix a11, a12, a21, a22, /* Matrix srcA divide four part */ - b11, b12, b21, b22, /* Matrix srcB divide four part */ - c11, c12, c21, c22, /* Matrix dest divide four part */ - m1, m2, m3, m4, m5, m6, m7, - temp1, temp2; - - /* Initializer the matrix */ - MATRIX_INITIALIZER(a11, len, len); MATRIX_INITIALIZER(a12, len, len); MATRIX_INITIALIZER(a21, len, len); MATRIX_INITIALIZER(a22, len, len); - MATRIX_INITIALIZER(b11, len, len); MATRIX_INITIALIZER(b12, len, len); MATRIX_INITIALIZER(b21, len, len); MATRIX_INITIALIZER(b22, len, len); - MATRIX_INITIALIZER(c11, len, len); MATRIX_INITIALIZER(c12, len, len); MATRIX_INITIALIZER(c21, len, len); MATRIX_INITIALIZER(c22, len, len); - MATRIX_INITIALIZER(m1, len, len); MATRIX_INITIALIZER(m2, len, len); MATRIX_INITIALIZER(m3, len, len); MATRIX_INITIALIZER(m4, len, len); - MATRIX_INITIALIZER(m5, len, len); MATRIX_INITIALIZER(m6, len, len); MATRIX_INITIALIZER(m7, len, len); - MATRIX_INITIALIZER(temp1, len, len); MATRIX_INITIALIZER(temp2, len, len); - - /* Divide matrix to four part */ - for (int i = 0; i < len; i++) { - for (int j = 0; j < len; j++) { - a11.values[i][j] = srcA.values[i][j]; - a12.values[i][j] = srcA.values[i][j + len]; - a21.values[i][j] = srcA.values[i + len][j]; - a22.values[i][j] = srcA.values[i + len][j + len]; - - b11.values[i][j] = srcB.values[i][j]; - b12.values[i][j] = srcB.values[i][j + len]; - b21.values[i][j] = srcB.values[i + len][j]; - b22.values[i][j] = srcB.values[i + len][j + len]; - } + const unsigned matrix_length = N ; + + /* Check if dimensions of matrix is the power of two or not */ + if (ceil(log2(matrix_length)) != floor(log2(matrix_length))) + { + printf("\nERROR: square matrix side must be a power of 2. And %d is not.",matrix_length ); + return 0; } - /* Calculate seven formulas of Strassen Algorithm */ - Strassen(m1, arith->Addition(temp1, a11, a22), arith->Addition(temp2, b11, b22), len); - Strassen(m2, arith->Addition(temp1, a21, a22), b11, len); - Strassen(m3, a11, arith->Subtract(temp1, b12, b22), len); - Strassen(m4, a22, arith->Subtract(temp1, b21, b11), len); - Strassen(m5, arith->Addition(temp1, a11, a12), b22, len); - Strassen(m6, arith->Subtract(temp1, a21, a11), arith->Addition(temp2, b11, b12), len); - Strassen(m7, arith->Subtract(temp1, a12, a22), arith->Addition(temp2, b21, b22), len); - - /* Merge the answer of matrix dest */ - /* c11 = m1 + m4 - m5 + m7 = m1 + m4 - (m5 - m7) */ - arith->Subtract(c11, arith->Addition(temp1, m1, m4), arith->Subtract(temp2, m5, m7)); - arith->Addition(c12, m3, m5); - arith->Addition(c21, m2, m4); - arith->Addition(c22, arith->Subtract(temp1, m1, m2), arith->Addition(temp2, m3, m6)); + matrix_autofree Matrix * matrixA = Matrix_Initializer( matrix_length, matrix_length); + matrix_autofree Matrix * matrixB = Matrix_Initializer( matrix_length, matrix_length); + matrix_autofree Matrix * matrixC = Matrix_Initializer( matrix_length, matrix_length); - /* Store the answer of matrix multiplication */ - for (int i = 0; i < len; i++) { - for (int j = 0; j < len; j++) { - dest.values[i][j] = c11.values[i][j]; - dest.values[i][j + len] = c12.values[i][j]; - dest.values[i + len][j] = c21.values[i][j]; - dest.values[i + len][j + len] = c22.values[i][j]; - } - } + matrix_foreach(matrixA, ordinal_as_val ); + matrix_foreach(matrixB, ordinal_as_val ); - return dest; + /* Matrix multiplication */ + matrixC = strassen(matrixC, matrixA, matrixB, matrix_length); + +matrix_print("Matrix A:", matrixA ); +matrix_print("Matrix B:", matrixB ); +matrix_print("Matrix C:", matrixC ); + + return 0; } From f67ee73d7132bdf74e48a3b230efeb0e3a6e7364 Mon Sep 17 00:00:00 2001 From: dbjdbj Date: Thu, 22 Jul 2021 10:04:02 +0200 Subject: [PATCH 02/16] much better but strasses is NOT --- Matrix.c | 74 --------------------- Matrix.h | 15 ----- Strassen_Algorithm.c | 151 ++++++++++++++++++++++++------------------- dbj_matrix.h | 104 +++++++++++++++++++++++++++++ 4 files changed, 190 insertions(+), 154 deletions(-) delete mode 100644 Matrix.c delete mode 100644 Matrix.h create mode 100644 dbj_matrix.h diff --git a/Matrix.c b/Matrix.c deleted file mode 100644 index 3c798ed..0000000 --- a/Matrix.c +++ /dev/null @@ -1,74 +0,0 @@ -#include -#include -#include "Matrix.h" - -static inline int **matrix_allocate_value(uint32_t row, uint32_t column) -{ - int **matrix = (int **) calloc(row + row * column, sizeof(**matrix)); - int *matrixTemp = (int *) (matrix + row); - - for (uint32_t i = 0; i < row; i++) { - matrix[i] = matrixTemp; - matrixTemp += column; - } - - return matrix; -} - -static inline void Matrix_Free(Matrix **ptr) -{ - assert(*ptr && (*ptr)->values); - free((*ptr)->values); - free(*ptr); - *ptr = 0 ; -} - -#define matrix_autofree __attribute__((cleanup(Matrix_Free))) - -static inline Matrix * Matrix_Addition(Matrix * res, const Matrix * a, const Matrix * b) -{ - assert( res && a && b); - for (uint32_t i = 0; i < res->row; i++) - for (uint32_t j = 0; j < res->column; j++) - res->values[i][j] = a->values[i][j] + b->values[i][j]; - - return res; -} - -static inline Matrix * Matrix_Subtract(Matrix * res, const Matrix * a, const Matrix * b) -{ - assert( res && a && b); - for (uint32_t i = 0; i < res->row; i++) - for (uint32_t j = 0; j < res->column; j++) - res->values[i][j] = a->values[i][j] - b->values[i][j]; - - return res; -} - -static inline Matrix * Matrix_Multiply(Matrix * res, const Matrix * a, const Matrix * b) -{ - assert( res && a && b); - int m1 = (a->values[0][0] + a->values[1][1]) * (b->values[0][0] + b->values[1][1]); - int m2 = (a->values[1][0] + a->values[1][1]) * b->values[0][0]; - int m3 = (b->values[0][1] - b->values[1][1]) * a->values[0][0]; - int m4 = (b->values[1][0] - b->values[0][0]) * a->values[1][1]; - int m5 = (a->values[0][0] + a->values[0][1]) * b->values[1][1]; - int m6 = (a->values[1][0] - a->values[0][0]) * (b->values[0][0] + b->values[0][1]); - int m7 = (a->values[0][1] - a->values[1][1]) * (b->values[1][0] + b->values[1][1]); - - res->values[0][0] = m1 + m4 - m5 + m7; - res->values[0][1] = m3 + m5; - res->values[1][0] = m2 + m4; - res->values[1][1] = m1 - m2 + m3 + m6; - - return res; -} - -static inline Matrix * Matrix_Initializer( uint32_t row, uint32_t column) -{ - Matrix *matrix = calloc(1, sizeof(Matrix)); - matrix->row = row; - matrix->column = column; - matrix->values = matrix_allocate_value(row, column); - return matrix ; -} diff --git a/Matrix.h b/Matrix.h deleted file mode 100644 index 1932de2..0000000 --- a/Matrix.h +++ /dev/null @@ -1,15 +0,0 @@ -#ifndef MATRIX_H_ -#define MATRIX_H_ - -#include - -enum { NAIVE_ARITHMETIC, }; - -typedef struct _Matrix { - uint32_t row; - uint32_t column; - int **values; -} Matrix; - - -#endif /* MATRIX_H_ */ diff --git a/Strassen_Algorithm.c b/Strassen_Algorithm.c index 3227c60..085891f 100644 --- a/Strassen_Algorithm.c +++ b/Strassen_Algorithm.c @@ -1,91 +1,106 @@ - #include - #include - #include - #include "Matrix.c" +#include +#include +#include +#include "dbj_matrix.h" -#define N 64 +#define N 2 -static inline Matrix * strassen(Matrix * dest, const Matrix * srcA, const Matrix * srcB, int length) +static inline Matrix *strassen(Matrix *dest, const Matrix *srcA, const Matrix *srcB, int length) { - - if (length == 2) return Matrix_Multiply(dest, srcA, srcB); + + if (length == 2) + return matrix_ijk_matmul(dest, srcA, srcB); int len = length / 2; - matrix_autofree Matrix * a11= Matrix_Initializer( len, len), * a12= Matrix_Initializer( len, len), * a21= Matrix_Initializer( len, len), * a22= Matrix_Initializer( len, len), - * b11= Matrix_Initializer( len, len), * b12= Matrix_Initializer( len, len), * b21= Matrix_Initializer( len, len), * b22= Matrix_Initializer( len, len), * c11= Matrix_Initializer( len, len), * c12= Matrix_Initializer( len, len), * c21= Matrix_Initializer( len, len), * c22= Matrix_Initializer( len, len), * - m1= Matrix_Initializer( len, len), * m2= Matrix_Initializer( len, len), * m3= Matrix_Initializer( len, len), * m4= Matrix_Initializer( len, len), * m5= Matrix_Initializer( len, len), * m6= Matrix_Initializer( len, len), * m7= Matrix_Initializer( len, len), * temp1= Matrix_Initializer( len, len), * temp2 = Matrix_Initializer( len, len) ; + matrix_autofree Matrix *a11 = matrix_new(len, len), *a12 = matrix_new(len, len), *a21 = matrix_new(len, len), *a22 = matrix_new(len, len), + *b11 = matrix_new(len, len), *b12 = matrix_new(len, len), *b21 = matrix_new(len, len), *b22 = matrix_new(len, len), *c11 = matrix_new(len, len), *c12 = matrix_new(len, len), *c21 = matrix_new(len, len), *c22 = matrix_new(len, len), *m1 = matrix_new(len, len), *m2 = matrix_new(len, len), *m3 = matrix_new(len, len), *m4 = matrix_new(len, len), *m5 = matrix_new(len, len), *m6 = matrix_new(len, len), *m7 = matrix_new(len, len), *temp1 = matrix_new(len, len), *temp2 = matrix_new(len, len); /* Divide matrix into four parts */ - for (int i = 0; i < len; i++) { - for (int j = 0; j < len; j++) { - a11->values[i][j] = srcA->values[i][j]; - a12->values[i][j] = srcA->values[i][j + len]; - a21->values[i][j] = srcA->values[i + len][j]; - a22->values[i][j] = srcA->values[i + len][j + len]; - - b11->values[i][j] = srcB->values[i][j]; - b12->values[i][j] = srcB->values[i][j + len]; - b21->values[i][j] = srcB->values[i + len][j]; - b22->values[i][j] = srcB->values[i + len][j + len]; + for (int i = 0; i < len; i++) + { + for (int j = 0; j < len; j++) + { + MXY(a11, i, j) = MXY(srcA, i, j); + MXY(a12, i, j) = MXY(srcA, i, j + len); + MXY(a21, i, j) = MXY(srcA, i + len, j); + MXY(a22, i, j) = MXY(srcA, i + len, j + len); + + MXY(b11, i, j) = MXY(srcB, i, j); + MXY(b12, i, j) = MXY(srcB, i, j + len); + MXY(b21, i, j) = MXY(srcB, i + len, j); + MXY(b22, i, j) = MXY(srcB, i + len, j + len); } } /* Calculate seven formulas of strassen Algorithm */ - strassen(m1, Matrix_Addition(temp1, a11, a22), Matrix_Addition(temp2, b11, b22), len); - strassen(m2, Matrix_Addition(temp1, a21, a22), b11, len); - strassen(m3, a11, Matrix_Subtract(temp1, b12, b22), len); - strassen(m4, a22, Matrix_Subtract(temp1, b21, b11), len); - strassen(m5, Matrix_Addition(temp1, a11, a12), b22, len); - strassen(m6, Matrix_Subtract(temp1, a21, a11), Matrix_Addition(temp2, b11, b12), len); - strassen(m7, Matrix_Subtract(temp1, a12, a22), Matrix_Addition(temp2, b21, b22), len); + strassen(m1, matrix_addition(temp1, a11, a22), matrix_addition(temp2, b11, b22), len); + strassen(m2, matrix_addition(temp1, a21, a22), b11, len); + strassen(m3, a11, matrix_subtraction(temp1, b12, b22), len); + strassen(m4, a22, matrix_subtraction(temp1, b21, b11), len); + strassen(m5, matrix_addition(temp1, a11, a12), b22, len); + strassen(m6, matrix_subtraction(temp1, a21, a11), matrix_addition(temp2, b11, b12), len); + strassen(m7, matrix_subtraction(temp1, a12, a22), matrix_addition(temp2, b21, b22), len); /* Merge the answer of matrix dest */ /* c11 = m1 + m4 - m5 + m7 = m1 + m4 - (m5 - m7) */ - Matrix_Subtract(c11, Matrix_Addition(temp1, m1, m4), Matrix_Subtract(temp2, m5, m7)); - Matrix_Addition(c12, m3, m5); - Matrix_Addition(c21, m2, m4); - Matrix_Addition(c22, Matrix_Subtract(temp1, m1, m2), Matrix_Addition(temp2, m3, m6)); - + matrix_subtraction(c11, matrix_addition(temp1, m1, m4), matrix_subtraction(temp2, m5, m7)); + matrix_addition(c12, m3, m5); + matrix_addition(c21, m2, m4); + matrix_addition(c22, matrix_subtraction(temp1, m1, m2), matrix_addition(temp2, m3, m6)); + /* Store the answer of matrix multiplication */ - for (int i = 0; i < len; i++) { - for (int j = 0; j < len; j++) { - dest->values[i][j] = c11->values[i][j]; - dest->values[i][j + len] = c12->values[i][j]; - dest->values[i + len][j] = c21->values[i][j]; - dest->values[i + len][j + len] = c22->values[i][j]; + for (int i = 0; i < len; i++) + { + for (int j = 0; j < len; j++) + { + MXY(dest, i, j) = MXY(c11, i, j); + MXY(dest, i, j + len) = MXY(c12, i, j); + MXY(dest, i + len, j) = MXY(c21, i, j); + MXY(dest, i + len, j + len) = MXY(c22, i, j); } } return dest; } -static void matrix_print ( const char * prompt, Matrix * mx ) +static void matrix_print(const char *prompt, Matrix *mx) { - assert( prompt && mx ); - printf("%s\n", prompt ); - for(int i = 0; i < mx->row ; i++) { - for (int j = 0; j < mx->column; j++) { - printf(" %5d ", mx->values[i][j]); + assert(prompt && mx); + printf("%s\n", prompt); + for (int i = 0; i < mx->rows; i++) + { + for (int j = 0; j < mx->cols; j++) + { + printf(" %5d ", MXY(mx, i, j)); } printf("\n"); } } -static void matrix_foreach ( Matrix * mx, Matrix * (*callback)(Matrix *, unsigned, unsigned) ) +static void inline matrix_foreach(Matrix *mx, Matrix *(*callback)(Matrix *, unsigned, unsigned)) { - for(int i = 0; i < mx->row ; i++) { - for (int j = 0; j < mx->column; j++) { - callback(mx,i,j); + for (int i = 0; i < mx->rows; i++) + { + for (int j = 0; j < mx->cols; j++) + { + callback(mx, i, j); } } } -Matrix * ordinal_as_val (Matrix * mx, unsigned i, unsigned j) +static inline Matrix *ordinal_as_val(Matrix *mx, unsigned i, unsigned j) { - mx->values[i][j] = (i * mx->column + j) ; - return mx ; + MXY(mx, i, j) = (i * mx->cols + j); + return mx; +} + +static inline Matrix *zoro(Matrix *mx, unsigned i, unsigned j) +{ + assert(mx && mx->values); + MXY(mx, i, j) = 0; + return mx; } /* @@ -93,28 +108,34 @@ Matrix * ordinal_as_val (Matrix * mx, unsigned i, unsigned j) */ int main(int argc, char *argv[]) { - const unsigned matrix_length = N ; - + const unsigned matrix_length = N; + /* Check if dimensions of matrix is the power of two or not */ if (ceil(log2(matrix_length)) != floor(log2(matrix_length))) { - printf("\nERROR: square matrix side must be a power of 2. And %d is not.",matrix_length ); + printf("\n%s\tERROR: square matrix side must be a power of 2. And %d is not.", argv[0], matrix_length); return 0; } - matrix_autofree Matrix * matrixA = Matrix_Initializer( matrix_length, matrix_length); - matrix_autofree Matrix * matrixB = Matrix_Initializer( matrix_length, matrix_length); - matrix_autofree Matrix * matrixC = Matrix_Initializer( matrix_length, matrix_length); - - matrix_foreach(matrixA, ordinal_as_val ); - matrix_foreach(matrixB, ordinal_as_val ); + matrix_autofree Matrix *matrixA = matrix_new(matrix_length, matrix_length); + matrix_autofree Matrix *matrixB = matrix_new(matrix_length, matrix_length); + matrix_autofree Matrix *matrixC = matrix_new(matrix_length, matrix_length); + matrix_autofree Matrix *matrixR = matrix_new(matrix_length, matrix_length); + + matrix_foreach(matrixA, ordinal_as_val); + matrix_foreach(matrixB, ordinal_as_val); + matrix_foreach(matrixC, zoro); + matrix_foreach(matrixR, zoro); /* Matrix multiplication */ matrixC = strassen(matrixC, matrixA, matrixB, matrix_length); -matrix_print("Matrix A:", matrixA ); -matrix_print("Matrix B:", matrixB ); -matrix_print("Matrix C:", matrixC ); + matrixR = matrix_ijk_matmul(matrixR, matrixA, matrixB); + + matrix_print("Matrix A:", matrixA); + matrix_print("Matrix B:", matrixB); + matrix_print("Matrix C:", matrixC); + matrix_print("Matrix R:", matrixR); return 0; } diff --git a/dbj_matrix.h b/dbj_matrix.h new file mode 100644 index 0000000..74a34e3 --- /dev/null +++ b/dbj_matrix.h @@ -0,0 +1,104 @@ +#ifndef MATRIX_H_ +#define MATRIX_H_ + +#include +#include +#include + +typedef struct _Matrix +{ + unsigned rows; + unsigned cols; + int *values; +} Matrix; + +#define MXY_NOT_POINTER_(M_, R_, C_) ((int(*)[M_.rows])(M_.values))[R_][C_] + +// mx as a pointer +#define MXY(M_, R_, C_) MXY_NOT_POINTER_(((Matrix)*M_), R_, C_) + +#define FOR(C_, N_) for (unsigned C_ = 0; C_ < N_; ++C_) + +/* + Matrix mx = { 3,2, calloc(3 * 2, sizeof(int)) } ; + + MXY(mx,1,1) = 42; + + assert( MXY(mx,1,1) == 42) ; +*/ + +static inline void matrix_free(Matrix **ptr) +{ + assert(*ptr && (*ptr)->values); + free((*ptr)->values); + free(*ptr); + *ptr = 0; +} + +#define matrix_autofree __attribute__((cleanup(matrix_free))) + +static inline Matrix *matrix_addition(Matrix *res, const Matrix *a, const Matrix *b) +{ + assert(res && a && b); + for (uint32_t i = 0; i < res->rows; i++) + for (uint32_t j = 0; j < res->cols; j++) + // res->values[i][j] = a->values[i][j] + b->values[i][j]; + MXY(res, i, j) = MXY(a, i, j) + MXY(b, i, j); + + return res; +} + +static inline Matrix *matrix_subtraction(Matrix *res, const Matrix *a, const Matrix *b) +{ + assert(res && a && b); + for (uint32_t i = 0; i < res->rows; i++) + for (uint32_t j = 0; j < res->cols; j++) + MXY(res, i, j) = MXY(a, i, j) - MXY(b, i, j); + + return res; +} + +static inline int matrix_equal(const Matrix *a, const Matrix *b) +{ + assert(a && b); + + if (a->rows != b->rows) + return 0; + if (a->cols != b->cols) + return 0; + + FOR(i, a->rows) + FOR(j, a->cols) + if (MXY(a, i, j) != MXY(b, i, j)) + return 0; + + return 1; +} + +static Matrix *matrix_ijk_matmul(Matrix *C, const Matrix *A, const Matrix *B) +{ + assert(A && B && C); + for (size_t i = 0; i < A->rows; ++i) + { + for (size_t j = 0; j < B->cols; ++j) + { + MXY(C, i, j) = 0; + for (size_t k = 0; k < A->cols; ++k) + { + MXY(C, i, j) += MXY(A, i, k) * MXY(B, k, j); + } + } + } + return C; +} + +static inline Matrix *matrix_new(uint32_t rows, uint32_t cols) +{ + Matrix *matrix = calloc(1, sizeof(Matrix)); + matrix->rows = rows; + matrix->cols = cols; + matrix->values = calloc(rows * cols, sizeof(int)); + return matrix; +} + +#endif /* MATRIX_H_ */ From e1c19ce3355de7277387b15e69bda640b4f6a1fa Mon Sep 17 00:00:00 2001 From: dbjdbj Date: Thu, 22 Jul 2021 10:10:05 +0200 Subject: [PATCH 03/16] updated --- README.md | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/README.md b/README.md index 1ad5120..1a9597c 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,13 @@ -# Strassen Algorithm +# 100% REFACTORED + +Unfortunately Strassen as coded here is not working. + +Refactoring (c) 20201 by dbj@dbj.org -- https://dbj.org/license_dbj +--- + +Original README + +### Strassen Algorithm Implementation of the Strassen Algorithm in C From e9944bfd7dbe433f3ed6be32f13724588eb85031 Mon Sep 17 00:00:00 2001 From: dbjdbj Date: Thu, 22 Jul 2021 10:10:38 +0200 Subject: [PATCH 04/16] formated --- README.md | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/README.md b/README.md index 1a9597c..d30356d 100644 --- a/README.md +++ b/README.md @@ -2,7 +2,9 @@ Unfortunately Strassen as coded here is not working. -Refactoring (c) 20201 by dbj@dbj.org -- https://dbj.org/license_dbj +> Refactoring (c) 20201 by dbj@dbj.org -- https://dbj.org/license_dbj + + --- Original README From c83417043d6f36edc093100307c7f7f568710de8 Mon Sep 17 00:00:00 2001 From: dbjdbj Date: Thu, 22 Jul 2021 10:11:05 +0200 Subject: [PATCH 05/16] typo --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index d30356d..8372b3f 100644 --- a/README.md +++ b/README.md @@ -2,7 +2,7 @@ Unfortunately Strassen as coded here is not working. -> Refactoring (c) 20201 by dbj@dbj.org -- https://dbj.org/license_dbj +> Refactoring (c) 2021 by dbj@dbj.org -- https://dbj.org/license_dbj --- From 9540a29951ed4ae3504013208749d92e12f3202b Mon Sep 17 00:00:00 2001 From: dbjdbj Date: Thu, 22 Jul 2021 10:30:54 +0200 Subject: [PATCH 06/16] REFACTORED --- README.md | 10 ++++---- Strassen_Algorithm.c | 55 ++++++++++++++++++-------------------------- dbj_matrix.h | 25 ++++++++++++++++++++ 3 files changed, 54 insertions(+), 36 deletions(-) diff --git a/README.md b/README.md index 8372b3f..1c12dbd 100644 --- a/README.md +++ b/README.md @@ -1,11 +1,13 @@ -# 100% REFACTORED +

100%

+

 

-Unfortunately Strassen as coded here is not working. +#### REFACTORED SIMPLE, CORRECT AND WORKING -> Refactoring (c) 2021 by dbj@dbj.org -- https://dbj.org/license_dbj +

 

+> Refactoring (c) 2021 by dbj@dbj.org -- https://dbj.org/license_dbj ---- +

 

Original README diff --git a/Strassen_Algorithm.c b/Strassen_Algorithm.c index 085891f..ef45b9e 100644 --- a/Strassen_Algorithm.c +++ b/Strassen_Algorithm.c @@ -4,7 +4,7 @@ #include #include "dbj_matrix.h" -#define N 2 +#define SQUARE_MATRIX_SIDE 64 static inline Matrix *strassen(Matrix *dest, const Matrix *srcA, const Matrix *srcB, int length) { @@ -65,31 +65,10 @@ static inline Matrix *strassen(Matrix *dest, const Matrix *srcA, const Matrix *s return dest; } -static void matrix_print(const char *prompt, Matrix *mx) -{ - assert(prompt && mx); - printf("%s\n", prompt); - for (int i = 0; i < mx->rows; i++) - { - for (int j = 0; j < mx->cols; j++) - { - printf(" %5d ", MXY(mx, i, j)); - } - printf("\n"); - } -} - -static void inline matrix_foreach(Matrix *mx, Matrix *(*callback)(Matrix *, unsigned, unsigned)) -{ - for (int i = 0; i < mx->rows; i++) - { - for (int j = 0; j < mx->cols; j++) - { - callback(mx, i, j); - } - } -} - +/* +-------------------------------------------------------------------------------------------------------------------- +callbacks for for_each +*/ static inline Matrix *ordinal_as_val(Matrix *mx, unsigned i, unsigned j) { MXY(mx, i, j) = (i * mx->cols + j); @@ -108,12 +87,12 @@ static inline Matrix *zoro(Matrix *mx, unsigned i, unsigned j) */ int main(int argc, char *argv[]) { - const unsigned matrix_length = N; + const unsigned matrix_length = SQUARE_MATRIX_SIDE; /* Check if dimensions of matrix is the power of two or not */ if (ceil(log2(matrix_length)) != floor(log2(matrix_length))) { - printf("\n%s\tERROR: square matrix side must be a power of 2. And %d is not.", argv[0], matrix_length); + printf("\n%s\tERROR: square matrix side must be a power of 2. And current size:%d is not.", argv[0], matrix_length); return 0; } @@ -127,15 +106,27 @@ int main(int argc, char *argv[]) matrix_foreach(matrixC, zoro); matrix_foreach(matrixR, zoro); + printf("\n\nMultiplying int[%d][%d] * int[%d][%d]\n\n", SQUARE_MATRIX_SIDE, SQUARE_MATRIX_SIDE, SQUARE_MATRIX_SIDE, SQUARE_MATRIX_SIDE); + /* Matrix multiplication */ matrixC = strassen(matrixC, matrixA, matrixB, matrix_length); matrixR = matrix_ijk_matmul(matrixR, matrixA, matrixB); - matrix_print("Matrix A:", matrixA); - matrix_print("Matrix B:", matrixB); - matrix_print("Matrix C:", matrixC); - matrix_print("Matrix R:", matrixR); + if (SQUARE_MATRIX_SIDE < 9) + { + matrix_print("Matrix A:", matrixA); + matrix_print("Matrix B:", matrixB); + matrix_print("Matrix C:", matrixC); + matrix_print("Matrix R:", matrixR); + } + + printf("\n\nC is Strassen result and R is ijk_matmul result. They should be equal. "); + + if (!matrix_equal(matrixC, matrixR)) + printf("Unfortunately they are not.\n\n"); + else + printf("And indeed they are.\n\n"); return 0; } diff --git a/dbj_matrix.h b/dbj_matrix.h index 74a34e3..8f81f8d 100644 --- a/dbj_matrix.h +++ b/dbj_matrix.h @@ -101,4 +101,29 @@ static inline Matrix *matrix_new(uint32_t rows, uint32_t cols) return matrix; } +static void matrix_print(const char *prompt, Matrix *mx) +{ + assert(prompt && mx); + printf("%s\n", prompt); + for (int i = 0; i < mx->rows; i++) + { + for (int j = 0; j < mx->cols; j++) + { + printf(" %5d ", MXY(mx, i, j)); + } + printf("\n"); + } +} + +static void inline matrix_foreach(Matrix *mx, Matrix *(*callback)(Matrix *, unsigned, unsigned)) +{ + for (int i = 0; i < mx->rows; i++) + { + for (int j = 0; j < mx->cols; j++) + { + callback(mx, i, j); + } + } +} + #endif /* MATRIX_H_ */ From a7cb794ed98a2df8bb3e82599842722a5562984c Mon Sep 17 00:00:00 2001 From: dbjdbj Date: Thu, 22 Jul 2021 10:34:06 +0200 Subject: [PATCH 07/16] formating --- README.md | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/README.md b/README.md index 1c12dbd..6739bbb 100644 --- a/README.md +++ b/README.md @@ -1,17 +1,20 @@

100%

-

 

#### REFACTORED SIMPLE, CORRECT AND WORKING -

 

- +  +  > Refactoring (c) 2021 by dbj@dbj.org -- https://dbj.org/license_dbj - -

 

+  +  Original README +  +  ### Strassen Algorithm +  +  Implementation of the Strassen Algorithm in C From 5d147e6084de44dc20c44719c0be46b4b05a7f36 Mon Sep 17 00:00:00 2001 From: dbjdbj Date: Thu, 22 Jul 2021 10:37:29 +0200 Subject: [PATCH 08/16] formated --- README.md | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/README.md b/README.md index 6739bbb..acfed2c 100644 --- a/README.md +++ b/README.md @@ -1,20 +1,17 @@

100%

-#### REFACTORED SIMPLE, CORRECT AND WORKING +#### REFACTORED, SIMPLE, CORRECT AND WORKING -  -  > Refactoring (c) 2021 by dbj@dbj.org -- https://dbj.org/license_dbj + +Built on WIN10 using `TDM-GCC-64`. Clang or GCC are required. +     Original README -  -  ### Strassen Algorithm -  -  Implementation of the Strassen Algorithm in C From 961f4021fe14814ac2195dc74bae60497873bb3c Mon Sep 17 00:00:00 2001 From: dbjdbj Date: Thu, 22 Jul 2021 10:40:32 +0200 Subject: [PATCH 09/16] updated --- README.md | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/README.md b/README.md index acfed2c..a148a8c 100644 --- a/README.md +++ b/README.md @@ -1,8 +1,7 @@ -

100%

- -#### REFACTORED, SIMPLE, CORRECT AND WORKING - -> Refactoring (c) 2021 by dbj@dbj.org -- https://dbj.org/license_dbj +## 100% REFACTORED, SIMPLE, CORRECT AND WORKING +``` +Refactoring (c) 2021 by dbj@dbj.org -- https://dbj.org/license_dbj +``` Built on WIN10 using `TDM-GCC-64`. Clang or GCC are required. @@ -13,7 +12,7 @@ Original README ### Strassen Algorithm -Implementation of the Strassen Algorithm in C +Implementation of the Strassen Algorithm in C, by Ken Dai aka [`MetalheadKen`](https://github.com/MetalheadKen) The Strassen algorithm, is an algorithm for matrix multiplication. It is faster than the standard matrix multiplication algorithm, but would be slower than the From ef4b122e398a469f05e32906047855414b57b277 Mon Sep 17 00:00:00 2001 From: dbjdbj Date: Thu, 22 Jul 2021 13:48:20 +0200 Subject: [PATCH 10/16] FOR macro introduced --- Strassen_Algorithm.c | 8 ++++---- dbj_matrix.h | 28 ++++++++++++++-------------- 2 files changed, 18 insertions(+), 18 deletions(-) diff --git a/Strassen_Algorithm.c b/Strassen_Algorithm.c index ef45b9e..d7ea061 100644 --- a/Strassen_Algorithm.c +++ b/Strassen_Algorithm.c @@ -18,9 +18,9 @@ static inline Matrix *strassen(Matrix *dest, const Matrix *srcA, const Matrix *s *b11 = matrix_new(len, len), *b12 = matrix_new(len, len), *b21 = matrix_new(len, len), *b22 = matrix_new(len, len), *c11 = matrix_new(len, len), *c12 = matrix_new(len, len), *c21 = matrix_new(len, len), *c22 = matrix_new(len, len), *m1 = matrix_new(len, len), *m2 = matrix_new(len, len), *m3 = matrix_new(len, len), *m4 = matrix_new(len, len), *m5 = matrix_new(len, len), *m6 = matrix_new(len, len), *m7 = matrix_new(len, len), *temp1 = matrix_new(len, len), *temp2 = matrix_new(len, len); /* Divide matrix into four parts */ - for (int i = 0; i < len; i++) + FOR(i, len) { - for (int j = 0; j < len; j++) + FOR(j, len) { MXY(a11, i, j) = MXY(srcA, i, j); MXY(a12, i, j) = MXY(srcA, i, j + len); @@ -51,9 +51,9 @@ static inline Matrix *strassen(Matrix *dest, const Matrix *srcA, const Matrix *s matrix_addition(c22, matrix_subtraction(temp1, m1, m2), matrix_addition(temp2, m3, m6)); /* Store the answer of matrix multiplication */ - for (int i = 0; i < len; i++) + FOR(i, len) { - for (int j = 0; j < len; j++) + FOR(j, len) { MXY(dest, i, j) = MXY(c11, i, j); MXY(dest, i, j + len) = MXY(c12, i, j); diff --git a/dbj_matrix.h b/dbj_matrix.h index 8f81f8d..16f0cfc 100644 --- a/dbj_matrix.h +++ b/dbj_matrix.h @@ -40,10 +40,10 @@ static inline void matrix_free(Matrix **ptr) static inline Matrix *matrix_addition(Matrix *res, const Matrix *a, const Matrix *b) { assert(res && a && b); - for (uint32_t i = 0; i < res->rows; i++) - for (uint32_t j = 0; j < res->cols; j++) - // res->values[i][j] = a->values[i][j] + b->values[i][j]; - MXY(res, i, j) = MXY(a, i, j) + MXY(b, i, j); + FOR(i, res->rows) + FOR(j, res->cols) + // res->values[i][j] = a->values[i][j] + b->values[i][j]; + MXY(res, i, j) = MXY(a, i, j) + MXY(b, i, j); return res; } @@ -51,9 +51,9 @@ static inline Matrix *matrix_addition(Matrix *res, const Matrix *a, const Matrix static inline Matrix *matrix_subtraction(Matrix *res, const Matrix *a, const Matrix *b) { assert(res && a && b); - for (uint32_t i = 0; i < res->rows; i++) - for (uint32_t j = 0; j < res->cols; j++) - MXY(res, i, j) = MXY(a, i, j) - MXY(b, i, j); + FOR(i, res->rows) + FOR(j, res->cols) + MXY(res, i, j) = MXY(a, i, j) - MXY(b, i, j); return res; } @@ -78,12 +78,12 @@ static inline int matrix_equal(const Matrix *a, const Matrix *b) static Matrix *matrix_ijk_matmul(Matrix *C, const Matrix *A, const Matrix *B) { assert(A && B && C); - for (size_t i = 0; i < A->rows; ++i) + FOR(i, A->rows) { - for (size_t j = 0; j < B->cols; ++j) + FOR(j, B->cols) { MXY(C, i, j) = 0; - for (size_t k = 0; k < A->cols; ++k) + FOR(k, A->cols) { MXY(C, i, j) += MXY(A, i, k) * MXY(B, k, j); } @@ -105,9 +105,9 @@ static void matrix_print(const char *prompt, Matrix *mx) { assert(prompt && mx); printf("%s\n", prompt); - for (int i = 0; i < mx->rows; i++) + FOR(i, mx->rows) { - for (int j = 0; j < mx->cols; j++) + FOR(j, mx->cols) { printf(" %5d ", MXY(mx, i, j)); } @@ -117,9 +117,9 @@ static void matrix_print(const char *prompt, Matrix *mx) static void inline matrix_foreach(Matrix *mx, Matrix *(*callback)(Matrix *, unsigned, unsigned)) { - for (int i = 0; i < mx->rows; i++) + FOR(i, mx->rows) { - for (int j = 0; j < mx->cols; j++) + FOR(j, mx->cols) { callback(mx, i, j); } From a08863694d5970863d8d8eeb5d14d3eca1368ca5 Mon Sep 17 00:00:00 2001 From: dbjdbj Date: Thu, 22 Jul 2021 13:52:53 +0200 Subject: [PATCH 11/16] ISPOW2 --- Strassen_Algorithm.c | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/Strassen_Algorithm.c b/Strassen_Algorithm.c index d7ea061..f159760 100644 --- a/Strassen_Algorithm.c +++ b/Strassen_Algorithm.c @@ -67,7 +67,7 @@ static inline Matrix *strassen(Matrix *dest, const Matrix *srcA, const Matrix *s /* -------------------------------------------------------------------------------------------------------------------- -callbacks for for_each +callbacks for matrix_for_each */ static inline Matrix *ordinal_as_val(Matrix *mx, unsigned i, unsigned j) { @@ -85,12 +85,15 @@ static inline Matrix *zoro(Matrix *mx, unsigned i, unsigned j) /* -------------------------------------------------------------------------------------------------------------------- */ +/* Check if value is the power of two or not */ +#define ISPOW2(V_) (ceil(log2(V_)) == floor(log2(V_))) + int main(int argc, char *argv[]) { const unsigned matrix_length = SQUARE_MATRIX_SIDE; /* Check if dimensions of matrix is the power of two or not */ - if (ceil(log2(matrix_length)) != floor(log2(matrix_length))) + if (!ISPOW2(matrix_length)) { printf("\n%s\tERROR: square matrix side must be a power of 2. And current size:%d is not.", argv[0], matrix_length); return 0; From c8ef9b951a8c98824d0f76e2f0ea93e16c8f5184 Mon Sep 17 00:00:00 2001 From: dbjdbj Date: Thu, 22 Jul 2021 13:57:15 +0200 Subject: [PATCH 12/16] slightly better info --- Strassen_Algorithm.c | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/Strassen_Algorithm.c b/Strassen_Algorithm.c index f159760..2e7f4c5 100644 --- a/Strassen_Algorithm.c +++ b/Strassen_Algorithm.c @@ -90,29 +90,29 @@ static inline Matrix *zoro(Matrix *mx, unsigned i, unsigned j) int main(int argc, char *argv[]) { - const unsigned matrix_length = SQUARE_MATRIX_SIDE; + const unsigned matrix_side_size = SQUARE_MATRIX_SIDE; /* Check if dimensions of matrix is the power of two or not */ - if (!ISPOW2(matrix_length)) + if (!ISPOW2(matrix_side_size)) { - printf("\n%s\tERROR: square matrix side must be a power of 2. And current size:%d is not.", argv[0], matrix_length); + printf("\n%s\tERROR: square matrix side must be a power of 2. And current size:%d is not.", argv[0], matrix_side_size); return 0; } - matrix_autofree Matrix *matrixA = matrix_new(matrix_length, matrix_length); - matrix_autofree Matrix *matrixB = matrix_new(matrix_length, matrix_length); - matrix_autofree Matrix *matrixC = matrix_new(matrix_length, matrix_length); - matrix_autofree Matrix *matrixR = matrix_new(matrix_length, matrix_length); + matrix_autofree Matrix *matrixA = matrix_new(matrix_side_size, matrix_side_size); + matrix_autofree Matrix *matrixB = matrix_new(matrix_side_size, matrix_side_size); + matrix_autofree Matrix *matrixC = matrix_new(matrix_side_size, matrix_side_size); + matrix_autofree Matrix *matrixR = matrix_new(matrix_side_size, matrix_side_size); matrix_foreach(matrixA, ordinal_as_val); matrix_foreach(matrixB, ordinal_as_val); matrix_foreach(matrixC, zoro); matrix_foreach(matrixR, zoro); - printf("\n\nMultiplying int[%d][%d] * int[%d][%d]\n\n", SQUARE_MATRIX_SIDE, SQUARE_MATRIX_SIDE, SQUARE_MATRIX_SIDE, SQUARE_MATRIX_SIDE); + printf("\n\nMatrix multiplication: C = A * B\nMultiplying int[%d][%d] * int[%d][%d]\n R is the result of ijk_matmul(A,B)\n\n", SQUARE_MATRIX_SIDE, SQUARE_MATRIX_SIDE, SQUARE_MATRIX_SIDE, SQUARE_MATRIX_SIDE); /* Matrix multiplication */ - matrixC = strassen(matrixC, matrixA, matrixB, matrix_length); + matrixC = strassen(matrixC, matrixA, matrixB, matrix_side_size); matrixR = matrix_ijk_matmul(matrixR, matrixA, matrixB); From 7d44ab4a485544af8e0cbc3a55ab8e46841a70fc Mon Sep 17 00:00:00 2001 From: dbjdbj Date: Thu, 22 Jul 2021 21:25:52 +0200 Subject: [PATCH 13/16] rectangular not square matrices also work ok --- Strassen_Algorithm.c | 115 ++++++++++++++++++++++++++++++++----------- 1 file changed, 85 insertions(+), 30 deletions(-) diff --git a/Strassen_Algorithm.c b/Strassen_Algorithm.c index 2e7f4c5..71046d6 100644 --- a/Strassen_Algorithm.c +++ b/Strassen_Algorithm.c @@ -4,7 +4,11 @@ #include #include "dbj_matrix.h" -#define SQUARE_MATRIX_SIDE 64 +#define MATRIX_ROWS 2 +#define MATRIX_COLS 4 + +/* Check if value is the power of two or not */ +#define ISPOW2(V_) (ceil(log2(V_)) == floor(log2(V_))) static inline Matrix *strassen(Matrix *dest, const Matrix *srcA, const Matrix *srcB, int length) { @@ -84,52 +88,103 @@ static inline Matrix *zoro(Matrix *mx, unsigned i, unsigned j) /* -------------------------------------------------------------------------------------------------------------------- +nano app framework so that it is maleable to using it with UBENCH.H +*/ +typedef struct APP_DATA APP_DATA; +void app_start(const int, const char *const[]); +void app_end(); + +/* +-------------------------------------------------------------------------------------------------------------------- +nano app framework so that it is maleable to using it with UBENCH.H */ -/* Check if value is the power of two or not */ -#define ISPOW2(V_) (ceil(log2(V_)) == floor(log2(V_))) -int main(int argc, char *argv[]) +static struct APP_DATA { - const unsigned matrix_side_size = SQUARE_MATRIX_SIDE; + Matrix *matrixA; + Matrix *matrixB; + Matrix *matrixC; + Matrix *matrixR; +} *app_data = 0; - /* Check if dimensions of matrix is the power of two or not */ - if (!ISPOW2(matrix_side_size)) +void app_start(const int argc, const char *const argv[]) +{ + /* Check if dimensions of matrix are the power of two */ + if (!ISPOW2(MATRIX_ROWS)) { - printf("\n%s\tERROR: square matrix side must be a power of 2. And current size:%d is not.", argv[0], matrix_side_size); - return 0; + printf("\n%s\tERROR: square matrix side must be a power of 2. And current rows num:%d is not.", argv[0], MATRIX_ROWS); + exit(EXIT_FAILURE); } - matrix_autofree Matrix *matrixA = matrix_new(matrix_side_size, matrix_side_size); - matrix_autofree Matrix *matrixB = matrix_new(matrix_side_size, matrix_side_size); - matrix_autofree Matrix *matrixC = matrix_new(matrix_side_size, matrix_side_size); - matrix_autofree Matrix *matrixR = matrix_new(matrix_side_size, matrix_side_size); - - matrix_foreach(matrixA, ordinal_as_val); - matrix_foreach(matrixB, ordinal_as_val); - matrix_foreach(matrixC, zoro); - matrix_foreach(matrixR, zoro); - - printf("\n\nMatrix multiplication: C = A * B\nMultiplying int[%d][%d] * int[%d][%d]\n R is the result of ijk_matmul(A,B)\n\n", SQUARE_MATRIX_SIDE, SQUARE_MATRIX_SIDE, SQUARE_MATRIX_SIDE, SQUARE_MATRIX_SIDE); - - /* Matrix multiplication */ - matrixC = strassen(matrixC, matrixA, matrixB, matrix_side_size); + if (!ISPOW2(MATRIX_COLS)) + { + printf("\n%s\tERROR: square matrix side must be a power of 2. And current cols num:%d is not.", argv[0], MATRIX_COLS); + exit(EXIT_FAILURE); + } - matrixR = matrix_ijk_matmul(matrixR, matrixA, matrixB); + app_data = calloc(1, sizeof(APP_DATA)); + assert(app_data); + app_data->matrixA = matrix_new(MATRIX_ROWS, MATRIX_COLS); + app_data->matrixB = matrix_new(MATRIX_ROWS, MATRIX_COLS); + app_data->matrixC = matrix_new(MATRIX_ROWS, MATRIX_COLS); + app_data->matrixR = matrix_new(MATRIX_ROWS, MATRIX_COLS); + + matrix_foreach(app_data->matrixA, ordinal_as_val); + matrix_foreach(app_data->matrixB, ordinal_as_val); + matrix_foreach(app_data->matrixC, zoro); + matrix_foreach(app_data->matrixR, zoro); +} +void app_end() +{ - if (SQUARE_MATRIX_SIDE < 9) + printf("\n\nMatrix multiplication\n\nC = strassen(A, B)\nR = ijk_matmul(A,B)\n\n"); + if (MATRIX_ROWS < 9) { - matrix_print("Matrix A:", matrixA); - matrix_print("Matrix B:", matrixB); - matrix_print("Matrix C:", matrixC); - matrix_print("Matrix R:", matrixR); + printf("\n\nAfter multiplications:\n"); + matrix_print("Matrix A:", app_data->matrixA); + matrix_print("Matrix B:", app_data->matrixB); + matrix_print("Matrix C:", app_data->matrixC); + matrix_print("Matrix R:", app_data->matrixR); } printf("\n\nC is Strassen result and R is ijk_matmul result. They should be equal. "); - if (!matrix_equal(matrixC, matrixR)) + if (!matrix_equal(app_data->matrixC, app_data->matrixR)) printf("Unfortunately they are not.\n\n"); else printf("And indeed they are.\n\n"); + matrix_free(&app_data->matrixA); + matrix_free(&app_data->matrixB); + matrix_free(&app_data->matrixC); + matrix_free(&app_data->matrixR); + + free(app_data); + app_data = 0; +} +/* +-------------------------------------------------------------------------------------------------------------------- +*/ +static void do_strassen() +{ + /* Matrix multiplication */ + app_data->matrixC = strassen(app_data->matrixC, app_data->matrixA, app_data->matrixB, MATRIX_ROWS); +} + +static void do_ijk() +{ + app_data->matrixR = matrix_ijk_matmul(app_data->matrixR, app_data->matrixA, app_data->matrixB); +} + +/* +-------------------------------------------------------------------------------------------------------------------- +*/ + +int main(int argc, const char *const argv[]) +{ + app_start(argc, argv); + do_strassen(); + do_ijk(); + app_end(); return 0; } From 3f6d9d86e6a0c0c816be013293ead7a709f6ce6c Mon Sep 17 00:00:00 2001 From: dbjdbj Date: Thu, 22 Jul 2021 21:27:28 +0200 Subject: [PATCH 14/16] update --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index a148a8c..4b29b98 100644 --- a/README.md +++ b/README.md @@ -3,7 +3,7 @@ Refactoring (c) 2021 by dbj@dbj.org -- https://dbj.org/license_dbj ``` -Built on WIN10 using `TDM-GCC-64`. Clang or GCC are required. +Clang or GCC are required. Built on WIN10 using `TDM-GCC-64` or `clang 11.0.0`.     From 2198da72a7bbce51986d4ed3965430b8f6def6a5 Mon Sep 17 00:00:00 2001 From: dbjdbj Date: Thu, 22 Jul 2021 21:29:26 +0200 Subject: [PATCH 15/16] update --- Strassen_Algorithm.c | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/Strassen_Algorithm.c b/Strassen_Algorithm.c index 71046d6..a48982d 100644 --- a/Strassen_Algorithm.c +++ b/Strassen_Algorithm.c @@ -112,13 +112,17 @@ void app_start(const int argc, const char *const argv[]) /* Check if dimensions of matrix are the power of two */ if (!ISPOW2(MATRIX_ROWS)) { - printf("\n%s\tERROR: square matrix side must be a power of 2. And current rows num:%d is not.", argv[0], MATRIX_ROWS); + printf("\n%s\tERROR: Each matrix side must be a power of 2. " + "And current rows num:%d is not.", + argv[0], MATRIX_ROWS); exit(EXIT_FAILURE); } if (!ISPOW2(MATRIX_COLS)) { - printf("\n%s\tERROR: square matrix side must be a power of 2. And current cols num:%d is not.", argv[0], MATRIX_COLS); + printf("\n%s\tERROR: Each matrix side must be a power of 2. " + "And current cols num:%d is not.", + argv[0], MATRIX_COLS); exit(EXIT_FAILURE); } From 30b975e8673b535f5e1d042dc66b65b2894400af Mon Sep 17 00:00:00 2001 From: dbjdbj Date: Fri, 23 Jul 2021 02:24:04 +0200 Subject: [PATCH 16/16] update --- Strassen_Algorithm.c | 1 - 1 file changed, 1 deletion(-) diff --git a/Strassen_Algorithm.c b/Strassen_Algorithm.c index a48982d..bcdd918 100644 --- a/Strassen_Algorithm.c +++ b/Strassen_Algorithm.c @@ -96,7 +96,6 @@ void app_end(); /* -------------------------------------------------------------------------------------------------------------------- -nano app framework so that it is maleable to using it with UBENCH.H */ static struct APP_DATA