Matrix-multiplication using BLAS from Common Lisp

问题:

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

回答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 <[e-mail]>

(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))

回答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