使用 Common Lisp 的 BLAS 进行矩阵乘法

Question

假设我有两个矩阵（以 Common Lisp 数组的形式）foo 和 bar 这样：

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

我想使用 BLAS 不使用 使用 Matlisp、GSLL、LLA 等包装器执行矩阵乘法。这样我得到一个结果数组：

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

我应该采取哪些步骤来执行这样的操作？

我的理解是，我应该从 REPL 中调用 BLAS 矩阵乘法函数，并将我的参数 foo 和 bar 传递给它。

在 R 中，我可以很容易地这样做：

foo %*% bar

如何在 Common Lisp 中实现？

免责声明： 1) 我使用 SBCL 2) 我不是经验丰富的计算机科学家

Answer 1

在 R 中，您使用的是 R 包装器。您无法避免使用 "wrapper"。所以你应该使用最适合你的。

抱歉，如果这没有太大帮助，但事情就是这样。

马可

Answer 2

这是我一直在寻找的完美答案。感谢布拉格查理大学的 Miroslav Urbanek。

"这是基本思路。我从中找到了一个我想使用的函数 BLAS/LAPACK。在矩阵乘法的情况下，它是 DGEMM。 "D" 站对于双浮点数，"GE" 代表一般矩阵（没有特殊的形状如对称、三角形等），"MM"代表矩阵乘法。文档在这里：

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

然后我使用 SBCL FFI 定义了一个外星例程。我通过 Lisp 数组直接使用一些特殊的 SBCL 函数。 Lisp 数组必须是使用选项创建 :element-type 'double-float.

重要的一点是SBCL以行优先存储数组元素顺序，与 C 类似。Fortran 使用列优先顺序。这个有效地对应于转置矩阵。矩阵的顺序因此在从调用 DGEMM 时必须更改它们的尺寸口齿不清。"

;; Matrix multiplication in SBCL using BLAS
;; Miroslav Urbanek <mu@miroslavurbanek.com>

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

使用 Common Lisp 的 BLAS 进行矩阵乘法

Matrix-multiplication using BLAS from Common Lisp

common-lisp

matrix-multiplication

blas

sbcl