I hope this post does not offend the creator of gsll, as this is not
the intent. There are two ways to read this post: One is that it may
raise a valid point, and the other that I am missing something already
provided, or a mode of using the routines that would make this comment
unnecessary.
So, in effect, I think (as a newbie) I am missing something. But what?
I find the naming convention in gsll inconvenient. For example the
routine solve-tridiagonal is a wrapper for gsl_linalg_solve_tridiag.
To use it to solve a system defined by vectors d e f b, one needs to
use a letm block to this effect:
(letm ((d* (vector-double-float d))
(e* (vector-double-float e))
(f* (vector-double-float f))
(x* (vector-double-float (make-array (length d) :element-type 'double
:initial-element 0d0))
(b* (vector-double-float b)))
(solve-tridiagonal d* e* f* b* x*)
(data x*))
Using the routine solve-tridiagonal is still not convenient, because
one still needs to convert the argument vectors via the
vector-double-float, obscuring the code. If I have multiple calls, I
will have to repeat these, or define my own interface function. And
here is the problem: *the very nice name solve-tridiagonal is taken*.
I could name it solve-tridiag, but sooner or later, I will make a
mistake and call the wrong function with the incorrect arguments.
Ideally, once we have gsll, one would just need to do
(solve-tridiagonal d e f b). To accomplish that, gsll would need to
use the following convention:
- the c-code wrapper would be named the same as the c-code, except
with the underscores replaced with dashes. Thus we would have
(defmfun gsl-linalg-solve-tridiag (diag e f b x)
"gsl_linalg_solve_tridiag"
(((pointer diag) gsl-vector-c) ((pointer e) gsl-vector-c)
((pointer f) gsl-vector-c) ((pointer b) gsl-vector-c)
((pointer x) gsl-vector-c))
:invalidate (x)... etc)
- Users/contributors can define interface routines that would link to
the gsll routine:
(defun solve-tridiag (d e f b)
(letm ((d* (vector-double-float d))
(e* (vector-double-float e))
(f* (vector-double-float f))
(x* (vector-double-float (make-array (length d) :element-type 'double
:initial-element 0d0))
(b* (vector-double-float b)))
(solve-tridiagonal d* e* f* b* x*)
(data x*))
(I have not tried it, but maybe macros should be better to remove the
pass-by-value cost of calling the routine).
So, what am I missing?
Thanks,
Mirko
