Add naive Racket matrix-multiplication

implementation/README.md

@@ -70,13 +70,33 @@ The following benchmarking results were obtained using the following devices:
 
 1. 12th Gen Intel Core i7-1260P running Ubuntu 22.04.4, Python 3.10.12 and NumPy 1.25.5.
 2. Apple M1 Pro running macOS 13.6.3, Python 3.12.2 and NumPy 1.24.2.
+3. Apple M2 (powered), macOS 14.4, Python 3.12.2, numpy 1.26.4
 
-The times taken are to perform 16777216 matrix multiply operations.
+The times taken are to perform 16,777,216 matrix multiply operations.
 
-| Device   | Python (s) | C (s) | Python (ops/s) | C (ops/s)   |
-|:---------|-----------:|------:|---------------:|------------:|
-| Intel i7 | 39.92      | 1.44  | 420289.64      | 11666383.42 |
-| Apple M1 | 28.70      | 1.54  | 584562.19      | 10873615.05 |
+### Python
 
-This makes the simple C implementation 27.72 times faster on the Intel processor and 18.64 times faster on the M1 processor, compared to the equivalent NumPy implementation running on the same machine.
+| Device       | Python (s) | Python (ops/s) |
+|:-------------|-----------:|---------------:|
+| Intel i7     |      39.92 |     420,289.64 |
+| Apple M1 Pro |      28.70 |     584,562.19 |
+| Apple M2     |      29.37 |        571,142 |
+
+### C
+
+| Device       | C (s) |     C (ops/s) |
+|:-------------|------:|--------------:|
+| Intel i7     |  1.44 | 11,666,383.42 |
+| Apple M1 Pro |  1.54 | 10,873,615.05 |
+| Apple M2     |  1.02 |    16,525,369 |
+
+### Racket
+
+Naive implementation, boxed floats.
+
+| Device   |     s |   ops/s |
+|:---------|------:|--------:|
+|          |       |         |
+|          |       |         |
+| Apple M2 | 39.36 | 426,272 |
 

implementation/matmul-racket/load.rkt

@@ -0,0 +1,69 @@
+#lang racket
+
+(require binaryio)
+
+;; Returns (as bs cs) where xs is a list of "arrays"
+;; Each array is a pair (n . v), where v is a vector
+;; and n the length of the rows.
+(provide load-float64-arrays-from-dir)
+
+(struct npy-header
+  (magic-string
+   major-version-number
+   minor-version-number
+   length
+   format)
+  #:transparent)
+
+(define (load-float64-arrays-from-dir dir)
+  (define filenames (directory-list dir))
+  ;; each filename is X_999.npy, where X is a, b, or c and 999 is a number->string
+
+  (define n-matrices
+    (+ 1
+       (apply max
+              (map
+               (λ (entry)
+                 (string->number
+                  (cadr (regexp-match #rx"[abc]_([0-9]+).npy" (path->string entry)))))
+               filenames))))
+
+  ;; returns (values as bs cs)
+  (define-values (as bs cs)
+    (for/lists (as bs cs)
+               ([i (in-range n-matrices)])
+      (when (zero? (modulo i (quotient n-matrices 8)))
+        (display "."))
+      (apply values
+             (for/list ([x '("a" "b" "c")])
+               (call-with-input-file
+                 (build-path dir (string-append x "_" (number->string i) ".npy"))
+                 load-float64-array-from-npy)))))
+
+  (list as bs cs))
+
+
+(define (load-float64-array-from-npy in)
+  ;; Read header
+  (define the-header 
+    (let* ([mgc (read-bytes 6 in)]
+           [maj (read-integer 1 #f in)]
+           [min (read-integer 1 #f in)]
+           [len (read-integer 2 #f in #f)]
+           [fmt (bytes->string/latin-1 (read-bytes len in))])
+      (npy-header mgc maj min len fmt)))
+
+  ;; Parse shape from a Python dictionary (!)
+  (define the-shape
+    (let ([shape (regexp-match #rx"'shape': \\(([0-9]+), ([0-9]+)\\)"
+                               (npy-header-format the-header))])
+      (list (string->number (cadr shape)) 
+            (string->number (caddr shape)))))
+
+  ;; Read array and return it and the length of a row
+  (cons (cadr the-shape)
+        (for/vector ([_ (in-range (* (car the-shape) (cadr the-shape)))])
+          (read-float 8 in #f))))
+
+
+

implementation/matmul-racket/matmul.rkt

@@ -0,0 +1,95 @@
+#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))
+  )
+
+;; ------------------------------------------------------------
+;; main
+
+
+(display "Loading test data ")
+
+(match-define (list as bs cs)
+  (map
+   (λ (x)
+     (map
+      (λ (m)
+        (let ([ncols (car m)]
+              [vs    (cdr m)])
+          (matrix ncols vs)))
+      x))
+   (load-float64-arrays-from-dir "../testdata/matrices")))
+
+(display "done \n")
+
+(display "Checking test cases...")
+
+(unless
+    (andmap
+     (λ (m1 m2) (matrix-equal-within? 0.01 m1 m2))
+     cs (map matrix-mul as bs))
+  (error "something went wrong!"))
+
+(displayln " done.")
+
+(define *repeats* 32768)
+
+(printf "Trying ~a matrix multipications.\n" (* *repeats* (length as)))
+
+(display
+ (time
+  (for ([i (in-range *repeats*)])
+    (when (zero? (modulo i 1024))
+      (display "."))
+    (map matrix-mul as bs))))
+
+