M*M

mmtest.lisp

(defpackage #:matrix-multiplication-test
  (:nicknames #:mm-test)
  (:use #:cl
	#:clml.hjs.matrix))

(in-package #:matrix-multiplication-test)

(setf *read-default-float-format* 'double-float)
(defparameter *n* 150)

(defun make-matrix (rows cols)
  (make-dmat rows cols))

(defun randomize (m)
  (let ((rows (nrow m))
	(cols (ncol m)))
    (dotimes (row rows)
      (dotimes (col cols)
	(setf (aref m row col) (random 1.0)))))
  m)

(defmacro time&cycle (form)
  `(format t
	   "~A~%"
	   (subseq (with-output-to-string (*trace-output*)
		     (time ,form))
		   19 24)))

(defun run (gemm &optional (n 1000))
  (let ((ma (randomize (make-matrix *n* *n*)))
	(mb (randomize (make-matrix *n* *n*))))
    (time&cycle
     (loop :repeat n
	   :do (funcall gemm ma mb)
	   :finally (return nil)))))

(defmacro dotimes-unroll ((i n unroll) &body body)
  (let ((n_      (gensym "n")))
    `(let ((,n_ ,n))
       (do ((,i 0))
           ((< ,n_ (the fixnum (+ ,unroll ,i)))
            (do ((,i ,i (the fixnum (1+ ,i))))
                ((< ,n_ (the fixnum (1+ ,i))))
              ,@body
              ))
         ,@(loop :repeat unroll
		 :append (append body `((setq ,i (the fixnum (1+ ,i))))))))))

(defun on-register-gemm (ma mb)
  (declare (optimize (speed 3) (debug 0) (safety 0)))
  (declare (type (simple-array double-float (* *)) ma mb))
  (let ((rows (array-dimension ma 0))
        (cols (array-dimension mb 1)))
    (declare (type fixnum rows cols))
    (let ((result (make-matrix rows cols)))
      (declare (type (simple-array double-float (* *)) result))

      (dotimes (row rows)
        (dotimes (col cols)
          (let ((res (aref result row col)))
           (dotimes-unroll (k cols 16)
             (setf res
                   (the double-float (+ res (* (aref ma row k) (aref mb k col))))))
           (setf (aref result row col) res))))
      result)))

(defun caching-gemm (ma mb)
  (declare (optimize (speed 3) (debug 0) (safety 0)))
  (declare (type (simple-array double-float (* *)) ma mb))
  (let ((rows (array-dimension ma 0))
        (cols (array-dimension mb 1)))
    (declare (type fixnum rows cols))
    (let ((result (make-matrix rows cols)))
      (declare (type (simple-array double-float (* *)) result))

      (dotimes (row rows)
        (dotimes (col cols)
          (let ((cell (aref result row col)))
           (dotimes (k cols)
             (setf (aref result row col)
              (+ cell (* (aref ma row k) (aref mb k col)))))
           )))
      result)))

(defun simple-gemm (ma mb)
  (declare (optimize (speed 3) (debug 0) (safety 0)))
  (declare (type (simple-array double-float (* *)) ma mb))
  (let ((rows (array-dimension ma 0))
        (cols (array-dimension mb 1)))
    (declare (type fixnum rows cols))
    (let ((result (make-matrix rows cols)))
      (declare (type (simple-array double-float (* *)) result))
      (dotimes (row rows)
        (dotimes (col cols)
          (dotimes (k cols)
            (incf (aref result row col)
                  (* (aref ma row k) (aref mb k col))))))
      result)))

(defun test ()
  (loop :for *n* :from 50 :to 400 :by 50
	:do (progn
	      (format t "~%*N*=~A~%" *n*)
	      (format t "SIMPLE-GEMM~18T")
	      (run #'simple-gemm 100)
	      (format t "M*M(CLML)~18T")
	      (run #'m*m 100)
	      (format t "ON-REGISTER-GEMM~18T")
	      (run #'on-register-gemm 100)
	      (format t "CACHING-GEMM~18T")
	      (run #'caching-gemm 100))))

REPL

MM-TEST> (test)

*N*=50
SIMPLE-GEMM       0.062
M*M(CLML)         0.069
ON-REGISTER-GEMM  0.031
CACHING-GEMM      0.069

*N*=100
SIMPLE-GEMM       0.497
M*M(CLML)         0.545
ON-REGISTER-GEMM  0.237
CACHING-GEMM      0.387

*N*=150
SIMPLE-GEMM       1.660
M*M(CLML)         1.796
ON-REGISTER-GEMM  0.811
CACHING-GEMM      1.283

*N*=200
SIMPLE-GEMM       4.116
M*M(CLML)         4.246
ON-REGISTER-GEMM  2.111
CACHING-GEMM      3.299

*N*=250
SIMPLE-GEMM       8.291
M*M(CLML)         8.139
ON-REGISTER-GEMM  4.552
CACHING-GEMM      7.111

*N*=300
SIMPLE-GEMM       14.21
M*M(CLML)         13.99
ON-REGISTER-GEMM  7.655
CACHING-GEMM      13.12

*N*=350
SIMPLE-GEMM       24.80
M*M(CLML)         23.03
ON-REGISTER-GEMM  13.66
CACHING-GEMM      21.99

*N*=400
SIMPLE-GEMM       69.86
M*M(CLML)         39.30
ON-REGISTER-GEMM  23.87
CACHING-GEMM      32.65