Re: [Gsl-shell-info] [gsl-shell] mail list not working (#3)

[Francesco Abbate]

2012/10/21 sonoro1234 <address@hidden>
>
> I did post to address@hidden after subscribing and checked at 
> http://lists.nongnu.org/archive/html/gsl-shell-info/ but there was not my 
> email. It was:
>
> Should I write ffi interface for using gsl_poly_complex_solve or is there
> something better to do?
>
> Is there any documentation about templates?

Hi Victor,

sorry about the post, there was a problem with the list. Now it should
be find if you send other messages, they will be posted directly
without moderator approval required.

> In the meanwhile I did
>
> '''
> local ffi = require 'ffi'
> local gsl = require 'gsl'
> ffi.cdef[[
> typedef long gsl_poly_complex_workspace;
> typedef double* gsl_complex_packed_ptr;
> gsl_poly_complex_workspace * gsl_poly_complex_workspace_alloc (size_t n);
> void gsl_poly_complex_workspace_free (gsl_poly_complex_workspace * w);
> int gsl_poly_complex_solve (const double * a, size_t n, 
> gsl_poly_complex_workspace * w, gsl_complex_packed_ptr z);
> ]]
>
> function solvePoly(t)
> local n = #t
> local coefs = ffi.new("double[?]", n)
> local roots = ffi.new("double[?]", 2(n-1))
> for i=1,n do coefs[i-1]=t[i] end
> local ws = gsl.gsl_poly_complex_workspace_alloc (n)
> if gsl.GSL_SUCCESS~=gsl.gsl_poly_complex_solve (coefs,n,ws, roots) then
> print("no convergence")
> end
> gsl.gsl_poly_complex_workspace_free(ws)
> local res = {}
> for i=1,n-1 do
> res[i]=complex.new(roots[2(i-1)],roots[2*(i-1)+1])
> end
> return res
> end
> '''

Seems to be a very good starting point, there are just a couple of
trivial errors when you write 2(i-1) instead of 2*(i-1).

What you are doing is correct, you can use the FFI module to access to
gsl functions directly. The only reason the gsl_poly_* functions does
not exists in GSL Shell is lack of time and your contribution can be
welcome.

To allocate the arrays I see that you are using

> ffi.new('double[?]', n)

This is fine but the preferred method is to call "m = matrix.alloc(n,
1)". In this way you have the array in the "m.data" field. You can use
matrix.calloc to allocate a complex matrix and still the .data field
will be a double* as required in your implementation.

> but I am afraid that if it was not in gsl-shell it is because there is a 
> better alternative (roots.lua??) and is inacurate in multiple roots.

The module "root" can be interesting but it does work only for real
functions and is not restricted to polynomials. Unfortunately I never
made the effort of add more functions to make it useful in a wider set
of use cases and documenting it.

> I still dont undestand templates.

Don't worry, you don't need to understand templates :-)
You can think of them as a ugly hack to produce on the fly tight code
so that the JIT compiler is more effective. It was only useful for the
ODE integration module. For all the other modules its use is optional.
In the case of the SF module it was indeed useful to generate easily
all the needed stub code.

> I am also working in a library for digital filter design and I dont know if 
> there is any interest in
> incorporating the library

I'm definitely interested in the implementation of a wrapper for the
gsl_poly functions.

About the library for digital filter this is definitely very
interesting. I would be interested to review your work for inclusion
in GSL Shell. Other people already expressed interest for GSL shell in
the area of signal processing.

Best regards,
Francesco