Re: calling Fortran code from Lisp

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