Matrix-multiplication using BLAS from Common Lisp

Question

Let's say I have two matrices (in the form of a Common Lisp array) foo and bar such that:

(defvar foo #2A((2 1 6) (7 3 4)))
(defvar bar #2A((3 1) (6 5) (2 3)))

I would like to perform a matrix multiplication using BLAS without using wrappers such as Matlisp, GSLL, LLA, & co. so that I get an array with the result:

#2A((24 25) (47 34))

Which steps should I take to perform such operation?

My understanding is that I should call the BLAS matrix multiplication function from the REPL and pass it my arguments foo and bar.

In R, I can easily do it like this:

foo %*% bar

How can I do it in Common Lisp?

Disclaimer: 1) I use SBCL 2) I am not a seasoned computer scientist

Podesta

Answer 1

Here's the perfect answer I was looking for. Credits to Miroslav Urbanek from Charles University in Prague.

"Here's the basic idea. I find a function I want to use from BLAS/LAPACK. In case of matrix multiplication, it's DGEMM. "D" stands for double float, "GE" stands for general matrices (without a special shape like symmetric, triangular, etc.), and "MM" stands for matrix multiplication. The documentation is here:

http://www.netlib.org/lapack/explore-html/d7/d2b/dgemm_8f.html

Then I define an alien routine using SBCL FFI. I pass Lisp array directly using some special SBCL functions. The Lisp arrays must be created with an option :element-type 'double-float.

An important point is that SBCL stores array elements in row-major order, similarly to C. Fortran uses column-major order. This effectively corresponds to transposed matrices. The order of matrices and their dimensions must be therefore changed when calling DGEMM from Lisp."

;; Matrix multiplication in SBCL using BLAS
;; Miroslav Urbanek <[email protected]>

(load-shared-object "libblas.so.3")

(declaim (inline dgemm))

(define-alien-routine ("dgemm_" dgemm) void
  (transa c-string)
  (transb c-string)
  (m int :copy)
  (n int :copy)
  (k int :copy)
  (alpha double :copy)
  (a (* double))
  (lda int :copy)
  (b (* double))
  (ldb int :copy)
  (beta double :copy)
  (c (* double))
  (ldc int :copy))

(defun pointer (array)
  (sap-alien (sb-sys:vector-sap (array-storage-vector array)) (* double)))

(defun mm (a b)
  (unless (= (array-dimension a 1) (array-dimension b 0))
    (error "Matrix dimensions do not match."))
  (let* ((m (array-dimension a 0))
     (n (array-dimension b 1))
     (k (array-dimension a 1))
     (c (make-array (list m n) :element-type 'double-float)))
    (sb-sys:with-pinned-objects (a b c)
      (dgemm "n" "n" n m k 1d0 (pointer b) n (pointer a) k 0d0 (pointer c) n))
    c))

(defparameter a (make-array '(2 3) :element-type 'double-float :initial-contents '((2d0 1d0 6d0) (7d0 3d0 4d0))))
(defparameter b (make-array '(3 2) :element-type 'double-float :initial-contents '((3d0 1d0) (6d0 5d0) (2d0 3d0))))

(format t "a = ~A~%b = ~A~%" a b)

(defparameter c (mm a b))

Podesta

Answer 2

In R you are using the R wrapper. You cannot avoid using a "wrapper". So you should use that best suits you.

Sorry if this isn't much helpful, but that's how things are.

Marco

Reflexive

My current understanding is that you need some foreign function interface binding to ensure the link between the BLAS function and its call from Common Lisp. If my understanding is accurate, in that sense, yes, you do need a wrapper. Rephrasing my question with some more precision would result in: how can I directly use a ffi without using any other higher level Lisp library such as Matlisp, etc? I am trying to understand how I can perform a matrix multiplication on two CL arrays foo and bar or compute a covariance matrix on some array baz by using already-existing BLAS functions. – Podesta