Matrix Multiplication in Common Lisp

Question

I am writing a program in CL (with SBCL 1.2.15) that uses linear algebra. At run time, it often multiplies the matrix by a vector.

The profiler showed that most of the time (80%) the program spends this way, multiplying the matrix by the vector. It also shows that this function does a lot of consing (80,000,000 for ~ 100 calls for 100x100 matrices), which also starts the GC. Having done something similar in F # (which has the advantage of static type checking, but otherwise usually does not exceed SBCL), the F # program runs 10 times faster.

Am I doing it wrong?

(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 ij) (aref vector j)))))))) 

PS. I also tried using DOTIMES instead of internal LOOP - no difference

SFC. The next parameter can use BLAS from CL, but (i) I do not expect it to work on Windows, (ii) it will require matrix / vector marshaling back and forth.

Answer 1

A common problem is that the float arithmetic, here with double floats (regardless of the surrounding code, such as matrix multiplication), agrees.

Typically, the first thing to do with SBCL in such cases is:

Put the code in a file and compile it

The compiler then raises many optimization problems, given that one compiles for speed. Then you will need to study the notes and see what you can do.

Here, for example, the LOOP sum does not contain type information.

In fact, there is LOOP syntax for declaring a variable sum type. I do not know if SBCL uses the fact that:

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

SBCL 1.3.0 on 32-bit ARM for your code:

 * (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 IJ) ; --> 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 IJ) (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 IJ) (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 IJ) (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. 

eight)

 ; --> BLOCK LET SB-LOOP::WITH-SUM-COUNT LET SB-LOOP::LOOP-BODY TAGBODY SETQ THE ; ==> ; (+ #:LOOP-SUM-8 (THE DOUBLE-FLOAT (* (AREF MATRIX IJ) (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