Матричное умножение в Common Lisp

Question

Я пишу программу в CL (с SBCL 1.2.15), которая использует линейную алгебру. В ходе выполнения он часто умножает матрицу на вектор.

Профилировщик показал, что большую часть времени (80%) программа проводит именно то, что умножает матрицу на вектор. Он также показывает, что эта функция делает много consing (80 000 000 для ~ 100 вызовов для матриц 100x100), которая также запускает GC. Сделав что-то подобное в F # (которое имеет преимущество статической проверки типов, но в противном случае обычно не превосходит SBCL), программа F # выполняется в 10 раз быстрее.

Я делаю это неправильно?

(defun matrix-mul (matrix vector dest) 
    "Multiply MATRIX by VECTOR putting the result into DEST. 
Optimized for DOUBLE-FLOAT vectors and matrices" 
    (declare (type (array double-float (* *)) matrix) 
      (type (vector double-float *) vector dest) 
      (optimize (speed 3) (debug 0) (safety 0))) 
    (destructuring-bind (rows cols) (array-dimensions matrix) 
    (dotimes (i rows) 
    (setf (aref dest i) 
      (loop for j below cols sum (the double-float 
              (* (aref matrix i j) (aref vector j)))))))) 

PS. Я также пытался использовать DOTIMES вместо внутреннего LOOP - без разницы

PPS. Следующий параметр может использовать BLAS из CL, но (i) я не ожидаю, что он будет работать в Windows, (ii) потребуется маршалинг матриц/векторов взад и вперед.

Answer 1

Обычная проблема заключается в том, что арифметика с плавающей точкой, здесь с двойными поплавками (независимо от окружающего кода, такого как матричное умножение).

Обычно первое, что нужно сделать с SBCL в таких случаях:

Поместите код в файл и скомпилировать его

Компилятор будет выводить много задач оптимизации, учитывая, что один компилирует для скорости. Затем вам нужно будет изучить заметки и посмотреть, что вы можете сделать.

Здесь, например, информация о LOOP отсутствует информация о типе.

На самом деле существует синтаксис LOOP, чтобы объявить тип переменной суммы. Я не знаю, пользуется ли SBCL тем, что:

(loop repeat 10 
     sum 1.0d0 of-type double-float) 

SBCL 1.3.0 на 32-битный ARM для вашего кода:

* (compile-file "/tmp/test.lisp") 

; compiling file "/tmp/test.lisp" (written 13 DEC 2015 11:34:26 AM): 
; compiling (DEFUN MATRIX-MUL ...) 
; file: /tmp/test.lisp 

1)

; in: DEFUN MATRIX-MUL 
;  (SETF (AREF DEST I) 
;    (LOOP FOR J BELOW COLS 
;     SUM (THE DOUBLE-FLOAT (* # #)))) 
; --> LET* FUNCALL SB-C::%FUNCALL (SETF AREF) 
; ==> 
; (SB-KERNEL:HAIRY-DATA-VECTOR-SET ARRAY SB-INT:INDEX SB-C::NEW-VALUE) 
; 
; note: unable to 
; avoid runtime dispatch on array element type 
; due to type uncertainty: 
; The first argument is a (VECTOR DOUBLE-FLOAT), not a SIMPLE-ARRAY. 

2)

;  (AREF MATRIX I J) 
; --> LET* 
; ==> 
; (SB-KERNEL:HAIRY-DATA-VECTOR-REF ARRAY SB-INT:INDEX) 
; 
; note: unable to 
; avoid runtime dispatch on array element type 
; due to type uncertainty: 
; The first argument is a (ARRAY DOUBLE-FLOAT (* *)), not a SIMPLE-ARRAY. 

3)

;  (AREF VECTOR J) 
; ==> 
; (SB-KERNEL:HAIRY-DATA-VECTOR-REF ARRAY SB-INT:INDEX) 
; 
; note: unable to 
; avoid runtime dispatch on array element type 
; due to type uncertainty: 
; The first argument is a (VECTOR DOUBLE-FLOAT), not a SIMPLE-ARRAY. 

4)

;  (LOOP FOR J BELOW COLS 
;   SUM (THE DOUBLE-FLOAT (* (AREF MATRIX I J) (AREF VECTOR J)))) 
; --> BLOCK LET SB-LOOP::WITH-SUM-COUNT LET SB-LOOP::LOOP-BODY TAGBODY SETQ THE 
; ==> 
; (+ #:LOOP-SUM-8 (THE DOUBLE-FLOAT (* (AREF MATRIX I J) (AREF VECTOR J)))) 
; 
; note: unable to 
; optimize 
; due to type uncertainty: 
; The first argument is a NUMBER, not a (COMPLEX SINGLE-FLOAT). 
; 
; note: unable to 
; optimize 
; due to type uncertainty: 
; The first argument is a NUMBER, not a (COMPLEX DOUBLE-FLOAT). 

5)

; --> BLOCK LET SB-LOOP::WITH-SUM-COUNT LET SB-LOOP::LOOP-BODY TAGBODY WHEN IF 
; --> >= OR LET IF OR THE = IF 
; ==> 
; (= SB-C::X SB-C::Y) 
; 
; note: unable to open code because: The operands might not be the same type. 

6)

;  (DOTIMES (I ROWS) 
;  (SETF (AREF DEST I) 
;    (LOOP FOR J BELOW COLS 
;      SUM (THE DOUBLE-FLOAT #)))) 
; --> DO BLOCK LET TAGBODY UNLESS IF >= IF 
; ==> 
; (< SB-C::X SB-C::Y) 
; 
; note: forced to do static-fun Two-arg-< (cost 53) 
;  unable to do inline fixnum comparison (cost 4) because: 
;  The second argument is a INTEGER, not a FIXNUM. 
;  unable to do inline (signed-byte 32) comparison (cost 6) because: 
;  The second argument is a INTEGER, not a (SIGNED-BYTE 32). 
;  etc. 

7)

;  (LOOP FOR J BELOW COLS 
;   SUM (THE DOUBLE-FLOAT (* (AREF MATRIX I J) (AREF VECTOR J)))) 
; --> BLOCK LET SB-LOOP::WITH-SUM-COUNT LET SB-LOOP::LOOP-BODY TAGBODY WHEN IF 
; --> >= OR LET > IF 
; ==> 
; (> SB-C::X SB-C::Y) 
; 
; note: forced to do static-fun Two-arg-> (cost 53) 
;  unable to do inline fixnum comparison (cost 4) because: 
;  The second argument is a REAL, not a FIXNUM. 
;  unable to do inline (signed-byte 32) comparison (cost 6) because: 
;  The second argument is a REAL, not a (SIGNED-BYTE 32). 
;  etc. 

8)

; --> BLOCK LET SB-LOOP::WITH-SUM-COUNT LET SB-LOOP::LOOP-BODY TAGBODY SETQ THE 
; ==> 
; (+ #:LOOP-SUM-8 (THE DOUBLE-FLOAT (* (AREF MATRIX I J) (AREF VECTOR J)))) 
; 
; note: forced to do static-fun Two-arg-+ (cost 53) 
;  unable to do inline float arithmetic (cost 2) because: 
;  The first argument is a NUMBER, not a DOUBLE-FLOAT. 
;  The result is a (VALUES NUMBER &OPTIONAL), not a (VALUES DOUBLE-FLOAT 
;                &REST T). 
; 
; note: doing float to pointer coercion (cost 13), for: 
;  the second argument of static-fun Two-arg-+ 
; 
; compilation unit finished 
; printed 10 notes