From 19e04dbf6d39f059f517973bd3c6e2bde448cd93 Mon Sep 17 00:00:00 2001
From: Nicholas Noll <nbnoll@eml.cc>
Date: Thu, 30 Apr 2020 13:25:02 -0700
Subject: feat: added more level 2 functions

---
 sys/libmath/blas.c | 96 ++++++++++++++++++++++++++++++++++++++----------------
 1 file changed, 68 insertions(+), 28 deletions(-)

(limited to 'sys/libmath/blas.c')

diff --git a/sys/libmath/blas.c b/sys/libmath/blas.c
index a008ffb..b42a3b8 100644
--- a/sys/libmath/blas.c
+++ b/sys/libmath/blas.c
@@ -326,21 +326,50 @@ blas·gemv(int nrow, int ncol, double a, double *m, double *x, double b, double
 
 /* 
  * rank one addition
- * M = axy + M
+ * M = ax(y^T) + M
  */
 
-#define NITER 50
-#define NCOL  1000
-#define NROW  1000
+void
+blas·ger(int nrow, int ncol, double a, double *x, double *y, double *m)
+{
+    int i, j;
+
+    for (i = 0; i < nrow; i++) {
+        for (j = 0; j < ncol; j++) {
+            m[i+ncol*j] += a * x[i] * y[j];
+        }
+    }
+}
+
+/* 
+ * symmetric rank one addition
+ * M = ax(x^T) + M
+ */
+
+void
+blas·syr(int nrow, int ncol, double a, double *x, double *m)
+{
+    int i, j;
+
+    for (i = 0; i < nrow; i++) {
+        for (j = 0; j < ncol; j++) {
+            m[i+ncol*j] += a * x[i] * x[j];
+        }
+    }
+}
+
+#define NITER 500
+#define NCOL  1007
+#define NROW  1007
 
 error
 main()
 {
-    int i;
+    int i, n;
     clock_t t;
-    double  res;
 
     double *x, *y, *m;
+    double res[2], tprof[2];
 
     openblas_set_num_threads(1);
 
@@ -348,33 +377,44 @@ main()
     y = malloc(sizeof(*x)*NCOL);
     m = malloc(sizeof(*x)*NCOL*NROW);
 
-    for (i = 0; i < NCOL; i++) {
-        y[i] = i;
-    }
+    tprof[0] = 0;
+    tprof[1] = 0;
 
-    t = clock();
-    for (i = 0; i < NITER; i++) {
-        cblas_dgemv(CblasRowMajor, CblasNoTrans, NROW, NCOL, 1.5, m, NCOL, x, 1, 2.5, y, 1);
-    }
-    t   = clock() - t;
-    res = cblas_ddot(NROW, y, 1, y, 1);
+    for (n = 0; n < NITER; n++) {
+        for (i = 0; i < NCOL; i++) {
+            x[i] = i*i+1;
+            y[i] = i+1;
+        }
 
-    printf("the result is %f\n", res);
-    printf("time elapsed (blas): %fms\n", 1000.*t/CLOCKS_PER_SEC);
+        for (i = 0; i < NCOL*NROW; i++) {
+            m[i] = .1;
+        }
 
-    for (i = 0; i < NCOL; i++) {
-        y[i] = i;
-    }
+        t = clock();
+        cblas_dger(CblasRowMajor, NROW, NCOL, 1.2, x, 1, y, 1, m, NCOL);
+        t = clock() - t;
+        tprof[1] += 1000.*t/CLOCKS_PER_SEC;
+        res[1]    = cblas_ddot(NROW*NCOL, m, 1, m, 1);
 
-    t = clock();
-    for (i = 0; i < NITER; i++) {
-        blas·gemv(NROW, NCOL, 1.5, m, x, 2.5, y);
-    }
-    t   = clock() - t;
-    res = blas·dot(NCOL, y, y);
+        for (i = 0; i < NCOL; i++) {
+            x[i] = i*i+1;
+            y[i] = i+1;
+        }
 
-    printf("the dot product is %f\n", res);
-    printf("time elapsed (naive): %fms\n", 1000.*t/CLOCKS_PER_SEC);
+        for (i = 0; i < NCOL*NROW; i++) {
+            m[i] = .1;
+        }
+
+        t = clock();
+        blas·ger(NROW, NCOL, 1.2, x, y, m);
+        t = clock() - t;
+        tprof[0] += 1000.*t/CLOCKS_PER_SEC;
+        res[0]    = blas·dot(NROW*NCOL, m, m);
+    }
+    printf("avg time elapsed/iteration (naive): %fms\n", tprof[0]/NITER);
+    printf("--> result (naive): %f\n", res[0]);
+    printf("avg time elapsed/iteration (oblas): %fms\n", tprof[1]/NITER);
+    printf("--> result (oblas): %f\n", res[1]);
 
     return 0;
 }
-- 
cgit v1.2.1