Testing CBLAS_dgemm

Author: audictheprogrammer

src/data/Matrix_A_1D.txt

@@ -0,0 +1,4 @@
+3
+[ 1, 2, 1 ]
+[ 3, 4, 1 ]
+[ 0, 0, 1 ]

src/data/Matrix_B_1D.txt

@@ -0,0 +1,4 @@
+3
+[ 0, 1, 3 ]
+[ 0, 0, 0 ]
+[ 3, 0, 3]

src/main_test.c

@@ -25,10 +25,10 @@ int main(int argc, char** argv){
     const int TEST1 = 0;
     const int TEST2 = 0;
     const int TEST3 = 0;
-    const int TEST4 = 1;
+    const int TEST4 = 0;
     const int TEST5 = 0;
-    const int TEST6 = 1;
-    // const int TEST7 = 1;
+    const int TEST6 = 0;
+    const int TEST7 = 1;
 
 
     if (TEST1){
@@ -215,8 +215,8 @@ int main(int argc, char** argv){
         int bitsize_p = 26;
         int n = 2;
 
-        int get = get_blocksize(bitsize_p, n);
-        printf("get = %d \n", get);
+        int blocksize = get_blocksize(bitsize_p, n);
+        printf("blocksize = %d \n", blocksize);
         double** A = random_matrix_2D(n, p);
         double** B = random_matrix_2D(n, p);
 
@@ -225,8 +225,7 @@ int main(int argc, char** argv){
         double* C = zero_matrix_1D(n*n);
         double** D = zero_matrix_2D(n);
 
-
-        cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,n,n,n, 1, A_1D,n, B_1D,n, 1, C,n);
+        cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,n,n,n, 1, A_1D, n-1, B_1D,n, 1, C,n);
         mp_ikj(A, B, D, n);
 
         print_matrix_2D(A, n);
@@ -249,6 +248,39 @@ int main(int argc, char** argv){
 
     }
 
+    if (TEST7){
+        // Testing cblas_dgemm's lda ldb and ldc parameters.
+        int n;
+
+        double** A = read_matrix("data/Matrix_A_1D.txt", &n);
+        double** B = read_matrix("data/Matrix_B_1D.txt", &n);
+        double* A_1D = convert_2D_to_1D(A, n);
+        double* B_1D = convert_2D_to_1D(A, n);
+        double* C = zero_matrix_1D(n*n);
+        double* D = zero_matrix_1D(n*n);
+        double** E = zero_matrix_2D(n);
+
+        cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,n,n,n, 1, A_1D, 900, B_1D, 3, 1, C,3);
+        cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,n,n,n, 1, A_1D, 6, B_1D, 3, 1, D,3);
+        mp_ijk(A, B, E, n);
+
+        print_matrix_1D(C, n);
+        printf("\n");
+        print_matrix_1D(D, n);
+        printf("\n");
+        print_matrix_2D(E, n);
+
+
+        delete_matrix_2D(&A, n);
+        delete_matrix_2D(&B, n);
+        delete_matrix_2D(&E, n);
+
+        delete_matrix_1D(&A_1D, n);
+        delete_matrix_1D(&B_1D, n);
+        delete_matrix_1D(&C, n);
+        delete_matrix_1D(&D, n);
+
+    }
 
 
 

src/matrix.c

@@ -255,7 +255,6 @@ double** read_matrix(char* filename, int* n){
 }
 
 int equals_matrix_2D_2D(double** A, double** B, int n){
-
     /* Check if A equals B. The return values are:
     0: different
     1: equals
@@ -281,7 +280,6 @@ int equals_matrix_2D_2D(double** A, double** B, int n){
 }
 
 int equals_matrix_2D_1D(double** A, double* B, int n){
-
     /* Check if A equals B. The return values are:
     0: different
     1: equals
@@ -301,6 +299,14 @@ int equals_matrix_file(char* filename1, char* filename2){
     int n2;
     double** mat1 = read_matrix(filename1, &n1);
     double** mat2 = read_matrix(filename2, &n2);
+    if (n1 != n2){
+        delete_matrix_2D(&mat1, n1);
+        delete_matrix_2D(&mat2, n2);
+        return 0;
+    }
+    int res = equals_matrix_2D_2D(mat1, mat2, n1);
+    delete_matrix_2D(&mat1, n1);
+    delete_matrix_2D(&mat2, n2);
 
-    return equals_matrix_2D_2D(mat1, mat2, n1);
+    return res;
 }

src/matrix.h

@@ -9,19 +9,27 @@
 // Basic operations
 double** zero_matrix_2D(int n);
 double* zero_matrix_1D(int n);
+
 double** random_matrix_2D(int n, double p);
 double* random_matrix_1D(int n, double p);
+
 double* convert_2D_to_1D(double** mat, int n);
 double** convert_1D_to_2D(double* mat, int n);
+
 void delete_matrix_2D(double*** mat, int n);
 void delete_matrix_1D(double** mat, int n);
+
 double** transpose_matrix(double** mat, int n);
+
 void print_matrix_2D(double** mat, int n);
 void print_matrix_1D(double* mat, int n);
+
 void write_matrix(double** mat, int n, char* filename);
 double** read_matrix(char* filename, int* n);
+
 int equals_matrix_2D_2D(double** A, double** B, int n);
 int equals_matrix_2D_1D(double** A, double* B, int n);
+
 int equals_matrix_file(char* filename1, char* filename2);
 
 

src/matrix_mul.c

@@ -6,7 +6,7 @@ double modulo_naive(double a, double p){
 }
 
 double modulo_SIMD1(double a, double p, double u){
-    /* Function .1 from SIMD article for floats.
+    /* Function 3.1 from SIMD article for floats.
     Hypothesis: Rounding mode = nearest and p < 2^26.
     */
     double b = a * u;
@@ -71,7 +71,7 @@ u_int32_t modulo_Barrett(u_int64_t a, u_int32_t p, u_int32_t u){
 }
 
 void mp_naive(double** A, double** B, double** C, int n, double p){
-    // Assert C is a zero matrix
+    // Assert C is a zero matrix.
     for (int i=0; i<n; i++){
         for (int k=0; k<n; k++){
             for (int j=0; j<n; j++){
@@ -166,9 +166,10 @@ void mp_Barrett(double** A, double** B, double** C, int n, double p, u_int32_t u
 
 }
 
-void mp_block(double** A, double** B, double** C, int n, double p, int bitsize_p){
+void mp_block(double* A, double* B, double* C, int n, double p, int bitsize_p){
+    /* Compute the matrix product A*B using blocks and OpenBLAS.
+    */
 
-
 }
 
 

src/matrix_mul.h

@@ -19,7 +19,7 @@ void mp_SIMD1(double** A, double** B, double** C, int n, double p, double u);
 void mp_SIMD2(double** A, double** B, double** C, int n, double p, double u);
 void mp_SIMD3(double** A, double** B, double** C, int n, double p, double u);
 void mp_Barrett(double** A, double** B, double** C, int n, double p, u_int32_t u);
-void mp_block(double** A, double** B, double** C, int n, double p, int bitsize_p);
+void mp_block(double* A, double* B, double* C, int n, double p, int bitsize_p);
 
 void mp_ijk(double** A, double** B, double** C, int n);
 void mp_kij(double** A, double** B, double** C, int n);