On Fri, Apr 06, 2007 at 05:11:27PM -0400, Raymond Toy wrote:
> >>>>> "Tamas" == Tamas K Papp <Tamas> writes:
>
> Tamas> Hi,
> Tamas> Is it possible to call a Fortran subroutine from Lisp? I have seen
> Tamas> the C language FFI in the docs, but I wonder if it would be
> possible
> Tamas> to call Fortran somehow. There is a subroutine for banded sparse
> Tamas> matrices in LAPACK I would like to use.
>
> As mentioned by others, there are already libraries that can do this.
> If you want to cook your own, you do need to exercise some care. All
> Fortran compilers I know of have, essentially, a C-compatible
> interface. The only thing you need to know about is how Fortran
> handles strings.
>
> If the existing libraries don't do what you want, please ask again
> here, and we can help with setting up a call to a Fortran routine.
Hi Raymond,
First, a short explanation why I am not using the existing libraries:
I need the ability to reshape a more than 2 dimensional array into a
matrix, do some matrix manipulation (eg matrix product, solving linear
system), reshape the result into an array, etc. I don't know how to
do this in either Matlisp or rfi's libraries.
I tried calling dgemm to do a simple matrix multiplication first.
Fortunately Atlas has a C interface, so I can use that. But I got
stuck because I don't know how to extract the pointer of a simple
array, the sys:vector-sap mentioned in the user's guide doesn't work
here. What should I use instead?
If you have any advice on what to do, or how to improve this code in
other ways, please tell me.
(ext:load-foreign "/usr/lib/sse2/libcblas.so") ; you might have to change this
(alien:def-alien-routine "cblas_dgemm" c-call:void
(order c-call:int :in) ; row- or column-major
(transa c-call:int :in) ; whether to transpose A (op)
(transb c-call:int :in) ; whether to transpose B (op)
(m c-call:int :in) ; number of rows in op(A)
(n c-call:int :in) ; number of columns in op(B)
(k c-call:int :in) ; common dimension of op(A) and op(B)
(alpha c-call:double :in) ; alpha
(a (* double-float) :in) ; A
(lda c-call:int :in) ; lda
(b (* double-float) :in) ; B
(ldb c-call:int :in) ; ldb
(beta c-call:double :in) ; beta
(c (* double-float) :in) ; C
(ldc c-call:int :in)) ; ldc
(defun matrix-multiply (a b)
(assert (typep a '(simple-array double-float (* *))))
(assert (typep b '(simple-array double-float (* *))))
(let* ((adim (array-dimensions a))
(bdim (array-dimensions b))
(m (first adim)) ; number of rows in A
(k (second adim)) ; number of columns in A
(n (second bdim)) ; number of columns in B
(c (make-array (list m n) :element-type 'double-float)))
(assert (= k (first bdim))) ; check that matrices are compatible
(sys:without-gcing
(let ((a-addr (sys:vector-sap a))
(b-addr (sys:vector-sap b))
(c-addr (sys:vector-sap c)))
(cblas-dgemm 102 111 111 ; column-major, don't transpose A or B
m n k ; dimensions
1d0 ; alpha
a-addr ; a
m ; number of rows in A
b-addr ; b
n ; number of cols in B
0d0 ; beta
c-addr ; c
m))) ; number of rows in C
c))
;; code to test
(defparameter *a* (make-array '(2 2) :element-type 'double-float
:initial-contents '((1d0 2d0)
(3d0 4d0))))
(defparameter *b* (make-array '(2 2) :element-type 'double-float
:initial-element 1d0))
(matrix-multiply *a* *b*)
Thanks,
Tamas
