Matrix-умножение с использованием BLAS из Common Lisp

Question

Допустим, у меня есть две матрицы (в виде Common Lisp массива) Foo и бар, что:

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

Я хотел бы выполнить умножение матриц с помощью BLAS без с использованием оберток, таких как Matlisp, GSLL, LLA, & co. так что я получаю массив с результатом:

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

Какие шаги следует предпринять для выполнения такой операции?

Мое понимание заключается в том, что я должен вызвать функцию умножения матрицы BLAS из REPL и передать ей свои аргументы foo и bar.

В R, я легко могу это сделать так:

foo %*% bar 

Как я могу сделать это в Common Lisp?

Отказ от ответственности: 1) Я использую SBCL 2) Я не опытный ученый

Answer 1

Вот идеальный ответ, который я искал. Кредиты Мирославу Урбанеку из Карлова университета в Праге.

«Вот основная идея. Я нахожу функцию я хочу использовать от BLAS/LAPACK. В случае умножения матриц, это DGEMM.„D“означает для двойного поплавка,„GE“обозначает общие матрицы (без специального формы, как симметричные, треугольные и т.д.), и «ММ» обозначает матрицу умножения документации здесь:.

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

Тогда я определяю инопланетную процедуру, используя SBCL FFI я прохожу. Lisp array непосредственно с использованием некоторых специальных функций SBCL. Массивы Lisp m ust be , созданный с опцией: double-float типа element-type.

Важным моментом является то, что 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)) 

Answer 2

В R вы используете R обертку. Вы не можете избежать использования «обертки». Поэтому вы должны использовать то, что вам подходит.

Извините, если это не очень полезно, но так оно и есть.

Marco