You are viewing lhealy

Liam Healy - Post a comment

Jan. 15th, 2009

06:15 pm - Passing a complex number to a foreign function

GSL has a number of functions that take complex scalars as arguments. It defines them with either

typedef struct
    double dat[2];

typedef struct
    float dat[2];

These complex numbers are always passed by value, not as a pointer. This presents a problem to Lisp users as many foreign function interfaces and CFFI do not allow passing structs directly.

I thought maybe there was a way around this. Could it be that the real and imaginary parts are simply passed sequentially as if they were two successive arguments to the function? In fact, this does work for numbers of type (complex double-float), but when I tried that on (complex single-float), it would return without error but with the result computed as if the imaginary part was 0. I thought maybe there was a padding issue, so I packed the two single-floats into a single 64-bit word that I then interpreted as a double-float. When passing this single argument to a function that expects a (complex single-float), the correct answer is computed.

For example, the BLAS function caxpy works on single-floats:

(FOREIGN-FUNCALL "gsl_blas_caxpy" 

where #'pack-complex-as-double does the packing on the complex scalar alpha described, and the function zaxpy works on double-floats
(FOREIGN-FUNCALL "gsl_blas_zaxpy" 

This has now been incorporated into GSLL. I assume it is not portable; I am on Debian amd64 but it works on x86 (32 bit) as well. It seems to be portable across CL implementations, at least it works in SBCL and CCL. So I have built in a test that tries a simple function; if it gets the wrong answer than all the functions that want a complex scalar will signal an error.

Tags: ,
Current Mood:

Leave a comment:

No HTML allowed in subject


Notice! This user has turned on the option that logs your IP address when posting. 

(will be screened)