[Gsll-devel] small modification to blas2.lisp

[Mirko Vukovic]

On Sat, May 17, 2008 at 8:56 PM, Mirko Vukovic <mirko.vukovic at gmail.com> wrote:
> On Sat, May 17, 2008 at 6:11 PM, Liam Healy <lhealy at common-lisp.net> wrote:
>> Mirko,
>>
>> Your error is that your function isn't defined in the GSL library;
>> that is what sb-kernel::undefined-alien-function-error means.  A quick
>> glance at "fft_complex_radix2_forward" makes me suspicious, as I
>> believe all GSL functions start with "gsl_" and this doesn't have
>> that.  Since you're on linux, use the shell script "list" which is
>> at top level in the GSLL files, thus:
>>  list | grep -i fft_complex_radix2_forward
>>  gsl_fft_complex_radix2_forward
>>

You were absolutely correct.  The name is gsl_fft...

I am now getting a memory fault:

Unhandled memory fault at #x0.
   [Condition of type SB-SYS:MEMORY-FAULT-ERROR]

Restarts:
  0: [ABORT-REQUEST] Abort handling SLIME request.
  1: [ABORT] Exit debugger, returning to top level.

Backtrace:
  0: (SB-SYS:MEMORY-FAULT-ERROR)
  1: (SB-SYS:MEMORY-FAULT-ERROR)
  2: ("foreign function: #x419F72")
  3: ("foreign function: #x41A040")
  4: ("foreign function: #x3919E4AEA4")
  5: ((LAMBDA (SB-PCL::.PV. SB-PCL::.NEXT-METHOD-CALL.
SB-PCL::.ARG0.)) #<error printing object>)

I have modified the defmfun to include a return declaration as :int
and also declared the function as a :function.  But that did not
change the above error.  Here is the new code.

(in-package :gsll)

(defmfun fft-c2f (x stride n)
  "gsl_fft_complex_radix2_forward"
;; for gsl doc and example see
;; http://www.gnu.org/software/gsl/manual/html_node/Radix_002d2-FFT-routines-for-complex-data.html
  (((pointer x) gsl-vector-c) (stride :int) (n :int))
  :type :function
  :c-return :int
  :documentation
  "Forward FFT for a complex double radix-2 vector")

;; test run
(let ((arg (make-array 4 :element-type 'complex :initial-element #c(0d0 0d0))))
  (setf (aref arg 2) #c(1d0 0d0))
  (let* ((dim (length arg))
	(double (make-array (* 2 dim)
			     :element-type 'double-float
			     :initial-element 0d0)))
    ;; repackaging complex as double -- is there a built-in?
    (loop
       for re-im across arg
       for i from 0 to (* 2 (1- dim)) by 2
       do (progn
	    (setf (aref double i) (realpart re-im))
	    (setf (aref double (1+ i)) (imagpart re-im))))

    (letm ((double* (vector-double-float double)))
      (fft-c2f double* 1 dim))))

I am essentially trying to replicate example3.lisp from nlisp's
distribution.  I am including it here:
;;; Low level wrapper for gsl fft routine
(cffi:defcfun "gsl_fft_complex_radix2_forward"
	      :int (a :pointer) (stride :int) (n :int))

;;; A high level NLISP style function
(defun fft (a)
  (let ((rn (.* a 1d0)))
    (with-accessors ((r val)) rn
       (sb-sys:without-gcing)
       (let ((r-addr (sb-sys:vector-sap r)))
	 (declare (type (simple-array (complex double-float)) r))
	 (gsl-fft-complex-radix2-forward r-addr 1 (array-dimension r 0)))) rn))

;;; Example program
(defun example3 (n)
  (let ((a (make-complex-double-array 128)))
    (setf (.aref a n) (complex 1d0 0d0))
    (let ((z (fft a)))
      (nplot (list (.realpart z) (.imagpart z)) nil
	     :legends (list "Real Part" "Imaginary Part")))))


Should I wait to for you to finish your transition to ffa-based arrays
before pursuing this further?

I could check the correctness of the call once I have something working right.

Mirko