[lnkForumImage]
TotalShareware - Download Free Software
Confronta i prezzi di migliaia di prodotti.
Asp Forum
 Home | Login | Register | Search 


 

Forums >

comp.lang.lisp

Call BLAS ddot routine from SBCL

Alexandre

7/17/2015 8:10:00 AM

I am trying to call the BLAS ddot routine from SBCL.

Based on the ddot documentation, its source code, some additional doc found on the web and a working example of calling the dgemm routine, I came up with the following script, which produces an error:

(load-shared-object "libblas.so.3")

(declaim (inline ddot))

(define-alien-routine ("ddot_" ddot) void
(n int :copy)
(dx (* double))
(incx int :copy)
(dy (* double))
(incy int :copy))

(defun pointer (array)
(sap-alien (sb-sys:vector-sap (array-storage-vector array)) (* double)))

(defun dot (dx dy)
(unless (= (length dx) (length dy))
(error "Vectors length does not match"))
(let ((n (length dx))
(result 0.0d0))
(sb-sys:with-pinned-objects (dx dy result)
(ddot n (pointer dx) 1 (pointer dy) 1))))

(defvar *a* (make-array 4 :initial-element 1.0d0 :element-type 'double-float))
(defvar *b* (make-array 4 :initial-element 2.0d0 :element-type 'double-float))
(dot *a* *b*)

which gets me the error below:

arithmetic error FLOATING-POINT-INVALID-OPERATION signalled
[Condition of type FLOATING-POINT-INVALID-OPERATION]

Here is a working example of calling the matrix multiplication function:

;; Matrix multiplication in SBCL using BLAS
;; Credits to Miroslav Urbanek, Charles University in Prague

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


Some links I referred to:
- the ddot documentation (http://www.netlib.org/lapack/explore-html/d5/df6/dd...),
- its source code (http://www.netlib.org/lapack/explore-html/d5/df6/ddot_8f_s...),
- some more doc (https://orion.math.iastate.edu/docs/cmlib...),

Any idea of where the problem is?
1 Answer

Alexandre

7/18/2015 7:32:00 AM

0

Le vendredi 17 juillet 2015 10:09:47 UTC+2, Alexandre Landi a écrit :

> Any idea of where the problem is?

Found it. Credits to Miroslav Urbanek from Charles University in Prague for the hint.

-(define-alien-routine ("ddot_" ddot) void
+(define-alien-routine ("ddot_" ddot) double

(defun dot (dx dy)
(unless (= (length dx) (length dy))
(error "Vectors length does not match"))
- (let ((n (length dx))
- (result 0.0d0))
- (sb-sys:with-pinned-objects (dx dy result)
+ (let ((n (length dx)))
+ (sb-sys:with-pinned-objects (dx dy)

The ddot routine is meant to return a double, not a void. And the result variable is useless. Things are so obvious *after* you realize them :-)