Modular function fix and installing OpenBLAS

Author: audictheprogrammer

src/data/Matrix_A.txt

@@ -1,7 +1,6 @@
-255 0.356913 
-511 3.051232 
-256 0.353928 
-512 3.022207 
-256 0.366828 
-512 3.066367 
-1024 27.525820 
+256 0.406682
+512 3.161513
+1024 28.813085
+256 0.360843 
+512 3.050652 
+1024 27.503507 

src/data/Matrix_B.txt

@@ -1,3 +1,3 @@
-256 0.447807 
-512 2.914416 
-1024 25.676863 
+256 0.320855 
+512 2.779970 
+1024 25.977462 

src/data/Matrix_C.txt

@@ -1,3 +1,3 @@
-256 0.393323 
-512 2.573454 
-1024 24.374175 
+256 0.353629 
+512 2.799502 
+1024 24.423494 

src/data/Matrix_D.txt

@@ -1,3 +1,3 @@
-256 0.308907 
-512 2.618251 
-1024 24.035276 
+256 0.312455 
+512 2.680338 
+1024 24.309112 

src/data/Matrix_E.txt

@@ -1,3 +1,3 @@
-256 0.276151 
-512 2.137076 
-1024 20.249978 
+256 0.236347 
+512 2.194848 
+1024 20.595512 

src/data/Matrix_F.txt

@@ -159,7 +159,7 @@ void benchmark_modulos(double p, double u, double u_overline, double u_b){
     /* Benchmarking different modulos.
     The most efficient one is IKJ.
     */
-    int m = 1;  // Executes m times each algo
+    int m = 2;  // Executes m times each algo
     for (int i=8; i<11; i++){
         int n = (int) pow(2, i);
         double sum_mod_naive = 0;

src/data/Matrix_G.txt

@@ -0,0 +1,18 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <cblas.h>
+
+
+int main(){
+    int i=0;
+    double A[6] = {1.0,2.0,1.0,-3.0,4.0,-1.0};
+    double B[6] = {1.0,2.0,1.0,-3.0,4.0,-1.0};
+    double C[9] = {.5,.5,.5,.5,.5,.5,.5,.5,.5};
+    cblas_dgemm(CblasColMajor, CblasNoTrans, CblasTrans,3,3,2,1,A, 3, B, 3,2,C,3);
+
+    for(i=0; i<9; i++){
+        printf("%lf ", C[i]);
+    }
+    printf("\n");
+    return 0;
+}

src/data/benchmark_modulo_Barrett.txt

@@ -24,9 +24,9 @@ int main(int argc, char** argv){
     const int TEST1 = 0;
     const int TEST2 = 0;
     const int TEST3 = 0;
-    const int TEST4 = 0;
+    const int TEST4 = 1;
     const int TEST5 = 0;
-    const int TEST6 = 1;
+    // const int TEST6 = 0;
     // const int TEST7 = 1;
 
 
@@ -35,11 +35,11 @@ int main(int argc, char** argv){
         double p = pow(2, 26) - 5;
 
         // Precomputed constants for Modular functions
-        fesetround(FE_DOWNWARD);
+        fesetround(FE_TONEAREST);
         double u = 1.0 / p;  // Constant for SIMD
         fesetround(FE_UPWARD);
         double u_overline = 1.0 / p;  // Constant for SIMD2 and SIMD3
-        fesetround(FE_DOWNWARD);
+        fesetround(FE_TONEAREST);
         u_int32_t u_b = (int) (pow(2, 56) / p);  // Constant for Barrett
 
         double a = pow(2, 26);
@@ -64,11 +64,11 @@ int main(int argc, char** argv){
         double p = pow(2, 26) - 5;
 
         // Precomputed constants for Modular functions
-        fesetround(FE_DOWNWARD);
+        fesetround(FE_TONEAREST);
         double u = 1.0 / p;  // Constant for SIMD
         fesetround(FE_UPWARD);
         double u_overline = 1.0 / p;  // Constant for SIMD2 and SIMD3
-        fesetround(FE_DOWNWARD);
+        fesetround(FE_TONEAREST);
         u_int32_t u_b = (int) (pow(2, 56) / p);  // Constant for Barrett
 
         double a = (p-1) * (p-1) * 32;
@@ -93,11 +93,11 @@ int main(int argc, char** argv){
         double p = pow(2, 26) - 5;
 
         // Precomputed constants for Modular functions
-        fesetround(FE_DOWNWARD);
+        fesetround(FE_TONEAREST);
         double u = 1.0 / p;  // Constant for SIMD
         fesetround(FE_UPWARD);
         double u_overline = 1.0 / p;  // Constant for SIMD2 and SIMD3
-        fesetround(FE_DOWNWARD);
+        fesetround(FE_TONEAREST);
         u_int32_t u_b = (int) (pow(2, 56) / p);  // Constant for Barrett
 
         double a = (p-1) * (p-1) * 16;  // a > 2^(25+25+4) = 2^54
@@ -123,16 +123,16 @@ int main(int argc, char** argv){
         double p = pow(2, 26) - 5;
 
         // Precomputed constants for Modular functions
-        fesetround(FE_DOWNWARD);
+        fesetround(FE_TONEAREST);
         double u = 1.0 / p;  // Constant for SIMD
         fesetround(FE_UPWARD);
         double u_overline = 1.0 / p;  // Constant for SIMD2 and SIMD3
-        fesetround(FE_DOWNWARD);
         u_int32_t u_b = (int) (pow(2, 56) / p);  // Constant for Barrett
+        fesetround(FE_TONEAREST);
 
-        int n = 200;
+        int n = 1024;
 
-        for (int i=0; i<10; i++){
+        for (int i=0; i<1; i++){
 
             double**A = random_matrix(n, p);
             double**B = random_matrix(n, p);
@@ -163,6 +163,16 @@ int main(int argc, char** argv){
             int nb3 = equals_matrix(C, F, n);
             int nb4 = equals_matrix(C, G, n);
 
+            delete_matrix(&A, n);
+            delete_matrix(&B, n);
+            delete_matrix(&C, n);
+            delete_matrix(&D, n);
+            delete_matrix(&E, n);
+            delete_matrix(&F, n);
+            delete_matrix(&G, n);
+
+
+            printf("i=%d \n", i);
             assert(nb1==1);
             assert(nb2==1);
             assert(nb3==1);
@@ -177,49 +187,24 @@ int main(int argc, char** argv){
         /* Testing the equivalence between
         if (d < 0) return d+p; else return d;
         And
-        return d + p()
+        return d + (-(d<0) & (int)p)
+        And
+        return d + p * (d<0)
         */
 
         double d1 = -2;
         double d2 = 2;
         double p = 100.0;
 
-        double a = 1024.5;
-        double b = 0b1111111111111111111111111111111111111111111111111111111111111111;
-        double temp = b && a;
-
-        // 0b 1111 1111 ... 1111 for double = ...
-
-        double res1 = d1 + ((-(d1<0)) && p);
-        double res2 = d2 + ((-(d2<0)) && p);
+        double res1 = d1 + ((-(d1<0)) & (u_int64_t)p);
+        double res2 = d2 + ((-(d2<0)) & (u_int64_t)p);
+        double res3 = d1 + p * (d1<0);
+        double res4 = d2 + p * (d2<0);
 
         printf("res1 = %f \n", res1);
         printf("res2 = %f \n", res2);
-        printf("b = %f \n", b);
-        printf("temp = %f \n", temp);
-    }
-
-    if (TEST6){
-        /* Testing FE_DOWNWARD and FE_UPWARD */
-        double p = pow(2, 26) - 5;
-        // double u1 = (p-1) * (p-1) * (p-1);
-        double u1 = -(p-1) * (p-1) * (p-1) * 64;
-        fesetround(FE_DOWNWARD);
-        double u2 = -(p-1) * (p-1) * (p-1) * 64;
-        fesetround(FE_UPWARD);
-        double u3 = -(p-1) * (p-1) * (p-1) * 64;
-        fesetround(FE_TOWARDZERO);
-        double u4 = -(p-1) * (p-1) * (p-1) * 64;
-        fesetround(FE_TONEAREST);
-        double u5 = -(p-1) * (p-1) * (p-1) * 64;
-
-        printf("u1=%f \n", u1);
-        printf("u2=%f \n", u2);
-        printf("u3=%f \n", u3);
-        printf("u4=%f \n", u3);
-        printf("u5=%f \n", u3);
-
-
+        printf("res3 = %f \n", res3);
+        printf("res4 = %f \n", res4);
     }
 
     return 0;

src/data/benchmark_modulo_SIMD1.txt

@@ -1,5 +1,5 @@
 # Targets
-all: main main_test
+all: main main_test main_blas
 
 # Headers
 matrix.o: matrix.c
@@ -12,7 +12,9 @@ main: main.c matrix.o matrix_mul.o
 	gcc -o main main.c matrix.o matrix_mul.o -lm
 main_test: main_test.c matrix.o matrix_mul.o
 	gcc -o main_test main_test.c matrix.o matrix_mul.o -lm
+main_blas: main_blas.c
+	gcc -o main_blas main_blas.c -lopenblas
 
 clean:
 	rm *.o
-	rm -rf main main_test
+	rm -rf main main_test main_blas

src/data/benchmark_modulo_SIMD2.txt

@@ -40,11 +40,6 @@ double** random_matrix(int size, double p){
             // int rand_int = rand();
             // mat[i][j] = modulo_SIMD1(rand_int, p);
             mat[i][j] = (double) (rand() % ((int) p));
-
-            // mat[i][j] = (double)rand()/(double)(RAND_MAX/p);
-
-
-            // mat[i][j] = rand() % mod;
         }
     }
 

src/data/benchmark_modulo_SIMD3.txt

@@ -12,8 +12,18 @@ double modulo_SIMD1(double a, double p, double u){
     double b = a * u;
     double c = (double)(int)b;
     double d = a - c * p;
-    if (d >= p) return d-p;
-    if (d < 0) return d+p;
+    // if (d >= p) return d-p;
+    if (d >= p) {
+        printf("SIMD1: if1 \n");
+        return d-p;
+    }
+    // printf("SIMD1: if1 \n");
+    if (d < 0){
+        printf("SIMD2: if2 \n");
+        return d+p;
+    }
+
+    // if (d < 0) return d+p;
     return d;
 }
 
@@ -24,7 +34,11 @@ double modulo_SIMD2(double a, double p, double u){
     double b = a * u;
     double c = (double)(int)b;
     double d = a - c * p;
-    if (d < 0) return d+p;
+    // if (d < 0) return d+p;
+    if (d < 0){
+        printf("d < 0 \n");
+        return d+p;
+    }
     return d;
 }
 
@@ -35,7 +49,8 @@ double modulo_SIMD3(double a, double p, double u){
     double b = a * u;
     double c = (double)(int)b;
     double d = a - c * p;
-    return d + (-(d<0) && p);
+    // return d + (-(d<0) & (int64_t) p);
+    return d + p * (d<0);
 }
 
 u_int32_t modulo_Barrett(u_int64_t a, u_int32_t p, u_int32_t u){

src/data/benchmark_modulo_naive.txt

@@ -19,6 +19,6 @@ replot "../data/benchmark_order_ikj.txt" using 1:2 title "No Mod" with linespoin
 replot "../data/benchmark_modulo_SIMD1.txt" using 1:2 title "SIMD1" with linespoints ls 3
 replot "../data/benchmark_modulo_SIMD2.txt" using 1:2 title "SIMD2" with linespoints ls 4
 replot "../data/benchmark_modulo_SIMD3.txt" using 1:2 title "SIMD3" with linespoints ls 5
-replot "../data/benchmark_modulo_Barrett.txt" using 1:2 title "Bararett" with linespoints ls 6
+# replot "../data/benchmark_modulo_Barrett.txt" using 1:2 title "Barrett" with linespoints ls 6
 
 # Graphe Order of loops