On Wed, May 14, 2008 at 5:02 PM, Liam Healy lhealy@common-lisp.net wrote:

On Wed, May 14, 2008 at 2:23 PM, Mirko Vukovic mirko.vukovic@gmail.com wrote:

...

Hi,

I was playing with gsll and blas. I modified the cblas-transpose definition in blas2.lisp (see top of file) to capture a few more values:

Here is the new version:

(cffi:defcenum cblas-transpose "CBLAS_TRANSPOSE from /usr/include/gsl/gsl_cblas.h." (:notrans 111) (:trans 112) (:conjtrans 113))

I fail to see the difference with what is in there now: (cffi:defcenum cblas-transpose "CBLAS_TRANSPOSE from /usr/include/gsl/gsl_cblas.h." (:notrans 111) :trans :conjtrans) See http://common-lisp.net/project/cffi/manual/html_node/defcenum.html#defcenum "If value is omitted its value will either be 0, if it's the first entry, or it it will continue the progression from the last specified value." so :trans gets 112, :conjtrans 113 by default.

I was wondering why the other values were omitted.

...

I am posting it here because I was not able to download with svn so I don't have versioning for gsll.

Also, here is a code snippet that I used to test for tridiagonal matrix solving that also uses gsll-s blas2 routine gemv to double check that things work.

(in-package :gsll) ;; solve a tridiagonal system and then verify result by matrix multiplication ;; uses solve-tridiagonal and gemv (let* ((dim 6) (dim- (- dim 1))) (let ((d (make-array dim :element-type 'double-float :initial-element -2.d0)) (e (make-array dim- :element-type 'double-float :initial-element 0.2d0)) (f (make-array dim- :element-type 'double-float :initial-element 1.d0)) (x (make-array dim :element-type 'double-float :initial-element 0.d0)) (b (make-array dim :element-type 'double-float :initial-element 0.5d0)) (mat (make-array (list dim dim) :element-type 'double-float :initial-element 0.d0)))

;; load matrix
(loop for i from 0 below dim-
   do (progn
        (setf (aref mat i i) (aref d i))
        (setf (aref mat (+ i 1) i) (aref f i))
        (setf (aref mat i (+ i 1)) (aref e i))))
(setf (aref mat dim- dim-) (aref d dim-))

;; gsl operations
(letm ((d* (vector-double-float d))
       (e* (vector-double-float e))
       (f* (vector-double-float f))
       (x* (vector-double-float x))
       (b* (vector-double-float b))
       (v2* (vector-double-float
                   (make-array dim :element-type 'double-float
                                   :initial-element 0.d0)))
       (mat* (matrix-double-float mat)))
  (solve-tridiagonal d* e* f* b* x*)
  ;; gemv returns the result in the last argument
  (gemv :notrans 1.d0 mat* x* 0.d0 v2*)
  (data v2*))))

I am sure this code betrays my lack of lisp experience. So feel free to clean it up if you want it included in the verification part.

Mirko

I might make a few simplifications; #'1- to do decrementing and I don't think there's a need to make arrays in CL and then remake them on the foreign side. It is possibly a candidate for example/ regression tests.

Liam

I agree for regression that there is no need to duplicate the arrays. But,

1) I will be using the arrays outside of (letm ...) to plot them etc 2) I am still feeling my way around. Even the above point may be un-necessary, with, for example the (data ... function)

Mirko

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

In the transition FFA-based arrays, I am incoporating have support for all GSL array types, including complex. That, with some valid examples of how to do FFT might allow for its addition to GSLL.

Liam

On Fri, May 16, 2008 at 8:52 AM, Mirko Vukovic mirko.vukovic@gmail.com wrote:

Hello.

I am trying to use gsl's fft routines and I am having trouble. Here is the sample code where I have defined the interface followed by a little driver.

...

From my reading of the documentation and the example (see

http://www.gnu.org/software/gsl/manual/html_node/Radix_002d2-FFT-routines-fo...) , the fft routine needs the complex vector repackaged as a double vector, a stride index, and the size. That is what I tried to do below.

I get the error sb-kernel::undefined-alien-function-error. I am using sbcl 1.0.14 on linux.

The code is at the end of the message.

Thanks,

Mirko

(in-package :gsll)

(defmfun fft-c2f (x stride n) "fft_complex_radix2_forward" ;; for gsl doc and example see ;; http://www.gnu.org/software/gsl/manual/html_node/Radix_002d2-FFT-routines-fo... (((pointer x) gsl-vector-c) (stride :int) (n :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)) (letm ((dim (length arg)) (double* (vector-double-float ;; repackaging complex as double -- is there a built-in? (let ((double (make-array (* 2 dim) :element-type 'double-float :initial-element 0d0))) (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)))) double)))) (fft-c2f double* 1 dim) (data double*))) _______________________________________________ Gsll-devel mailing list Gsll-devel@common-lisp.net http://common-lisp.net/cgi-bin/mailman/listinfo/gsll-devel

On Sat, May 17, 2008 at 8:56 PM, Mirko Vukovic mirko.vukovic@gmail.com wrote:

On Sat, May 17, 2008 at 6:11 PM, Liam Healy lhealy@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-fo... (((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

On Tue, May 20, 2008 at 6:04 PM, Liam Healy lhealy@common-lisp.net wrote:

Mirko,

Thanks for posting your port of the FFT code. As I say on the documentation page, the reason I held off on the FFT port was mainly because I had no working example to guide me, and I've not ever done FFTs myself. I did see the NLISP work on FFTs later and thought that would be a good place to start when I got around to thinking about FFTs again. So your work is very valuable.

As far as waiting for the ffa transition, it is up to you. Though I've been working steadily on it, it is not close to completion. And, I will be travelling the next couple weeks, so not much is going to happen in the near future anyway.

Keep us posted on further improvements.

Liam

Liam,

I am sorry if I caused confusion, but my code crashes with an unhandled memory fault. If you/anyone could help me find the fault, I could work on the "benchmark" - I am familiar with dft's.

In my long email, I have included the example from nlisp that succesfully implements a call to the same gsl's fft routine. So, it is a matter of finding where my invocation of gsll's facilities is causing this problem.

Mirko

Mirko