Re: Where a generalized Richardson extrapolation routine would fit in GSL?

[Rhys Ulerich <rhys dot ulerich at gmail dot com>]

> Also, are the "exact/normtable" arguments essential or just for
> convenience? If they can be computed easily by the user it would be
> good to have only a minimal number of arguments in the base function.

They're just for convenience.  Computing the norms doesn't require
difficult coding, but doing it in a storage and efficient way for the
vector case requires computing those norms during the extrapolation
process. It isn't difficult but it does require some index futzing.

> Ideally I'd like to start with a scalar version using normal C arrays,
> similar to the gsl_sum functions, for implementation simplicity.

For API simplicity?  Or simplicity of the underlying implementation?
There's just enough to the argument processing and general normtable
handling that it's easier to get good test coverage if there's one
implementation that does the whole enchilada and then everything else
sits as convenience wrappers.  That prevents changes to one
convenience function from accidentally changing the interface's
semantics (especially for the 'k' handling).

All that said, I'm happy to make the interface conform.  Would you be
willing to take draft public API in gsl_extrap.h (inline below) and
change the functions and their signatures to resemble what you'd like?
 That'll save me from needing to update the unit tests and
documentation as we iterate to convergence.

- Rhys

The declarations from that last patch look like:

int
gsl_extrap_richardson_vector_step(
        gsl_vector * Ah,
        const gsl_vector * Aht,
        double t,
        double ki);

int
gsl_extrap_richardson_step(
        double * Ah,
        const double * Aht,
        double t,
        double ki);

int
gsl_extrap_richardson_vector(
        gsl_matrix * A,
        double t,
        const gsl_vector * k,
        gsl_matrix * normtable,
        const gsl_vector * exact);

int
gsl_extrap_richardson(
        gsl_vector * A,
        double t,
        const gsl_vector * k,
        gsl_matrix * normtable,
        double exact);