calling Fortran code from Lisp

Hi,

Is it possible to call a Fortran subroutine from Lisp?  I
have seen
the C language FFI in the docs, but I wonder if it would be
possible
to call Fortran somehow.  There is a subroutine for banded
sparse
matrices in LAPACK I would like to use.

Thanks,

Tamas

Possible, and easy.
I actually have an autogenerated wrapper for all of LAPACK
using CFFI:
http://middleangle.com/rif/derifatives
/article/21/cl-blapack-alpha-release

I have not personally run this on cmucl (only sbcl), but it
should be
relatively straightforward.  

There is also MATLISP, which is a higher-level and therefore
more
pleasant interface, although I don't know if it includes
the
functionality you wnat.  I have added new LAPACK functions
to MATLISP in
the past, but it's not an enjoyable activity.

rif

Tamas K Papp writes:
> Is it possible to call a Fortran subroutine from Lisp? 
I have seen
> the C language FFI in the docs, but I wonder if it
would be possible
> to call Fortran somehow.  There is a subroutine for
banded sparse
> matrices in LAPACK I would like to use.

If it's possible to call Fortran from C, yes.
Just declare the function with the FFI like you'd do it in
C.

If the call protocol of your Fortran compiler is different
than that
of your C compiler, then something might still be possible,
using a
library like ffcall/avcall (used by clisp, I don't know if
cmucl uses
it too, but you could write a FFI to that library to then
write the
FFI to the Fortran call protocol).

-- 
__Pascal Bourguignon__
http://www.informatimago
.com
http://pjb.ogamita.org

>>>>> "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.

Ray

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

Tamas K Papp wrote:
> 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.
>
>   
Matlisp only handles 2-D matrices.  You can mix and match
Lisp arrays 
and Matlisp matrices, with some care.
> 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?
>
>   
vector-sap only works on specialized arrays.  For other
array types, you 
need to get to the actual storage area that hold the data. 

(%array-data? or with-array-data?)  Once you have that
object, 
vector-sap will do what you want.
> If you have any advice on what to do, or how to improve
this code in
> other ways, please tell me.
>
>
>   
[snip code]

What happens when you run this code?  I can't easily run
this myself 
right now.
One note, however.  LAPACK matrices are in column-major
order but Lisp 
is row-major order.  You might need to transpose your
matrices before 
calling LAPACK routines.  This is why Matlisp stores its
matrices as a 
1D vector internally.  Perhaps the C interface doesn't need
this.

Ray