[Gsll-devel] small modification to blas2.lisp

[Mirko Vukovic]

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