Refactor somewhat

implementation/README.md

@@ -64,6 +64,13 @@ Equivalent to 11666383.42 operations per second
 
 As with the Python implementation the benchmarks use 64-bit floating-point arithmetic (doubles). 
 
+## Racket implementation
+
+Unzip `matrices.npz` to a directory called `matrices` in the same
+directory. (I can't get Racket's `unzip` library to work.) Eg, on
+MacOS:
+`unzip matrices.npz -d matrices`
+
 ## Benchmarking Results
 
 The following benchmarking results were obtained using the following devices:
@@ -76,6 +83,8 @@ The times taken are to perform 16,777,216 matrix multiply operations.
 
 ### Python
 
+Using Numpy.
+
 | Device       | Python (s) | Python (ops/s) |
 |:-------------|-----------:|---------------:|
 | Intel i7     |      39.92 |     420,289.64 |
@@ -94,9 +103,8 @@ The times taken are to perform 16,777,216 matrix multiply operations.
 
 Naive implementation, boxed floats.
 
-| Device   |     s |   ops/s |
-|:---------|------:|--------:|
-|          |       |         |
-|          |       |         |
-| Apple M2 | 39.36 | 426,272 |
-
+| Device       |     s |   ops/s |
+|:-------------|------:|--------:|
+| Intel i7     |       |         |
+| Apple M1 Pro |       |         |
+| Apple M2     | 39.36 | 426,272 |

implementation/matmul-racket/matmul.rkt

@@ -1,56 +1,7 @@
 #lang racket
 
-(require "load.rkt")
-
-
-;; ------------------------------------------------------------
-;; Represent a matrix as a struct
-
-(struct matrix (ncols vec) #:transparent)
-
-(define (make-matrix nrows ncols [v 0])
-  (matrix ncols (make-vector (* nrows ncols) v)))
-
-(define (matrix-nrows m)
-  (quotient (vector-length (matrix-vec m)) (matrix-ncols m)))
-
-;; Return element in the ith row and jth column
-;; Index from 0
-(define (matrix-ref m i j)
-  (vector-ref (matrix-vec m) (+ j (* i (matrix-ncols m)))))
-
-(define (matrix-set! m i j v)
-  (vector-set! (matrix-vec m) (+ j (* i (matrix-ncols m))) v))
-
-(define (matrix-equal-within? delta m1 m2)
-  (and
-   (= (matrix-ncols m1) (matrix-ncols m2))
-   (= (matrix-nrows m1) (matrix-nrows m2))
-   (andmap
-    (λ (v1 v2)
-      (<= (abs (- v1 v2)) delta))
-    (vector->list (matrix-vec m1))
-    (vector->list (matrix-vec m2)))))
-
-;; ------------------------------------------------------------
-;; Mathematics
-
-(define (matrix-mul m1 m2)
-  (let ([I (matrix-nrows m1)]
-        [J (matrix-ncols m2)]
-        [K (matrix-ncols m1)])
-    (unless (= K (matrix-nrows m2))
-      (error "Matrices not compatible"))
-    (define vec
-      (for*/vector #:length (* I J)
-                   ([i (in-range I)]
-                    [j (in-range J)])
-        (for/sum ([k (in-range K)])
-          (* (matrix-ref m1 i k)
-             (matrix-ref m2 k j)))))
-
-    (matrix J vec))
-  )
+(require "load.rkt"
+         "matrix-naive.rkt")
 
 ;; ------------------------------------------------------------
 ;; main
@@ -63,13 +14,11 @@
    (λ (x)
      (map
       (λ (m)
-        (let ([ncols (car m)]
-              [vs    (cdr m)])
-          (matrix ncols vs)))
+        (numpy->matrix m))
       x))
    (load-float64-arrays-from-dir "../testdata/matrices")))
 
-(display "done \n")
+(displayln " done.")
 
 (display "Checking test cases...")
 
@@ -83,13 +32,11 @@
 
 (define *repeats* 32768)
 
-(printf "Trying ~a matrix multipications.\n" (* *repeats* (length as)))
+(printf "Trying ~a matrix multiplications.\n" (* *repeats* (length as)))
 
 (display
  (time
-  (for ([i (in-range *repeats*)])
-    (when (zero? (modulo i 1024))
-      (display "."))
+  (for ([_ (in-range *repeats*)])
     (map matrix-mul as bs))))
 
 

implementation/matmul-racket/matrix-naive.rkt

@@ -0,0 +1,64 @@
+#lang racket/base
+
+(provide
+ (struct-out matrix)
+ make-matrix
+ matrix-nrows
+ matrix-ref
+ matrix-set!
+ matrix-equal-within?
+ matrix-mul
+ numpy->matrix
+ )
+
+;; ------------------------------------------------------------
+;; Represent a matrix as a struct
+
+(struct matrix (ncols vec) #:transparent)
+
+(define (numpy->matrix m)
+  (matrix (car m) (cdr m)))
+
+(define (make-matrix nrows ncols [v 0])
+  (matrix ncols (make-vector (* nrows ncols) v)))
+
+(define (matrix-nrows m)
+  (quotient (vector-length (matrix-vec m)) (matrix-ncols m)))
+
+;; Return element in the ith row and jth column
+;; Index from 0
+(define (matrix-ref m i j)
+  (vector-ref (matrix-vec m) (+ j (* i (matrix-ncols m)))))
+
+(define (matrix-set! m i j v)
+  (vector-set! (matrix-vec m) (+ j (* i (matrix-ncols m))) v))
+
+(define (matrix-equal-within? delta m1 m2)
+  (and
+   (= (matrix-ncols m1) (matrix-ncols m2))
+   (= (matrix-nrows m1) (matrix-nrows m2))
+   (andmap
+    (λ (v1 v2)
+      (<= (abs (- v1 v2)) delta))
+    (vector->list (matrix-vec m1))
+    (vector->list (matrix-vec m2)))))
+
+;; ------------------------------------------------------------
+;; Mathematics
+
+(define (matrix-mul m1 m2)
+  (let ([I (matrix-nrows m1)]
+        [J (matrix-ncols m2)]
+        [K (matrix-ncols m1)])
+    (unless (= K (matrix-nrows m2))
+      (error "Matrices not compatible"))
+    (define vec
+      (for*/vector #:length (* I J)
+                   ([i (in-range I)]
+                    [j (in-range J)])
+        (for/sum ([k (in-range K)])
+          (* (matrix-ref m1 i k)
+             (matrix-ref m2 k j)))))
+
+    (matrix J vec))
+  )