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

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

[Rhys Ulerich]

> I've got a generalized Richardson extrapolation .m routine [1] that I'm going to port to C/GSL.  I'd like to contribute the finished code, but I'm not sure where in the GSL the functionality would fit.  I've looked but have found no similar sequence acceleration code.
>
> Does any place in the source tree jump out at anyone?
>
> [1] http://www.mathworks.com/matlabcentral/fileexchange/24388

The routine I mentioned in my note from August has been complete for
quite some time now.  I'd like to contribute it to GSL, but I'm still
unsure where it fits within GSL's build/layout system.  Any
suggestions for where it should go?

- Rhys

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

[Brian Gough]

At Sun, 29 Nov 2009 17:38:45 -0600,
Rhys Ulerich wrote:
> The routine I mentioned in my note from August has been complete for
> quite some time now.  I'd like to contribute it to GSL, but I'm still
> unsure where it fits within GSL's build/layout system.  Any
> suggestions for where it should go?

The closest routines are in sum/ but that's not an ideal name so we
can just create a new extrap/ directory for it.

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

[Rhys Ulerich]

[-- Attachment #1: Type: text/plain, Size: 494 bytes --]

Attached is a patch adding Richardson extrapolation routines, tests,
documentation, and an example.  Please let me know if I can answer any
questions.

> The closest routines are in sum/ but that's not an ideal name so we
> can just create a new extrap/ directory for it.

The routines sit in gsl_extrap.h and build under extrap/.  I included
the documentation under sum.texi as part of the Series Acceleration
chapter.  The technique's math fits there more closely than anywhere
else.

- Rhys

[-- Attachment #2: 0001-Added-Richardson-extrapolation-tests-and-docs.patch --]
[-- Type: application/octet-stream, Size: 55292 bytes --]

From f8f191fc06ddbc2a3f951568a13beb3f975fa41b Mon Sep 17 00:00:00 2001
From: Rhys Ulerich <rhys.ulerich@gmail.com>
Date: Sun, 13 Dec 2009 17:45:05 -0600
Subject: [PATCH] Added Richardson extrapolation, tests, and docs

Created build directory extrap/
Added Richardson extrapolation header, routines, and test cases
Added documentation to sum.texi, including example code
---
 Makefile.am                 |    4 +-
 configure.ac                |    2 +-
 doc/Makefile.am             |    4 +-
 doc/examples/richardson.c   |   84 ++++++++
 doc/examples/richardson.out |   17 ++
 doc/sum.texi                |  179 +++++++++++++++--
 extrap/ChangeLog            |    3 +
 extrap/Makefile.am          |   20 ++
 extrap/demo_richardson.c    |   84 ++++++++
 extrap/gsl_extrap.h         |   56 +++++
 extrap/richardson.c         |  184 +++++++++++++++++
 extrap/test_richardson.c    |  474 +++++++++++++++++++++++++++++++++++++++++++
 12 files changed, 1087 insertions(+), 24 deletions(-)
 create mode 100644 doc/examples/richardson.c
 create mode 100644 doc/examples/richardson.out
 create mode 100644 extrap/ChangeLog
 create mode 100644 extrap/Makefile.am
 create mode 100644 extrap/demo_richardson.c
 create mode 100644 extrap/gsl_extrap.h
 create mode 100644 extrap/richardson.c
 create mode 100644 extrap/test_richardson.c

diff --git a/Makefile.am b/Makefile.am
index 45b7c15..a512604 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -2,9 +2,9 @@
 
 # AUTOMAKE_OPTIONS = readme-alpha
 
-SUBDIRS = gsl utils sys test err const complex cheb block vector matrix permutation combination multiset sort ieee-utils cblas blas linalg eigen specfunc dht qrng rng randist fft poly fit multifit statistics siman sum integration interpolation histogram ode-initval roots multiroots min multimin monte ntuple diff deriv cdf wavelet bspline doc
+SUBDIRS = gsl utils sys test err const complex cheb block vector matrix permutation combination multiset sort ieee-utils cblas blas linalg eigen specfunc dht qrng rng randist fft poly fit multifit statistics siman sum integration interpolation histogram ode-initval roots multiroots min multimin monte ntuple diff deriv cdf wavelet bspline extrap doc
 
-SUBLIBS = block/libgslblock.la blas/libgslblas.la bspline/libgslbspline.la complex/libgslcomplex.la cheb/libgslcheb.la dht/libgsldht.la diff/libgsldiff.la deriv/libgslderiv.la eigen/libgsleigen.la err/libgslerr.la fft/libgslfft.la fit/libgslfit.la histogram/libgslhistogram.la ieee-utils/libgslieeeutils.la integration/libgslintegration.la interpolation/libgslinterpolation.la linalg/libgsllinalg.la matrix/libgslmatrix.la min/libgslmin.la monte/libgslmonte.la multifit/libgslmultifit.la multimin/libgslmultimin.la multiroots/libgslmultiroots.la ntuple/libgslntuple.la ode-initval/libgslodeiv.la permutation/libgslpermutation.la combination/libgslcombination.la multiset/libgslmultiset.la poly/libgslpoly.la qrng/libgslqrng.la randist/libgslrandist.la rng/libgslrng.la roots/libgslroots.la siman/libgslsiman.la sort/libgslsort.la specfunc/libgslspecfunc.la statistics/libgslstatistics.la sum/libgslsum.la sys/libgslsys.la test/libgsltest.la utils/libutils.la vector/libgslvector.la cdf/libgslcdf.la wavelet/libgslwavelet.la
+SUBLIBS = block/libgslblock.la blas/libgslblas.la bspline/libgslbspline.la complex/libgslcomplex.la cheb/libgslcheb.la dht/libgsldht.la diff/libgsldiff.la deriv/libgslderiv.la eigen/libgsleigen.la err/libgslerr.la fft/libgslfft.la fit/libgslfit.la histogram/libgslhistogram.la ieee-utils/libgslieeeutils.la integration/libgslintegration.la interpolation/libgslinterpolation.la linalg/libgsllinalg.la matrix/libgslmatrix.la min/libgslmin.la monte/libgslmonte.la multifit/libgslmultifit.la multimin/libgslmultimin.la multiroots/libgslmultiroots.la ntuple/libgslntuple.la ode-initval/libgslodeiv.la permutation/libgslpermutation.la combination/libgslcombination.la multiset/libgslmultiset.la poly/libgslpoly.la qrng/libgslqrng.la randist/libgslrandist.la rng/libgslrng.la roots/libgslroots.la siman/libgslsiman.la sort/libgslsort.la specfunc/libgslspecfunc.la statistics/libgslstatistics.la sum/libgslsum.la sys/libgslsys.la test/libgsltest.la utils/libutils.la vector/libgslvector.la cdf/libgslcdf.la wavelet/libgslwavelet.la extrap/libgslextrap.la
 
 pkginclude_HEADERS = gsl_math.h gsl_pow_int.h gsl_nan.h gsl_machine.h gsl_mode.h gsl_precision.h gsl_types.h gsl_version.h gsl_minmax.h gsl_inline.h
 
diff --git a/configure.ac b/configure.ac
index f9f5b04..6915dea 100644
--- a/configure.ac
+++ b/configure.ac
@@ -543,5 +543,5 @@ AH_VERBATIM([GSL_DISABLE_DEPRECATED],
 #define GSL_DISABLE_DEPRECATED 1])
 
 dnl
-AC_CONFIG_FILES([gsl_version.h gsl.spec gsl/Makefile test/Makefile err/Makefile sys/Makefile utils/Makefile const/Makefile min/Makefile multimin/Makefile ieee-utils/Makefile fft/Makefile specfunc/Makefile dht/Makefile fit/Makefile multifit/Makefile bspline/Makefile statistics/Makefile sum/Makefile roots/Makefile multiroots/Makefile ntuple/Makefile poly/Makefile qrng/Makefile rng/Makefile randist/Makefile siman/Makefile integration/Makefile interpolation/Makefile doc/Makefile block/Makefile vector/Makefile matrix/Makefile histogram/Makefile monte/Makefile ode-initval/Makefile cblas/Makefile blas/Makefile linalg/Makefile eigen/Makefile permutation/Makefile combination/Makefile multiset/Makefile sort/Makefile complex/Makefile diff/Makefile deriv/Makefile cheb/Makefile cdf/Makefile wavelet/Makefile Makefile])
+AC_CONFIG_FILES([gsl_version.h gsl.spec gsl/Makefile test/Makefile err/Makefile sys/Makefile utils/Makefile const/Makefile min/Makefile multimin/Makefile ieee-utils/Makefile fft/Makefile specfunc/Makefile dht/Makefile fit/Makefile multifit/Makefile bspline/Makefile statistics/Makefile sum/Makefile roots/Makefile multiroots/Makefile ntuple/Makefile poly/Makefile qrng/Makefile rng/Makefile randist/Makefile siman/Makefile integration/Makefile interpolation/Makefile doc/Makefile block/Makefile vector/Makefile matrix/Makefile histogram/Makefile monte/Makefile ode-initval/Makefile cblas/Makefile blas/Makefile linalg/Makefile eigen/Makefile permutation/Makefile combination/Makefile multiset/Makefile sort/Makefile complex/Makefile diff/Makefile deriv/Makefile cheb/Makefile cdf/Makefile wavelet/Makefile extrap/Makefile Makefile])
 AC_OUTPUT
diff --git a/doc/Makefile.am b/doc/Makefile.am
index c172250..a20c111 100644
--- a/doc/Makefile.am
+++ b/doc/Makefile.am
@@ -8,9 +8,9 @@ man_MANS = gsl.3 gsl-config.1 gsl-randist.1 gsl-histogram.1
 
 figures = multimin.eps siman-test.eps siman-energy.eps 12-cities.eps initial-route.eps final-route.eps fft-complex-radix2-f.eps fft-complex-radix2-t.eps fft-complex-radix2.eps fft-real-mixedradix.eps roots-bisection.eps roots-false-position.eps roots-newtons-method.eps roots-secant-method.eps histogram.eps histogram2d.eps min-interval.eps fit-wlinear.eps fit-wlinear2.eps fit-exp.eps ntuple.eps qrng.eps cheb.eps vdp.eps interp2.eps rand-beta.tex rand-binomial.tex rand-cauchy.tex rand-chisq.tex rand-erlang.tex rand-exponential.tex rand-fdist.tex rand-flat.tex rand-gamma.tex rand-gaussian.tex rand-geometric.tex rand-laplace.tex rand-logarithmic.tex rand-logistic.tex rand-lognormal.tex rand-pareto.tex rand-poisson.tex rand-hypergeometric.tex rand-nbinomial.tex rand-pascal.tex rand-bivariate-gaussian.tex rand-rayleigh.tex rand-rayleigh-tail.tex rand-tdist.tex rand-weibull.tex random-walk.tex randplots.gnp rand-exppow.tex rand-levy.tex rand-levyskew.tex rand-gumbel.tex rand-bernoulli.tex rand-gaussian-tail.tex rand-gumbel1.tex rand-gumbel2.tex landau.dat rand-landau.tex dwt-orig.eps dwt-samp.eps interpp2.eps bspline.eps
 
-examples_src = examples/blas.c examples/block.c examples/cblas.c examples/cdf.c examples/cheb.c examples/combination.c examples/multiset.c examples/const.c examples/demo_fn.c examples/diff.c examples/eigen.c examples/fft.c examples/fftmr.c examples/fftreal.c examples/fitting.c examples/fitting2.c examples/fitting3.c examples/histogram.c examples/histogram2d.c examples/ieee.c examples/ieeeround.c examples/integration.c examples/interp.c examples/intro.c examples/linalglu.c examples/matrix.c examples/matrixw.c examples/min.c examples/monte.c examples/ntupler.c examples/ntuplew.c examples/ode-initval.c examples/odefixed.c examples/permseq.c examples/permshuffle.c examples/polyroots.c examples/qrng.c examples/randpoisson.c examples/randwalk.c examples/rng.c examples/rngunif.c examples/rootnewt.c examples/roots.c examples/siman.c examples/sortsmall.c examples/specfun.c examples/specfun_e.c examples/stat.c examples/statsort.c examples/sum.c examples/vector.c examples/vectorr.c examples/vectorview.c examples/vectorw.c examples/demo_fn.h examples/dwt.c examples/expfit.c examples/nlfit.c examples/interpp.c examples/eigen_nonsymm.c examples/bspline.c examples/multimin.c examples/multiminfn.c examples/nmsimplex.c
+examples_src = examples/blas.c examples/block.c examples/cblas.c examples/cdf.c examples/cheb.c examples/combination.c examples/multiset.c examples/const.c examples/demo_fn.c examples/diff.c examples/eigen.c examples/fft.c examples/fftmr.c examples/fftreal.c examples/fitting.c examples/fitting2.c examples/fitting3.c examples/histogram.c examples/histogram2d.c examples/ieee.c examples/ieeeround.c examples/integration.c examples/interp.c examples/intro.c examples/linalglu.c examples/matrix.c examples/matrixw.c examples/min.c examples/monte.c examples/ntupler.c examples/ntuplew.c examples/ode-initval.c examples/odefixed.c examples/permseq.c examples/permshuffle.c examples/polyroots.c examples/qrng.c examples/randpoisson.c examples/randwalk.c examples/rng.c examples/rngunif.c examples/rootnewt.c examples/roots.c examples/siman.c examples/sortsmall.c examples/specfun.c examples/specfun_e.c examples/stat.c examples/statsort.c examples/sum.c examples/vector.c examples/vectorr.c examples/vectorview.c examples/vectorw.c examples/demo_fn.h examples/dwt.c examples/expfit.c examples/nlfit.c examples/interpp.c examples/eigen_nonsymm.c examples/bspline.c examples/multimin.c examples/multiminfn.c examples/nmsimplex.c examples/richardson.c
 
-examples_out = examples/blas.out examples/block.out examples/cblas.out examples/cdf.out examples/combination.out examples/multiset.out examples/const.out examples/diff.out examples/integration.out examples/intro.out examples/linalglu.out examples/min.out examples/polyroots.out examples/randpoisson.2.out examples/randpoisson.out examples/rng.out examples/rngunif.2.out examples/rngunif.out examples/sortsmall.out examples/specfun.out examples/specfun_e.out examples/stat.out examples/statsort.out examples/sum.out examples/vectorview.out examples/ecg.dat examples/dwt.dat examples/multimin.out examples/nmsimplex.out
+examples_out = examples/blas.out examples/block.out examples/cblas.out examples/cdf.out examples/combination.out examples/multiset.out examples/const.out examples/diff.out examples/integration.out examples/intro.out examples/linalglu.out examples/min.out examples/polyroots.out examples/randpoisson.2.out examples/randpoisson.out examples/rng.out examples/rngunif.2.out examples/rngunif.out examples/sortsmall.out examples/specfun.out examples/specfun_e.out examples/stat.out examples/statsort.out examples/sum.out examples/vectorview.out examples/ecg.dat examples/dwt.dat examples/multimin.out examples/nmsimplex.out examples/richardson.out
 
 noinst_DATA = $(examples_src) $(examples_out) $(figures)
 
diff --git a/doc/examples/richardson.c b/doc/examples/richardson.c
new file mode 100644
index 0000000..caa76bc
--- /dev/null
+++ b/doc/examples/richardson.c
@@ -0,0 +1,84 @@
+#include <stdio.h>
+#include <math.h>
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_extrap.h>
+
+const double h = 0.4;
+const double t = 2.00;
+const double x = 0.75;
+
+/* Method has truncation error like O(h^2) + O(h^4) + ... */
+double estimate(double h)
+{
+  return (sin(x + h) - sin(x - h)) / (2*h);
+}
+
+int
+main (void)
+{
+  printf("\nPerforming a single Richardson extrapolation step...\n");
+  {
+    gsl_vector *Ah  = gsl_vector_alloc(1);
+    gsl_vector *Aht = gsl_vector_alloc(1);
+
+    gsl_vector_set(Ah,  0, estimate(h));
+    gsl_vector_set(Aht, 0, estimate(h / t));
+
+    printf("Exact value for (d/dx) sin(%4f) =  %14.12f\n", x, cos(x));
+    printf("Approximate value for h   = %4f is %14.12f with error %14e\n",
+           h, gsl_vector_get(Ah, 0), fabs(cos(x) - gsl_vector_get(Ah, 0)));
+    printf("Approximate value for h/t = %4f is %14.12f with error %14e\n",
+           h / t, gsl_vector_get(Aht, 0), fabs(cos(x) - gsl_vector_get(Aht, 0)));
+
+    gsl_extrap_richardson_step(Ah, Aht, t, 2.0);
+
+    printf("After one extrapolation step value is   %14.12f with error %14e\n",
+           gsl_vector_get(Ah, 0), fabs(cos(x) - gsl_vector_get(Ah, 0)));
+
+    gsl_vector_free(Aht);
+    gsl_vector_free(Ah);
+  }
+
+  printf("\nPerforming iterated Richardson extrapolation...\n");
+  {
+    gsl_matrix *A         = gsl_matrix_alloc(1, 4);
+    gsl_vector *k         = gsl_vector_alloc(2);
+    gsl_matrix *normtable = gsl_matrix_alloc(A->size2, A->size2);
+    gsl_vector *exact     = gsl_vector_alloc(1);
+    int i, j;
+
+    gsl_vector_set(exact, 0, cos(x));
+
+    printf("Computing %d initial estimates of (d/dx) sin(%f)\n",
+           (int) A->size2, x);
+    for (i = 0; i < A->size2; ++i)
+      {
+        gsl_matrix_set(A, 0, i, estimate(h / pow(t, i)));
+      }
+
+    printf("Using leading error terms like O(h^2), O(h^4), O(h^6), O(h^8)...\n");
+    gsl_vector_set(k, 0, 2);
+    gsl_vector_set(k, 1, 4);
+
+    gsl_extrap_richardson(A, t, k, normtable, exact);
+
+    printf("The l_2 error after each iterated extrapolation step is:\n");
+    for (i = 0; i < A->size2; ++i)
+      {
+        for (j = 0; j < A->size2; ++j)
+          {
+            printf(" %16e", gsl_matrix_get(normtable, i, j));
+          }
+        printf("\n");
+      }
+    printf("\nThe final estimate is %14.12f\n", gsl_matrix_get(A, 0, 0));
+
+    gsl_vector_free(exact);
+    gsl_matrix_free(normtable);
+    gsl_vector_free(k);
+    gsl_matrix_free(A);
+  }
+
+  return 0;
+}
diff --git a/doc/examples/richardson.out b/doc/examples/richardson.out
new file mode 100644
index 0000000..b055ff2
--- /dev/null
+++ b/doc/examples/richardson.out
@@ -0,0 +1,17 @@
+
+Performing a single Richardson extrapolation step...
+Exact value for (d/dx) sin(0.750000) =  0.731688868874
+Approximate value for h   = 0.400000 is 0.712332666006 with error   1.935620e-02
+Approximate value for h/t = 0.200000 is 0.726820689647 with error   4.868179e-03
+After one extrapolation step value is   0.731650030860 with error   3.883801e-05
+
+Performing iterated Richardson extrapolation...
+Computing 4 initial estimates of (d/dx) sin(0.750000)
+Using leading error terms like O(h^2), O(h^4), O(h^6), O(h^8)...
+The l_2 error after each iterated extrapolation step is:
+     1.935620e-02              nan              nan              nan
+     4.868179e-03     3.883801e-05              nan              nan
+     1.218872e-03     2.436061e-06     9.264229e-09              nan
+     3.048323e-04     1.523898e-07     1.450697e-10     3.211875e-13
+
+The final estimate is 0.731688868873
diff --git a/doc/sum.texi b/doc/sum.texi
index f06dc44..4c920b3 100644
--- a/doc/sum.texi
+++ b/doc/sum.texi
@@ -4,21 +4,34 @@
 @cindex u-transform for series
 @cindex Levin u-transform
 @cindex convergence, accelerating a series
+@cindex extrapolation, Richardson
+@cindex Richardson extrapolation
 
-The functions described in this chapter accelerate the convergence of a
-series using the Levin @math{u}-transform.  This method takes a small number of
-terms from the start of a series and uses a systematic approximation to
-compute an extrapolated value and an estimate of its error.  The
+The functions described in this chapter accelerate the convergence of a series.
+Two methods are provided.
+
+The first method uses the Levin @math{u}-transform.  This method takes a small
+number of terms from the start of a series and uses a systematic approximation
+to compute an extrapolated value and an estimate of its error.  The
 @math{u}-transform works for both convergent and divergent series, including
 asymptotic series.
 
-These functions are declared in the header file @file{gsl_sum.h}.
+The second method is Richardson extrapolation.  This method combines estimates
+with known leading order errors to produce a result with higher order accuracy
+than any of the estimates.  It is well suited to problems where Taylor
+expansion-based error estimates are available, for example, finite difference
+techniques.
+
+These functions are declared in the header files @file{gsl_sum.h} and
+@file{gsl_extrap.h}.
 
 @menu
-* Acceleration functions::      
-* Acceleration functions without error estimation::  
-* Example of accelerating a series::  
-* Series Acceleration References::  
+* Acceleration functions::
+* Acceleration functions without error estimation::
+* Example of accelerating a series::
+* Richardson extrapolation::
+* Example of Richardson extrapolation::
+* Series Acceleration References::
 @end menu
 
 @node Acceleration functions
@@ -26,7 +39,7 @@ These functions are declared in the header file @file{gsl_sum.h}.
 
 The following functions compute the full Levin @math{u}-transform of a series
 with its error estimate.  The error estimate is computed by propagating
-rounding errors from each term through to the final extrapolation. 
+rounding errors from each term through to the final extrapolation.
 
 These functions are intended for summing analytic series where each term
 is known to high accuracy, and the rounding errors are assumed to
@@ -60,7 +73,7 @@ estimate of the absolute error stored in @var{abserr}.  The actual
 term-by-term sum is returned in @code{w->sum_plain}. The algorithm
 calculates the truncation error (the difference between two successive
 extrapolations) and round-off error (propagated from the individual
-terms) to choose an optimal number of terms for the extrapolation.  
+terms) to choose an optimal number of terms for the extrapolation.
 All the terms of the series passed in through @var{array} should be non-zero.
 @end deftypefun
 
@@ -111,7 +124,7 @@ truncation error, smoothing out any fluctuations.
 
 
 @node Example of accelerating a series
-@section Examples
+@section Example of accelerating a series
 
 The following code calculates an estimate of @math{\zeta(2) = \pi^2 / 6}
 using the series,
@@ -138,28 +151,120 @@ summation of the series converge slowly.
 @end example
 
 @noindent
-The output below shows that the Levin @math{u}-transform is able to obtain an 
-estimate of the sum to 1 part in 
+The output below shows that the Levin @math{u}-transform is able to obtain an
+estimate of the sum to 1 part in
 @c{$10^{10}$}
 @math{10^10} using the first eleven terms of the series.  The
 error estimate returned by the function is also accurate, giving
-the correct number of significant digits. 
+the correct number of significant digits.
 
 @example
-$ ./a.out 
+$ ./a.out
 @verbatiminclude examples/sum.out
 @end example
 
 @noindent
-Note that a direct summation of this series would require 
+Note that a direct summation of this series would require
 @c{$10^{10}$}
-@math{10^10} terms to achieve the same precision as the accelerated 
+@math{10^10} terms to achieve the same precision as the accelerated
 sum does in 13 terms.
 
+@node Richardson extrapolation
+@section Richardson extrapolation
+
+Richardson extrapolation estimates some value @math{A = lim A(h)} theoretically
+obtainable as @math{h} goes to @math{0}.  Here @math{A(h)} is an estimate of
+order @math{h^{k_i}} such that @math{A - A(h) = a_i*h^(k_i) + O(h^(k_(i+1)))}
+where @math{a_i} is nonzero and @math{m > n}.    Practically, one may not be
+able to compute this limit directly because it either requires large
+computational cost or because numerical truncation error dominates as @math{h}
+becomes small.  The technique performs series acceleration in the sense that
+from approximations @math{A_i(h)}, @math{A_i(h/t)}, @math{A_i(h/t^2)},
+@math{...} one may obtain an estimate of higher order than any known estimate
+in the sequence.
+
+In a Richardson extrapolation step two estimates @math{A_i(h)} and
+@math{A_i(h/t)} are combined to find an estimate of order @math{h^{k_(i+1)}}
+according to @math{A_(i+1)(h) = A_i(h/t) + (t^n - 1)^(-1)*(A_i(h/t) - A_i(h))}.
+Note that this did not require knowing the leading coefficient @math{a_i}, only
+the leading error order @math{k_i}.  One such step may be computed using
+@code{gsl_extrap_richardson_step}.  The process can be iterated to combine
+@math{A_(i+1)(h)} and @math{A_(i+1)(h/t)} to obtain @math{A_(i+2)(h)}.  The
+iterated process, as well as convenience mechanism to compute error norms given
+a known solution, can be performed using using @code{gsl_extrap_richardson}.
+
+@deftypefun {int} gsl_extrap_richardson_step ( gsl_vector * @var{Ah}, const gsl_vector * @var{Aht}, const double @var{t}, const double @var{ki})
+Given @math{A_i(h)}, an approximation of @math{A} such that @math{A - A(h) =
+a_i*h^k_i + a_(i+1)*h^k_(i+1) + ...} use Richardson extrapolation to
+find an approximation to leading order @math{k_(i+1)} from evaluations at
+@math{A_i(h)} and @math{A_i(h/t)} for @math{t > 0}.
+
+On entry, @var{Ah} contains the coarser approximation @math{A_i(h)} and
+@var{Aht} contains the finer approximation @math{A_i(h/t)}.  @var{ki} is
+@math{k_i}, the leading term order in the Tayor expansion, and @var{t} is the
+refinement factor between the two approximations.
+
+On success, the method returns @code{GSL_SUCCESS}.  On successful exit,
+@var{Ah} contains the extrapolated approximation @math{A_(i+1)(h)}.  On error,
+the method calls @code{gsl_error} and returns on of GSL's error constants.  The
+parameter @var{Ah} will be in an undefined state after any error.
+@end deftypefun
+
+@deftypefun {int} gsl_extrap_richardson ( gsl_matrix * const @var{A}, const double @var{t}, const gsl_vector * @var{k}, gsl_matrix * @var{normtable}, const gsl_vector * const @var{exact})
+Perform Richardson extrapolation on a sequence of refined approximations.
+The routine can perform multiple step extrapolation, e.g. using
+@math{A_0(h)}, @math{A_0(h/2)}, and @math{A_0(h/4)} to compute @math{A_2(h)}.
+See @code{gsl_extrap_richardson_step} for the terminology in use.
+
+On entry, column @math{i} of @var{A} contains @math{A_0(h/t^i)} while @var{t}
+contains the refinement factor between each pair of approximations.  The
+refinement factor must be fixed for all columns.  The leading error term orders
+are provided in @var{k}.  They are assumed to start at 1 if @var{k} is
+@code{NULL}.  They are assumed to increment by 1 if no second element is
+provided.  If the vector @var{k} is shorter than @code{A->size2}, any
+unspecified higher index entries are assumed to grow by
+@code{gsl_vector_get(k,k->size-1) - gsl_vector_get(k,k->size-2)}.
+
+If @var{normtable} is non-@code{NULL} on entry, the routine produces in
+@var{normtable} a table showing the @math{l_2} error at each step in the
+extrapolation process calculated against the provided exact solution in
+@var{exact}.  If @var{exact} is @code{NULL} then it is treated as as the zero
+vector and the resulting @var{normtable} simply contains the @math{l_2} norm of
+each extrapolation step.
+
+On success, the method returns @code{GSL_SUCCESS}.  On successful exit, column
+0 of @var{A} contains @math{A_(A->size2)(h)} and all other columns will have
+been overwritten with temporary results.  If provided, @var{normtable} will be
+modified as described above.  On error the method calls @code{gsl_error} and
+returns one of GSL's error codes.  The parameters @var{A} and @var{normtable}
+will be in an undefined state after any error.
+@end deftypefun
+
+@node Example of Richardson extrapolation
+@section Example of Richardson extrapolation
+
+Suppose, hypothetically, you wanted to compute the derivative of @math{sin(x)}
+accurate to a high order but you only knew the central difference formula
+@math{f'(x) = (f(x+h) - f(x-h))/(2*h) + O(h^2) + O(h^4) + ...}.  Using
+evaluations for @math{h}, @math{h/2}, @math{h/2^2}, @math{...} you can obtain a
+much higher order estimate of the derivative:
+
+@example
+@verbatiminclude examples/richardson.c
+@end example
+
+The following output is produced:
+
+@example
+$ ./a.out
+@verbatiminclude examples/richardson.out
+@end example
+
 @node Series Acceleration References
 @section References and Further Reading
 
-The algorithms used by these functions are described in the following papers,
+The algorithms used by the @math{u}-transform functions are described in
+the following papers,
 
 @itemize @asis
 @item
@@ -186,3 +291,39 @@ A review paper on the Levin Transform is available online,
 Herbert H. H. Homeier, Scalar Levin-Type Sequence Transformations,
 @uref{http://arxiv.org/abs/math/0005209}.
 @end itemize
+
+@noindent
+Information on Richardson extrapolation can be found in nearly any
+introductory numerical analysis text or online.
+One text providing information is
+@itemize asis
+@item
+David Kincaid and Ward Cheney, @cite{Numerical Analysis: Mathematics of
+Scientific Computing} (2009), AMS, ISBN 9780821847886.
+@end itemize
+
+@noindent
+The brief introduction given in this manual uses
+the following Wikipedia both for content and to fix notation:
+@itemize asis
+@item
+Wikipedia: Richardson extrapolation,
+@uref{http://en.wikipedia.org/wiki/Richardson_extrapolation}.
+@end itemize
+
+@noindent
+The original work by Richardson appeared in
+@itemize asis
+@item
+Richardson, L. F.,
+"The approximate arithmetical solution by finite differences of physical problems including differential equations, with an application to the stresses in a masonry dam",
+@cite{Philosophical Transactions of the Royal Society of London},
+Series A 210: 307--357, 1911,
+@uref{doi:10.1098/rsta.1911.0009}.
+@item
+Richardson, L. F.,
+"The deferred approach to the limit",
+@cite{Philosophical Transactions of the Royal Society of London},
+Series A 226: 299--349, 1927,
+@uref{doi:10.1098/rsta.1927.0008}.
+@end itemize
diff --git a/extrap/ChangeLog b/extrap/ChangeLog
new file mode 100644
index 0000000..9183ab4
--- /dev/null
+++ b/extrap/ChangeLog
@@ -0,0 +1,3 @@
+2009-12-13  Rhys Ulerich  <rhys.ulerich@gmail.com>
+
+	* added build structure and Richardson extrapolation
diff --git a/extrap/Makefile.am b/extrap/Makefile.am
new file mode 100644
index 0000000..855880d
--- /dev/null
+++ b/extrap/Makefile.am
@@ -0,0 +1,20 @@
+noinst_LTLIBRARIES = libgslextrap.la
+
+pkginclude_HEADERS = gsl_extrap.h
+
+INCLUDES = -I$(top_srcdir)
+
+libgslextrap_la_SOURCES = richardson.c
+
+noinst_HEADERS =
+
+TESTS = $(check_PROGRAMS)
+check_PROGRAMS = # Append below
+
+check_PROGRAMS += test_richardson
+test_richardson_SOURCES = test_richardson.c
+test_richardson_LDADD = libgslextrap.la ../vector/libgslvector.la  ../matrix/libgslmatrix.la  ../block/libgslblock.la ../ieee-utils/libgslieeeutils.la ../err/libgslerr.la ../blas/libgslblas.la ../cblas/libgslcblas.la ../test/libgsltest.la ../sys/libgslsys.la
+
+noinst_PROGRAMS = demo_richardson
+demo_richardson_SOURCES = demo_richardson.c
+demo_richardson_LDADD = libgslextrap.la ../vector/libgslvector.la  ../matrix/libgslmatrix.la  ../block/libgslblock.la ../ieee-utils/libgslieeeutils.la ../err/libgslerr.la ../blas/libgslblas.la ../cblas/libgslcblas.la ../test/libgsltest.la ../sys/libgslsys.la
diff --git a/extrap/demo_richardson.c b/extrap/demo_richardson.c
new file mode 100644
index 0000000..caa76bc
--- /dev/null
+++ b/extrap/demo_richardson.c
@@ -0,0 +1,84 @@
+#include <stdio.h>
+#include <math.h>
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_extrap.h>
+
+const double h = 0.4;
+const double t = 2.00;
+const double x = 0.75;
+
+/* Method has truncation error like O(h^2) + O(h^4) + ... */
+double estimate(double h)
+{
+  return (sin(x + h) - sin(x - h)) / (2*h);
+}
+
+int
+main (void)
+{
+  printf("\nPerforming a single Richardson extrapolation step...\n");
+  {
+    gsl_vector *Ah  = gsl_vector_alloc(1);
+    gsl_vector *Aht = gsl_vector_alloc(1);
+
+    gsl_vector_set(Ah,  0, estimate(h));
+    gsl_vector_set(Aht, 0, estimate(h / t));
+
+    printf("Exact value for (d/dx) sin(%4f) =  %14.12f\n", x, cos(x));
+    printf("Approximate value for h   = %4f is %14.12f with error %14e\n",
+           h, gsl_vector_get(Ah, 0), fabs(cos(x) - gsl_vector_get(Ah, 0)));
+    printf("Approximate value for h/t = %4f is %14.12f with error %14e\n",
+           h / t, gsl_vector_get(Aht, 0), fabs(cos(x) - gsl_vector_get(Aht, 0)));
+
+    gsl_extrap_richardson_step(Ah, Aht, t, 2.0);
+
+    printf("After one extrapolation step value is   %14.12f with error %14e\n",
+           gsl_vector_get(Ah, 0), fabs(cos(x) - gsl_vector_get(Ah, 0)));
+
+    gsl_vector_free(Aht);
+    gsl_vector_free(Ah);
+  }
+
+  printf("\nPerforming iterated Richardson extrapolation...\n");
+  {
+    gsl_matrix *A         = gsl_matrix_alloc(1, 4);
+    gsl_vector *k         = gsl_vector_alloc(2);
+    gsl_matrix *normtable = gsl_matrix_alloc(A->size2, A->size2);
+    gsl_vector *exact     = gsl_vector_alloc(1);
+    int i, j;
+
+    gsl_vector_set(exact, 0, cos(x));
+
+    printf("Computing %d initial estimates of (d/dx) sin(%f)\n",
+           (int) A->size2, x);
+    for (i = 0; i < A->size2; ++i)
+      {
+        gsl_matrix_set(A, 0, i, estimate(h / pow(t, i)));
+      }
+
+    printf("Using leading error terms like O(h^2), O(h^4), O(h^6), O(h^8)...\n");
+    gsl_vector_set(k, 0, 2);
+    gsl_vector_set(k, 1, 4);
+
+    gsl_extrap_richardson(A, t, k, normtable, exact);
+
+    printf("The l_2 error after each iterated extrapolation step is:\n");
+    for (i = 0; i < A->size2; ++i)
+      {
+        for (j = 0; j < A->size2; ++j)
+          {
+            printf(" %16e", gsl_matrix_get(normtable, i, j));
+          }
+        printf("\n");
+      }
+    printf("\nThe final estimate is %14.12f\n", gsl_matrix_get(A, 0, 0));
+
+    gsl_vector_free(exact);
+    gsl_matrix_free(normtable);
+    gsl_vector_free(k);
+    gsl_matrix_free(A);
+  }
+
+  return 0;
+}
diff --git a/extrap/gsl_extrap.h b/extrap/gsl_extrap.h
new file mode 100644
index 0000000..b15e97e
--- /dev/null
+++ b/extrap/gsl_extrap.h
@@ -0,0 +1,56 @@
+/* extrap/gsl_extrap.h
+ *
+ * Copyright (C) 2009 Rhys Ulerich
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 3 of the License, or (at
+ * your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+#ifndef __GSL_EXTRAP_H__
+#define __GSL_EXTRAP_H__
+
+#include <gsl/gsl_matrix_double.h>
+#include <gsl/gsl_vector_double.h>
+
+#undef __BEGIN_DECLS
+#undef __END_DECLS
+#ifdef __cplusplus
+# define __BEGIN_DECLS extern "C" {
+# define __END_DECLS }
+#else
+# define __BEGIN_DECLS /* empty */
+# define __END_DECLS /* empty */
+#endif
+
+__BEGIN_DECLS
+
+int
+gsl_extrap_richardson_step(
+        gsl_vector * Ah,
+        const gsl_vector * Aht,
+        const double t,
+        const double ki);
+
+int
+gsl_extrap_richardson(
+        gsl_matrix * const A,
+        const double t,
+        const gsl_vector * k,
+        gsl_matrix * normtable,
+        const gsl_vector * const exact);
+
+
+__END_DECLS
+
+#endif /* __GSL_EXTRAP_H__ */
diff --git a/extrap/richardson.c b/extrap/richardson.c
new file mode 100644
index 0000000..4b207cb
--- /dev/null
+++ b/extrap/richardson.c
@@ -0,0 +1,184 @@
+/* extrap/gsl_extrap.h
+ *
+ * Copyright (C) 2009 Rhys Ulerich
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 3 of the License, or (at
+ * your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+#include <gsl/gsl_blas.h>
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_nan.h>
+#include <gsl/gsl_errno.h>
+#include <gsl/gsl_extrap.h>
+#include <gsl/gsl_matrix_double.h>
+#include <gsl/gsl_vector_double.h>
+
+int
+gsl_extrap_richardson_step(
+  gsl_vector * Ah,
+  const gsl_vector * Aht,
+  const double t,
+  const double ki)
+{
+  if (ki == 0.0)
+    {
+      GSL_ERROR("ki == 0 invalid as it will cause a division-by-zero",
+                GSL_EDOM);
+    }
+
+  const double tki       = pow(t, ki);
+  const double inv_tkim1 = 1.0 / (tki - 1.0);
+
+  gsl_blas_dscal(-inv_tkim1, Ah);
+
+  int error = gsl_blas_daxpy(tki * inv_tkim1, Aht, Ah);
+  if (error) return error;
+
+  return GSL_SUCCESS;
+}
+
+int
+gsl_extrap_richardson(
+  gsl_matrix * const A,
+  const double t,
+  const gsl_vector * k,
+  gsl_matrix * normtable,
+  const gsl_vector * const exact)
+{
+  gsl_vector * scratch = NULL;
+  const size_t k_size = (k) ? k->size : 0;
+  int i, j;
+
+  if (exact)
+    {
+      if (!normtable)
+        {
+          GSL_ERROR("exact provided but normtable was not",
+                    GSL_EINVAL);
+        }
+      if (exact->size != A->size1)
+        {
+          GSL_ERROR("exact->size does not match A->size1",
+                    GSL_EINVAL);
+        }
+    }
+
+  if (normtable)
+    {
+      if (normtable->size1 != A->size2)
+        {
+          GSL_ERROR("normtable->size1 does not match A->size2",
+                    GSL_EINVAL);
+        }
+      if (normtable->size2 != A->size2)
+        {
+          GSL_ERROR("normtable->size2 does not match A->size2",
+                    GSL_EINVAL);
+        }
+
+      gsl_matrix_set_all(normtable, GSL_NAN);
+
+      /* Compute the first column in the norm table */
+      if (exact)
+        {
+          scratch = gsl_vector_alloc(A->size1);
+        }
+      for (i = 0; i < normtable->size2; ++i)
+        {
+          double norm;
+          if (exact)
+            {
+              gsl_matrix_get_col(scratch, A, i);
+              gsl_blas_daxpy(-1.0, exact, scratch);
+              norm = gsl_blas_dnrm2(scratch);
+            }
+          else
+            {
+              gsl_vector_view Ai = gsl_matrix_column(A, i);
+              norm = gsl_blas_dnrm2(&Ai.vector);
+            }
+          gsl_matrix_set(normtable, i, 0, norm);
+        }
+    }
+
+  for (i = 0; i < A->size2; ++i)
+    {
+
+      /* Provide automagic around the k parameter */
+      double ki;
+      if (k_size)
+        {
+          if (i < k_size)
+            {
+              ki = gsl_vector_get(k, i);
+            }
+          else if (k_size == 1)
+            {
+              ki = gsl_vector_get(k, 0) + i;
+            }
+          else
+            {
+              const double last_ki = gsl_vector_get(k, k_size - 1);
+              const double increment = last_ki - gsl_vector_get(k, k_size - 2);
+              ki = last_ki + (i - (k_size - 1)) * increment;
+            }
+        }
+      else
+        {
+          ki = i + 1;
+        }
+
+      /* Perform the extrapolation using leading term order ki */
+      for (j = 0; j < A->size2 - i - 1; ++j)
+        {
+          gsl_vector_view Aih  = gsl_matrix_column(A, j);
+          gsl_vector_view Aiht = gsl_matrix_column(A, j + 1);
+
+          const int error = gsl_extrap_richardson_step(
+                              &Aih.vector, &Aiht.vector, t, ki);
+          if (error)
+            {
+              if (scratch)
+                {
+                  gsl_vector_free(scratch);
+                }
+              return error;
+            }
+
+          if (normtable)
+            {
+              double norm;
+              if (exact)
+                {
+                  gsl_vector_memcpy(scratch, &Aih.vector);
+                  gsl_blas_daxpy(-1.0, exact, scratch);
+                  norm = gsl_blas_dnrm2(scratch);
+                }
+              else
+                {
+                  norm = gsl_blas_dnrm2(&Aih.vector);
+                }
+              gsl_matrix_set(normtable, i + j + 1, i + 1, norm);
+            }
+        }
+    }
+
+  if (scratch)
+    {
+      gsl_vector_free(scratch);
+    }
+
+  return GSL_SUCCESS;
+}
diff --git a/extrap/test_richardson.c b/extrap/test_richardson.c
new file mode 100644
index 0000000..903ccd7
--- /dev/null
+++ b/extrap/test_richardson.c
@@ -0,0 +1,474 @@
+#include <gsl/gsl_extrap.h>
+#include <gsl/gsl_ieee_utils.h>
+#include <gsl/gsl_machine.h>
+#include <gsl/gsl_sys.h>
+#include <gsl/gsl_test.h>
+
+void
+test_richardson_extrapolation_step()
+{
+    const int    ki           = 1;
+    const size_t n            = 1;
+    gsl_vector *Ah   = gsl_vector_alloc(n);
+    gsl_vector *Aht  = gsl_vector_alloc(n);
+
+    {
+        const double t = 2.0;
+        gsl_vector_set(Ah,  0, 1.0);
+        gsl_vector_set(Aht, 0, 2.0);
+        gsl_test(gsl_extrap_richardson_step(Ah, Aht, t, ki),
+                "Unexpected error reported in %s for t=%f", __func__, t);
+        gsl_test_abs(gsl_vector_get(Ah, 0), 3.0, GSL_DBL_EPSILON,
+                "Extrapolation in %s for ki=%d, t=%f", __func__, ki, t);
+    }
+
+    {
+        const double t = 3.0;
+        gsl_vector_set(Ah,  0, 1.0);
+        gsl_vector_set(Aht, 0, 2.0);
+        gsl_test(gsl_extrap_richardson_step(Ah, Aht, t, ki),
+                "Unexpected error reported in %s for t=%f", __func__, t);
+        gsl_test_abs(gsl_vector_get(Ah, 0), 5.0/2.0, GSL_DBL_EPSILON,
+                "Extrapolation in %s for ki=%d, t=%f", __func__, ki, t);
+    }
+
+    gsl_vector_free(Aht);
+    gsl_vector_free(Ah);
+}
+
+void
+test_richardson_extrapolation_defaults()
+{
+    gsl_matrix * data1 = gsl_matrix_alloc(1, 2);
+    gsl_matrix * data2 = gsl_matrix_alloc(data1->size1, data1->size2);
+    {
+        gsl_matrix_set(data1, 0, 0, 1.0);
+        gsl_matrix_set(data1, 0, 1, 2.0);
+
+        gsl_test(gsl_extrap_richardson(data1, 2, NULL, NULL, NULL),
+                "Unexpected error reported in %s");
+        gsl_test_abs(gsl_matrix_get(data1, 0, 0), 3.0, GSL_DBL_EPSILON,
+                "%s scalar correct result", __func__);
+    }
+
+    {
+        gsl_matrix_set(data1, 0, 0, 1.0);
+        gsl_matrix_set(data1, 0, 1, 2.0);
+
+        gsl_test(gsl_extrap_richardson(data1, 3, NULL, NULL, NULL),
+                "Unexpected error reported in %s");
+        gsl_test_abs(gsl_matrix_get(data1, 0, 0), 5.0/2.0, GSL_DBL_EPSILON,
+                "%s scalar correct result", __func__);
+    }
+
+    {
+        gsl_matrix_set(data1, 0, 0, 1.0);
+        gsl_matrix_set(data1, 0, 1, 2.0);
+        gsl_matrix_memcpy(data2, data1);
+
+        gsl_test(gsl_extrap_richardson(data1, 2, NULL, NULL, NULL),
+                "Unexpected error reported in %s");
+        const double result1 = gsl_matrix_get(data1, 0, 0);
+
+        gsl_test_abs(result1, 3.0, GSL_DBL_EPSILON,
+                "%s scalar correct result", __func__);
+
+        gsl_vector * k = gsl_vector_alloc(1);
+        gsl_vector_set(k, 0, 1);
+        gsl_test(gsl_extrap_richardson(data2, 2, k, NULL, NULL),
+                "Unexpected error reported in %s");
+        gsl_vector_free(k);
+        const double result2 = gsl_matrix_get(data2, 0, 0);
+
+        gsl_test_abs(result1, result2, GSL_DBL_EPSILON,
+                "%s with k == NULL", __func__);
+    }
+
+    {
+        int i, j;
+        gsl_matrix_set(data1, 0, 0, 1.0);
+        gsl_matrix_set(data1, 0, 1, 2.0);
+        gsl_matrix_memcpy(data2, data1);
+
+        gsl_matrix * normtable1 = gsl_matrix_alloc(data1->size2, data1->size2);
+        gsl_matrix * normtable2 = gsl_matrix_alloc(data1->size2, data1->size2);
+
+        gsl_test(gsl_extrap_richardson(
+                    data1, 2, NULL, normtable1, NULL),
+                "Unexpected error reported in %s");
+        const double result1 = gsl_matrix_get(data1, 0, 0);
+
+        gsl_test_abs(result1, 3.0, GSL_DBL_EPSILON,
+                "%s correct extrapolation answer");
+
+        for (i = 0; i < normtable1->size1 - 1; ++i) {
+            for (j = i+1; j < normtable1->size2; ++j) {
+                gsl_test(!gsl_isnan(gsl_matrix_get(normtable1, i, j)),
+                        "%s expected NAN in normtable1 at (%d, %d)",
+                        __func__, i, j);
+            }
+        }
+
+        gsl_test_abs(gsl_matrix_get(normtable1, 0, 0), 1.0, GSL_DBL_EPSILON,
+                "%s normtable1 (0,0) value", __func__);
+        gsl_test_abs(gsl_matrix_get(normtable1, 1, 0), 2.0, GSL_DBL_EPSILON,
+                "%s normtable1 (1,0) value", __func__);
+        gsl_test_abs(gsl_matrix_get(normtable1, 1, 1), 3.0, GSL_DBL_EPSILON,
+                "%s normtable1 (1,1) value", __func__);
+
+        gsl_vector * exact = gsl_vector_alloc(1);
+        gsl_vector_set_zero(exact);
+        gsl_test(gsl_extrap_richardson(
+                    data2, 2.0, NULL, normtable2, exact),
+                "Unexpected error reported in %s");
+        gsl_vector_free(exact);
+        const double result2 = gsl_matrix_get(data2, 0, 0);
+
+        for (i = 0; i < normtable2->size1 - 1; ++i) {
+            for (j = i+1; j < normtable2->size2; ++j) {
+                gsl_test(!gsl_isnan(gsl_matrix_get(normtable2, i, j)),
+                        "%s expected NAN in normtable2 at (%d, %d)",
+                        __func__, i, j);
+            }
+        }
+
+        gsl_test_abs(result1, result2, GSL_DBL_EPSILON,
+                "%s with exact == NULL", __func__);
+
+        for (i = 0; i < normtable1->size1; ++i) {
+            for (j = 0; j < i + 1; ++j) {
+                gsl_test_abs(gsl_matrix_get(normtable1, i, j),
+                             gsl_matrix_get(normtable2, i, j),
+                             GSL_DBL_EPSILON,
+                             "%s identical normtable results at (%d, %d)",
+                             __func__, i, j);
+            }
+        }
+
+        gsl_matrix_free(normtable2);
+        gsl_matrix_free(normtable1);
+    }
+
+    gsl_matrix_free(data2);
+    gsl_matrix_free(data1);
+}
+
+void
+test_richardson_extrapolation_twolevels()
+{
+    gsl_matrix * data2 = gsl_matrix_alloc(1, 2);
+    gsl_matrix * data3 = gsl_matrix_alloc(1, 3);
+    gsl_vector * k1 = gsl_vector_alloc(1);
+    gsl_vector * k2 = gsl_vector_alloc(2);
+
+    {
+        gsl_matrix_set(data2, 0, 0, 1.0);
+        gsl_matrix_set(data2, 0, 1, 2.0);
+        gsl_vector_set(k1, 0, 1);
+
+        gsl_test(gsl_extrap_richardson(data2, 2, k1, NULL, NULL),
+                "Unexpected error reported in %s");
+        gsl_test_abs(gsl_matrix_get(data2, 0, 0), 3.0, GSL_DBL_EPSILON,
+                "%s scalar correct result at %s:%d",
+                __func__, __FILE__, __LINE__);
+    }
+
+    {
+        gsl_matrix_set(data2, 0, 0, 2.0);
+        gsl_matrix_set(data2, 0, 1, 3.0);
+        gsl_vector_set(k1, 0, 1);
+
+        gsl_test(gsl_extrap_richardson(data2, 2, k1, NULL, NULL),
+                "Unexpected error reported in %s");
+        gsl_test_abs(gsl_matrix_get(data2, 0, 0), 4.0, GSL_DBL_EPSILON,
+                "%s scalar correct result at %s:%d",
+                __func__, __FILE__, __LINE__);
+    }
+
+    {
+        gsl_matrix_set(data2, 0, 0, 3.0);
+        gsl_matrix_set(data2, 0, 1, 4.0);
+        gsl_vector_set(k1, 0, 2);
+
+        gsl_test(gsl_extrap_richardson(data2, 2, k1, NULL, NULL),
+                "Unexpected error reported in %s");
+        gsl_test_abs(gsl_matrix_get(data2, 0, 0), 13.0/3.0, GSL_DBL_EPSILON,
+                "%s scalar correct result at %s:%d",
+                __func__, __FILE__, __LINE__);
+    }
+
+    {
+        gsl_matrix_set(data3, 0, 0, 1.0);
+        gsl_matrix_set(data3, 0, 1, 2.0);
+        gsl_matrix_set(data3, 0, 2, 3.0);
+        gsl_vector_set(k1, 0, 1);
+
+        gsl_test(gsl_extrap_richardson(data3, 2, k1, NULL, NULL),
+                "Unexpected error reported in %s");
+        gsl_test_abs(gsl_matrix_get(data3, 0, 0), 13.0/3.0, GSL_DBL_EPSILON,
+                "%s scalar correct result at %s:%d",
+                __func__, __FILE__, __LINE__);
+    }
+
+    {
+        int i, j;
+        gsl_matrix * normtable = gsl_matrix_alloc(data3->size2, data3->size2);
+
+        gsl_matrix_set(data3, 0, 0, 1.0);
+        gsl_matrix_set(data3, 0, 1, 2.0);
+        gsl_matrix_set(data3, 0, 2, 3.0);
+        gsl_vector_set(k2, 0, 1);
+        gsl_vector_set(k2, 1, 2);
+
+        gsl_test(gsl_extrap_richardson(data3, 2, k2, normtable, NULL),
+                "Unexpected error reported in %s");
+        gsl_test_abs(gsl_matrix_get(data3, 0, 0), 13.0/3.0, GSL_DBL_EPSILON,
+                "%s scalar correct result at %s:%d",
+                __func__, __FILE__, __LINE__);
+
+        for (i = 0; i < normtable->size1 - 1; ++i) {
+            for (j = i+1; j < normtable->size2; ++j) {
+                gsl_test(!gsl_isnan(gsl_matrix_get(normtable, i, j)),
+                        "%s expected NAN in normtable at (%d, %d)",
+                        __func__, i, j);
+            }
+        }
+
+        gsl_test_abs(gsl_matrix_get(normtable, 0, 0), 1.0, GSL_DBL_EPSILON,
+                "%s normtable (0,0) value", __func__);
+        gsl_test_abs(gsl_matrix_get(normtable, 1, 0), 2.0, GSL_DBL_EPSILON,
+                "%s normtable (1,0) value", __func__);
+        gsl_test_abs(gsl_matrix_get(normtable, 2, 0), 3.0, GSL_DBL_EPSILON,
+                "%s normtable (2,0) value", __func__);
+        gsl_test_abs(gsl_matrix_get(normtable, 1, 1), 3.0, GSL_DBL_EPSILON,
+                "%s normtable (1,1) value", __func__);
+        gsl_test_abs(gsl_matrix_get(normtable, 2, 1), 4.0, GSL_DBL_EPSILON,
+                "%s normtable (2,1) value", __func__);
+        gsl_test_abs(gsl_matrix_get(normtable, 2, 2), 13.0/3.0, GSL_DBL_EPSILON,
+                "%s normtable (2,2) value", __func__);
+
+        gsl_matrix_free(normtable);
+    }
+
+    {
+        gsl_matrix_set(data3, 0, 0, 1.0);
+        gsl_matrix_set(data3, 0, 1, 2.0);
+        gsl_matrix_set(data3, 0, 2, 3.0);
+        gsl_vector_set(k2, 0, 2);
+        gsl_vector_set(k2, 1, 3);
+
+        gsl_test(gsl_extrap_richardson(data3, 2, k2, NULL, NULL),
+                "Unexpected error reported in %s");
+        gsl_test_abs(gsl_matrix_get(data3, 0, 0), 73.0/21.0, GSL_DBL_EPSILON,
+                "%s scalar correct result at %s:%d",
+                __func__, __FILE__, __LINE__);
+    }
+
+    {
+        gsl_matrix_set(data3, 0, 0, 1.0);
+        gsl_matrix_set(data3, 0, 1, 2.0);
+        gsl_matrix_set(data3, 0, 2, 3.0);
+        gsl_vector_set(k1, 0, 2);
+
+        gsl_test(gsl_extrap_richardson(data3, 2, k1, NULL, NULL),
+                "Unexpected error reported in %s");
+        gsl_test_abs(gsl_matrix_get(data3, 0, 0), 73.0/21.0, GSL_DBL_EPSILON,
+                "%s scalar correct result at %s:%d",
+                __func__, __FILE__, __LINE__);
+    }
+
+
+    gsl_vector_free(k2);
+    gsl_vector_free(k1);
+    gsl_matrix_free(data3);
+    gsl_matrix_free(data2);
+
+}
+
+void
+test_richardson_extrapolation_multiplelevels()
+{
+    gsl_matrix * data2 = gsl_matrix_alloc(1, 2);
+    gsl_matrix * data4 = gsl_matrix_alloc(1, 4);
+    gsl_vector * k1 = gsl_vector_alloc(1);
+    gsl_vector * k2 = gsl_vector_alloc(2);
+    gsl_vector * k3 = gsl_vector_alloc(3);
+
+    {
+        int i, j;
+        gsl_vector_set(k1, 0, 3);
+
+        gsl_matrix_set(data2, 0, 0, 1.0);
+        gsl_matrix_set(data2, 0, 1, 2.0);
+        gsl_test(gsl_extrap_richardson(data2, 2, k1, NULL, NULL),
+                "Unexpected error reported in %s");
+        const double xx = gsl_matrix_get(data2, 0, 0);
+
+        gsl_matrix_set(data2, 0, 0, 2.0);
+        gsl_matrix_set(data2, 0, 1, 3.0);
+        gsl_test(gsl_extrap_richardson(data2, 2, k1, NULL, NULL),
+                "Unexpected error reported in %s");
+        const double xy = gsl_matrix_get(data2, 0, 0);
+
+        gsl_matrix_set(data2, 0, 0, 3.0);
+        gsl_matrix_set(data2, 0, 1, 4.0);
+        gsl_test(gsl_extrap_richardson(data2, 2, k1, NULL, NULL),
+                "Unexpected error reported in %s");
+        const double xz = gsl_matrix_get(data2, 0, 0);
+
+        gsl_vector_set(k1, 0, 6);
+
+        gsl_matrix_set(data2, 0, 0, xx);
+        gsl_matrix_set(data2, 0, 1, xy);
+        gsl_test(gsl_extrap_richardson(data2, 2, k1, NULL, NULL),
+                "Unexpected error reported in %s");
+        const double yx = gsl_matrix_get(data2, 0, 0);
+
+        gsl_matrix_set(data2, 0, 0, xy);
+        gsl_matrix_set(data2, 0, 1, xz);
+        gsl_test(gsl_extrap_richardson(data2, 2, k1, NULL, NULL),
+                "Unexpected error reported in %s");
+        const double yy = gsl_matrix_get(data2, 0, 0);
+
+        gsl_vector_set(k1, 0, 9);
+
+        gsl_matrix_set(data2, 0, 0, yx);
+        gsl_matrix_set(data2, 0, 1, yy);
+        gsl_test(gsl_extrap_richardson(data2, 2, k1, NULL, NULL),
+                "Unexpected error reported in %s");
+        const double zx = gsl_matrix_get(data2, 0, 0);
+
+        gsl_vector_set(k3, 0, 3);
+        gsl_vector_set(k3, 1, 6);
+        gsl_vector_set(k3, 2, 9);
+
+        gsl_matrix_set(data4, 0, 0, 1.0);
+        gsl_matrix_set(data4, 0, 1, 2.0);
+        gsl_matrix_set(data4, 0, 2, 3.0);
+        gsl_matrix_set(data4, 0, 3, 4.0);
+        gsl_test(gsl_extrap_richardson(data4, 2, k3, NULL, NULL),
+                "Unexpected error reported in %s");
+        gsl_test_abs(gsl_matrix_get(data4, 0, 0), zx, GSL_DBL_EPSILON,
+                "%s scalar correct result at %s:%d",
+                __func__, __FILE__, __LINE__);
+
+        gsl_matrix * normtable = gsl_matrix_alloc(data4->size2, data4->size2);
+
+        gsl_vector_set(k2, 0, 3);
+        gsl_vector_set(k2, 1, 6);
+
+        gsl_matrix_set(data4, 0, 0, 1.0);
+        gsl_matrix_set(data4, 0, 1, 2.0);
+        gsl_matrix_set(data4, 0, 2, 3.0);
+        gsl_matrix_set(data4, 0, 3, 4.0);
+        gsl_test(gsl_extrap_richardson(data4, 2, k2, normtable, NULL),
+                "Unexpected error reported in %s");
+        gsl_test_abs(gsl_matrix_get(data4, 0, 0), zx, GSL_DBL_EPSILON,
+                "%s scalar correct result at %s:%d",
+                __func__, __FILE__, __LINE__);
+
+        for (i = 0; i < normtable->size1 - 1; ++i) {
+            for (j = i+1; j < normtable->size2; ++j) {
+                gsl_test(!gsl_isnan(gsl_matrix_get(normtable, i, j)),
+                        "%s expected NAN in normtable at (%d, %d)",
+                        __func__, i, j);
+            }
+        }
+
+        gsl_test_abs(gsl_matrix_get(normtable, 0, 0), 1.0, GSL_DBL_EPSILON,
+                "%s normtable (0,0) value", __func__);
+        gsl_test_abs(gsl_matrix_get(normtable, 1, 0), 2.0, GSL_DBL_EPSILON,
+                "%s normtable (1,0) value", __func__);
+        gsl_test_abs(gsl_matrix_get(normtable, 2, 0), 3.0, GSL_DBL_EPSILON,
+                "%s normtable (2,0) value", __func__);
+        gsl_test_abs(gsl_matrix_get(normtable, 3, 0), 4.0, GSL_DBL_EPSILON,
+                "%s normtable (3,0) value", __func__);
+        gsl_test_abs(gsl_matrix_get(normtable, 1, 1), xx, GSL_DBL_EPSILON,
+                "%s normtable (1,1) value", __func__);
+        gsl_test_abs(gsl_matrix_get(normtable, 2, 1), xy, GSL_DBL_EPSILON,
+                "%s normtable (2,1) value", __func__);
+        gsl_test_abs(gsl_matrix_get(normtable, 3, 1), xz, GSL_DBL_EPSILON,
+                "%s normtable (3,1) value", __func__);
+        gsl_test_abs(gsl_matrix_get(normtable, 2, 2), yx, GSL_DBL_EPSILON,
+                "%s normtable (2,2) value", __func__);
+        gsl_test_abs(gsl_matrix_get(normtable, 3, 2), yy, GSL_DBL_EPSILON,
+                "%s normtable (3,2) value", __func__);
+        gsl_test_abs(gsl_matrix_get(normtable, 3, 3), zx, GSL_DBL_EPSILON,
+                "%s normtable (3,3) value", __func__);
+
+        gsl_matrix_free(normtable);
+    }
+
+    gsl_vector_free(k3);
+    gsl_vector_free(k2);
+    gsl_vector_free(k1);
+    gsl_matrix_free(data4);
+    gsl_matrix_free(data2);
+}
+
+void
+test_richardson_extrapolation_vectorinputs()
+{
+    {
+        gsl_matrix * data = gsl_matrix_alloc(2,2);
+        gsl_matrix_set(data, 0, 0, 1.0);
+        gsl_matrix_set(data, 0, 1, 2.0);
+        gsl_matrix_set(data, 1, 0, 2.0);
+        gsl_matrix_set(data, 1, 1, 3.0);
+
+        gsl_matrix * normtable = gsl_matrix_alloc(data->size2, data->size2);
+
+        gsl_vector * exact = gsl_vector_alloc(data->size2);
+        gsl_vector_set(exact, 0, 3.0);
+        gsl_vector_set(exact, 1, 4.0);
+
+        gsl_test(gsl_extrap_richardson(
+                    data, 2, NULL, normtable, exact),
+                "Unexpected error reported in %s");
+        gsl_test_abs(gsl_matrix_get(data, 0, 0), 3.0, GSL_DBL_EPSILON,
+                "%s data (0,0) value at %s:%d", __func__, __FILE__, __LINE__);
+        gsl_test_abs(gsl_matrix_get(data, 1, 0), 4.0, GSL_DBL_EPSILON,
+                "%s data (1,0) value at %s:%d", __func__, __FILE__, __LINE__);
+
+        gsl_test_abs(gsl_matrix_get(normtable, 1, 1), 0.0, GSL_DBL_EPSILON,
+                "%s normtable (1,1) value at %s:%d",
+                __func__, __FILE__, __LINE__);
+
+        gsl_vector_free(exact);
+        gsl_matrix_free(normtable);
+        gsl_matrix_free(data);
+    }
+
+    {
+        gsl_matrix * data = gsl_matrix_alloc(2,3);
+        gsl_matrix_set(data, 0, 0, 1.0);
+        gsl_matrix_set(data, 0, 1, 2.0);
+        gsl_matrix_set(data, 0, 2, 3.0);
+        gsl_matrix_set(data, 1, 0, 2.0);
+        gsl_matrix_set(data, 1, 1, 3.0);
+        gsl_matrix_set(data, 1, 2, 4.0);
+
+        gsl_test(gsl_extrap_richardson(data, 2, NULL, NULL, NULL),
+                "Unexpected error reported in %s");
+        gsl_test_abs(gsl_matrix_get(data, 0, 0), 13.0/3.0, GSL_DBL_EPSILON,
+                "%s data (0,0) value at %s:%d", __func__, __FILE__, __LINE__);
+        gsl_test_abs(gsl_matrix_get(data, 1, 0), 16.0/3.0, GSL_DBL_EPSILON,
+                "%s data (1,0) value at %s:%d", __func__, __FILE__, __LINE__);
+
+        gsl_matrix_free(data);
+    }
+}
+
+int
+main(int argc, char **argv)
+{
+    gsl_ieee_env_setup();
+
+    test_richardson_extrapolation_step();
+    test_richardson_extrapolation_defaults();
+    test_richardson_extrapolation_twolevels();
+    test_richardson_extrapolation_multiplelevels();
+    test_richardson_extrapolation_vectorinputs();
+
+    exit(gsl_test_summary());
+}
-- 
1.6.3.3

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

[Brian Gough]

At Sun, 13 Dec 2009 17:49:28 -0600,
Rhys Ulerich wrote:
> Attached is a patch adding Richardson extrapolation routines, tests,
> documentation, and an example.  Please let me know if I can answer any
> questions.

Thanks, I have tested it out.  One question: it looks like this
extrapolates a vector quantity.  For simplicity, would it make sense
to work with scalars as in the gsl_sum routines -- since presumably
that is the common case -- or is there some application where vector
extrapolation is unavoidable? (could each component can be
extrapolated independently?)

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

[Rhys Ulerich]

Thanks for looking through the patch Brian.

> One question: it looks like this
> extrapolates a vector quantity.  For simplicity, would it make sense
> to work with scalars as in the gsl_sum routines -- since presumably
> that is the common case -- or is there some application where vector
> extrapolation is unavoidable? (could each component can be
> extrapolated independently?)

You're right that each component could be extrapolated independently.
I wrote it using vectors because doing so allows using BLAS calls for
the linear algebra and speeds up many component use cases.  I'd prefer
to keep the code vector-capable under the covers under the theory that
people extrapolating only a scalar at a time aren't all that worried
about speed.

My vector use case arose from testing convergence orders for
timesteppers (non-GSL, but similar in interface) that manipulate
vectors of state variables.

If you'd like, I can submit an additional patch that 1) Renames
gsl_extrap_richardson to something like gsl_extrap_richardson_vector,
2) Provides  a scalar-only gsl_extrap_richardson that wraps up scalars
and calls the gsl_extrap_richardson_vector routines under the covers,
and 3) Exercises the scalar-only wrappers fully in the test suite.
That way if someone later does need a faster scalar-only version, he
or she can change over the underlying implementation without fear of
regression.

- Rhys

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

[Rhys Ulerich]

[-- Attachment #1: Type: text/plain, Size: 638 bytes --]

> ... I can submit an additional patch that 1) Renames
> gsl_extrap_richardson to something like gsl_extrap_richardson_vector,
> 2) Provides  a scalar-only gsl_extrap_richardson that wraps up scalars
> and calls the gsl_extrap_richardson_vector routines under the covers,
> and 3) Exercises the scalar-only wrappers fully in the test suite.
> That way if someone later does need a faster scalar-only version, he
> or she can change over the underlying implementation without fear of
> regression.

Attached is a second patch that does exactly what I described above.  It builds
atop the prior patch from this thread.

- Rhys

[-- Attachment #2: 0001-Created-scalar-version-of-Richardson-extrapolation.patch --]
[-- Type: text/x-patch, Size: 68285 bytes --]

From 0e6d3460d5cac6777a940b07a5f29288f2addab1 Mon Sep 17 00:00:00 2001
From: Rhys Ulerich <rhys.ulerich@gmail.com>
Date: Tue, 15 Dec 2009 22:21:28 -0600
Subject: [PATCH] Created scalar-version of Richardson extrapolation

Renames gsl_extrap_richardson to gsl_extrap_richardson_vector,
Provides  a scalar-only gsl_extrap_richardson that wraps scalars
   and calls gsl_extrap_richardson_vector under the covers,
Exercises the scalar-only wrappers fully in the test suite,
Documents the new scalar-only methods, and
Updates the example code.
---
 doc/examples/richardson.c   |   41 +-
 doc/examples/richardson.out |    4 +-
 doc/sum.texi                |   72 ++-
 extrap/ChangeLog            |    4 +
 extrap/demo_richardson.c    |   41 +-
 extrap/gsl_extrap.h         |   28 +-
 extrap/richardson.c         |   64 ++-
 extrap/test_richardson.c    | 1304 ++++++++++++++++++++++++++++---------------
 8 files changed, 1024 insertions(+), 534 deletions(-)

diff --git a/doc/examples/richardson.c b/doc/examples/richardson.c
index caa76bc..a942999 100644
--- a/doc/examples/richardson.c
+++ b/doc/examples/richardson.c
@@ -19,42 +19,34 @@ main (void)
 {
   printf("\nPerforming a single Richardson extrapolation step...\n");
   {
-    gsl_vector *Ah  = gsl_vector_alloc(1);
-    gsl_vector *Aht = gsl_vector_alloc(1);
-
-    gsl_vector_set(Ah,  0, estimate(h));
-    gsl_vector_set(Aht, 0, estimate(h / t));
+    double Ah        = estimate(h);
+    const double Aht = estimate(h / t);
 
     printf("Exact value for (d/dx) sin(%4f) =  %14.12f\n", x, cos(x));
     printf("Approximate value for h   = %4f is %14.12f with error %14e\n",
-           h, gsl_vector_get(Ah, 0), fabs(cos(x) - gsl_vector_get(Ah, 0)));
+           h, Ah, fabs(cos(x) - Ah));
     printf("Approximate value for h/t = %4f is %14.12f with error %14e\n",
-           h / t, gsl_vector_get(Aht, 0), fabs(cos(x) - gsl_vector_get(Aht, 0)));
+           h / t, Aht, fabs(cos(x) - Aht));
 
-    gsl_extrap_richardson_step(Ah, Aht, t, 2.0);
+    gsl_extrap_richardson_step(&Ah, &Aht, t, 2.0);
 
     printf("After one extrapolation step value is   %14.12f with error %14e\n",
-           gsl_vector_get(Ah, 0), fabs(cos(x) - gsl_vector_get(Ah, 0)));
-
-    gsl_vector_free(Aht);
-    gsl_vector_free(Ah);
+           Ah, fabs(cos(x) - Ah));
   }
 
   printf("\nPerforming iterated Richardson extrapolation...\n");
   {
-    gsl_matrix *A         = gsl_matrix_alloc(1, 4);
+    gsl_vector *A         = gsl_vector_alloc(4);
     gsl_vector *k         = gsl_vector_alloc(2);
-    gsl_matrix *normtable = gsl_matrix_alloc(A->size2, A->size2);
-    gsl_vector *exact     = gsl_vector_alloc(1);
+    gsl_matrix *normtable = gsl_matrix_alloc(A->size, A->size);
+    const double exact    = cos(x);
     int i, j;
 
-    gsl_vector_set(exact, 0, cos(x));
-
     printf("Computing %d initial estimates of (d/dx) sin(%f)\n",
-           (int) A->size2, x);
-    for (i = 0; i < A->size2; ++i)
+           (int) A->size, x);
+    for (i = 0; i < A->size; ++i)
       {
-        gsl_matrix_set(A, 0, i, estimate(h / pow(t, i)));
+        gsl_vector_set(A, i, estimate(h / pow(t, i)));
       }
 
     printf("Using leading error terms like O(h^2), O(h^4), O(h^6), O(h^8)...\n");
@@ -64,20 +56,19 @@ main (void)
     gsl_extrap_richardson(A, t, k, normtable, exact);
 
     printf("The l_2 error after each iterated extrapolation step is:\n");
-    for (i = 0; i < A->size2; ++i)
+    for (i = 0; i < A->size; ++i)
       {
-        for (j = 0; j < A->size2; ++j)
+        for (j = 0; j < A->size; ++j)
           {
             printf(" %16e", gsl_matrix_get(normtable, i, j));
           }
         printf("\n");
       }
-    printf("\nThe final estimate is %14.12f\n", gsl_matrix_get(A, 0, 0));
+    printf("\nThe final estimate is %14.12f\n", gsl_vector_get(A, 0));
 
-    gsl_vector_free(exact);
     gsl_matrix_free(normtable);
     gsl_vector_free(k);
-    gsl_matrix_free(A);
+    gsl_vector_free(A);
   }
 
   return 0;
diff --git a/doc/examples/richardson.out b/doc/examples/richardson.out
index b055ff2..a54ae94 100644
--- a/doc/examples/richardson.out
+++ b/doc/examples/richardson.out
@@ -11,7 +11,7 @@ Using leading error terms like O(h^2), O(h^4), O(h^6), O(h^8)...
 The l_2 error after each iterated extrapolation step is:
      1.935620e-02              nan              nan              nan
      4.868179e-03     3.883801e-05              nan              nan
-     1.218872e-03     2.436061e-06     9.264229e-09              nan
-     3.048323e-04     1.523898e-07     1.450697e-10     3.211875e-13
+     1.218872e-03     2.436061e-06     9.264230e-09              nan
+     3.048323e-04     1.523898e-07     1.450698e-10     3.212985e-13
 
 The final estimate is 0.731688868873
diff --git a/doc/sum.texi b/doc/sum.texi
index 4c920b3..c8af385 100644
--- a/doc/sum.texi
+++ b/doc/sum.texi
@@ -187,17 +187,19 @@ In a Richardson extrapolation step two estimates @math{A_i(h)} and
 @math{A_i(h/t)} are combined to find an estimate of order @math{h^{k_(i+1)}}
 according to @math{A_(i+1)(h) = A_i(h/t) + (t^n - 1)^(-1)*(A_i(h/t) - A_i(h))}.
 Note that this did not require knowing the leading coefficient @math{a_i}, only
-the leading error order @math{k_i}.  One such step may be computed using
-@code{gsl_extrap_richardson_step}.  The process can be iterated to combine
-@math{A_(i+1)(h)} and @math{A_(i+1)(h/t)} to obtain @math{A_(i+2)(h)}.  The
-iterated process, as well as convenience mechanism to compute error norms given
-a known solution, can be performed using using @code{gsl_extrap_richardson}.
-
-@deftypefun {int} gsl_extrap_richardson_step ( gsl_vector * @var{Ah}, const gsl_vector * @var{Aht}, const double @var{t}, const double @var{ki})
+the leading error order @math{k_i}.  One such step may be computed for a scalar
+using @code{gsl_extrap_richardson_step} or for a vector using
+@code{gsl_extrap_richardson_vector_step}.  The process can be iterated to
+combine @math{A_(i+1)(h)} and @math{A_(i+1)(h/t)} to obtain @math{A_(i+2)(h)}.
+The iterated process, as well as convenience mechanism to compute error norms
+given a known solution, can be performed using using either
+@code{gsl_extrap_richardson} or @code{gsl_extrap_richardson_vector}.
+
+@deftypefun {int} gsl_extrap_richardson_step ( double * @var{Ah}, const double * @var{Aht}, double @var{t}, double @var{ki})
 Given @math{A_i(h)}, an approximation of @math{A} such that @math{A - A(h) =
-a_i*h^k_i + a_(i+1)*h^k_(i+1) + ...} use Richardson extrapolation to
-find an approximation to leading order @math{k_(i+1)} from evaluations at
-@math{A_i(h)} and @math{A_i(h/t)} for @math{t > 0}.
+a_i*h^k_i + a_(i+1)*h^k_(i+1) + ...} use Richardson extrapolation to find an
+approximation to leading order @math{k_(i+1)} from evaluations at @math{A_i(h)}
+and @math{A_i(h/t)} for @math{t > 0}.
 
 On entry, @var{Ah} contains the coarser approximation @math{A_i(h)} and
 @var{Aht} contains the finer approximation @math{A_i(h/t)}.  @var{ki} is
@@ -210,34 +212,50 @@ the method calls @code{gsl_error} and returns on of GSL's error constants.  The
 parameter @var{Ah} will be in an undefined state after any error.
 @end deftypefun
 
-@deftypefun {int} gsl_extrap_richardson ( gsl_matrix * const @var{A}, const double @var{t}, const gsl_vector * @var{k}, gsl_matrix * @var{normtable}, const gsl_vector * const @var{exact})
-Perform Richardson extrapolation on a sequence of refined approximations.
-The routine can perform multiple step extrapolation, e.g. using
-@math{A_0(h)}, @math{A_0(h/2)}, and @math{A_0(h/4)} to compute @math{A_2(h)}.
-See @code{gsl_extrap_richardson_step} for the terminology in use.
+@deftypefun {int} gsl_extrap_richardson_vector_step ( gsl_vector * @var{Ah}, const gsl_vector * @var{Aht}, double @var{t}, double @var{ki})
+This method performs a fast, vector-valued Richardson extrapolation process
+equivalent to repeatedly calling @code{gsl_extrap_richardson_step} on each
+index in the vectors @var{Ah} and @var{Aht}.  All parameter processing and
+error handling follows @code{gsl_extrap_richardson_step}.
+@end deftypefun
+
+@deftypefun {int} gsl_extrap_richardson ( gsl_vector * @var{A}, double @var{t}, const gsl_vector * @var{k}, gsl_matrix * @var{normtable}, double @var{exact})
+Perform Richardson extrapolation on a sequence of refined approximations.  The
+routine can perform multiple step extrapolation, e.g. using @math{A_0(h)},
+@math{A_0(h/2)}, and @math{A_0(h/4)} to compute @math{A_2(h)}.  See
+@code{gsl_extrap_richardson_step} for the terminology in use.
 
-On entry, column @math{i} of @var{A} contains @math{A_0(h/t^i)} while @var{t}
+On entry, component @math{i} of @var{A} contains @math{A_0(h/t^i)} while @var{t}
 contains the refinement factor between each pair of approximations.  The
 refinement factor must be fixed for all columns.  The leading error term orders
 are provided in @var{k}.  They are assumed to start at 1 if @var{k} is
 @code{NULL}.  They are assumed to increment by 1 if no second element is
-provided.  If the vector @var{k} is shorter than @code{A->size2}, any
+provided.  If the vector @var{k} is shorter than @code{A->size}, any
 unspecified higher index entries are assumed to grow by
 @code{gsl_vector_get(k,k->size-1) - gsl_vector_get(k,k->size-2)}.
 
 If @var{normtable} is non-@code{NULL} on entry, the routine produces in
 @var{normtable} a table showing the @math{l_2} error at each step in the
 extrapolation process calculated against the provided exact solution in
-@var{exact}.  If @var{exact} is @code{NULL} then it is treated as as the zero
-vector and the resulting @var{normtable} simply contains the @math{l_2} norm of
-each extrapolation step.
-
-On success, the method returns @code{GSL_SUCCESS}.  On successful exit, column
-0 of @var{A} contains @math{A_(A->size2)(h)} and all other columns will have
-been overwritten with temporary results.  If provided, @var{normtable} will be
-modified as described above.  On error the method calls @code{gsl_error} and
-returns one of GSL's error codes.  The parameters @var{A} and @var{normtable}
-will be in an undefined state after any error.
+@var{exact}.  If @var{exact} is zero then the the resulting @var{normtable}
+simply contains the @math{l_2} norm of each extrapolation step.
+
+On success, the method returns @code{GSL_SUCCESS}.  On successful exit,
+component 0 of @var{A} contains @math{A_(A->size)(h)} and all other entries
+will have been overwritten with temporary results.  If provided,
+@var{normtable} will be modified as described above.  On error the method calls
+@code{gsl_error} and returns one of GSL's error codes.  The parameters @var{A}
+and @var{normtable} will be in an undefined state after any error.
+@end deftypefun
+
+@deftypefun {int} gsl_extrap_richardson_vector ( gsl_matrix * @var{A}, double @var{t}, const gsl_vector * @var{k}, gsl_matrix * @var{normtable}, const gsl_vector * @var{exact})
+This method performs a fast, vector-valued Richardson extrapolation process
+equivalent to repeatedly calling @code{gsl_extrap_richardson} on each row in
+@var{A}.  On entry, each column @math{i} of @var{A} represents the
+vector-valued @math{A_0(h/t^i)}.  If @var{exact} is @code{NULL} then it is
+treated as as the zero vector and any resulting @var{normtable} simply contains
+the @math{l_2} norm of each extrapolation step.  All parameter processing and
+error handling follows @code{gsl_extrap_richardson_vector}.
 @end deftypefun
 
 @node Example of Richardson extrapolation
diff --git a/extrap/ChangeLog b/extrap/ChangeLog
index 9183ab4..2a07a56 100644
--- a/extrap/ChangeLog
+++ b/extrap/ChangeLog
@@ -1,3 +1,7 @@
+2009-12-15  Rhys Ulerich  <rhys.ulerich@gmail.com>
+
+	* added scalar interface and tests by wrapping vector routines
+
 2009-12-13  Rhys Ulerich  <rhys.ulerich@gmail.com>
 
 	* added build structure and Richardson extrapolation
diff --git a/extrap/demo_richardson.c b/extrap/demo_richardson.c
index caa76bc..a942999 100644
--- a/extrap/demo_richardson.c
+++ b/extrap/demo_richardson.c
@@ -19,42 +19,34 @@ main (void)
 {
   printf("\nPerforming a single Richardson extrapolation step...\n");
   {
-    gsl_vector *Ah  = gsl_vector_alloc(1);
-    gsl_vector *Aht = gsl_vector_alloc(1);
-
-    gsl_vector_set(Ah,  0, estimate(h));
-    gsl_vector_set(Aht, 0, estimate(h / t));
+    double Ah        = estimate(h);
+    const double Aht = estimate(h / t);
 
     printf("Exact value for (d/dx) sin(%4f) =  %14.12f\n", x, cos(x));
     printf("Approximate value for h   = %4f is %14.12f with error %14e\n",
-           h, gsl_vector_get(Ah, 0), fabs(cos(x) - gsl_vector_get(Ah, 0)));
+           h, Ah, fabs(cos(x) - Ah));
     printf("Approximate value for h/t = %4f is %14.12f with error %14e\n",
-           h / t, gsl_vector_get(Aht, 0), fabs(cos(x) - gsl_vector_get(Aht, 0)));
+           h / t, Aht, fabs(cos(x) - Aht));
 
-    gsl_extrap_richardson_step(Ah, Aht, t, 2.0);
+    gsl_extrap_richardson_step(&Ah, &Aht, t, 2.0);
 
     printf("After one extrapolation step value is   %14.12f with error %14e\n",
-           gsl_vector_get(Ah, 0), fabs(cos(x) - gsl_vector_get(Ah, 0)));
-
-    gsl_vector_free(Aht);
-    gsl_vector_free(Ah);
+           Ah, fabs(cos(x) - Ah));
   }
 
   printf("\nPerforming iterated Richardson extrapolation...\n");
   {
-    gsl_matrix *A         = gsl_matrix_alloc(1, 4);
+    gsl_vector *A         = gsl_vector_alloc(4);
     gsl_vector *k         = gsl_vector_alloc(2);
-    gsl_matrix *normtable = gsl_matrix_alloc(A->size2, A->size2);
-    gsl_vector *exact     = gsl_vector_alloc(1);
+    gsl_matrix *normtable = gsl_matrix_alloc(A->size, A->size);
+    const double exact    = cos(x);
     int i, j;
 
-    gsl_vector_set(exact, 0, cos(x));
-
     printf("Computing %d initial estimates of (d/dx) sin(%f)\n",
-           (int) A->size2, x);
-    for (i = 0; i < A->size2; ++i)
+           (int) A->size, x);
+    for (i = 0; i < A->size; ++i)
       {
-        gsl_matrix_set(A, 0, i, estimate(h / pow(t, i)));
+        gsl_vector_set(A, i, estimate(h / pow(t, i)));
       }
 
     printf("Using leading error terms like O(h^2), O(h^4), O(h^6), O(h^8)...\n");
@@ -64,20 +56,19 @@ main (void)
     gsl_extrap_richardson(A, t, k, normtable, exact);
 
     printf("The l_2 error after each iterated extrapolation step is:\n");
-    for (i = 0; i < A->size2; ++i)
+    for (i = 0; i < A->size; ++i)
       {
-        for (j = 0; j < A->size2; ++j)
+        for (j = 0; j < A->size; ++j)
           {
             printf(" %16e", gsl_matrix_get(normtable, i, j));
           }
         printf("\n");
       }
-    printf("\nThe final estimate is %14.12f\n", gsl_matrix_get(A, 0, 0));
+    printf("\nThe final estimate is %14.12f\n", gsl_vector_get(A, 0));
 
-    gsl_vector_free(exact);
     gsl_matrix_free(normtable);
     gsl_vector_free(k);
-    gsl_matrix_free(A);
+    gsl_vector_free(A);
   }
 
   return 0;
diff --git a/extrap/gsl_extrap.h b/extrap/gsl_extrap.h
index b15e97e..0ba5d01 100644
--- a/extrap/gsl_extrap.h
+++ b/extrap/gsl_extrap.h
@@ -36,20 +36,34 @@
 __BEGIN_DECLS
 
 int
-gsl_extrap_richardson_step(
+gsl_extrap_richardson_vector_step(
         gsl_vector * Ah,
         const gsl_vector * Aht,
-        const double t,
-        const double ki);
+        double t,
+        double ki);
 
 int
-gsl_extrap_richardson(
-        gsl_matrix * const A,
-        const double t,
+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 * const exact);
+        const gsl_vector * exact);
 
+int
+gsl_extrap_richardson(
+        gsl_vector * A,
+        double t,
+        const gsl_vector * k,
+        gsl_matrix * normtable,
+        double exact);
 
 __END_DECLS
 
diff --git a/extrap/richardson.c b/extrap/richardson.c
index 4b207cb..602f276 100644
--- a/extrap/richardson.c
+++ b/extrap/richardson.c
@@ -26,9 +26,9 @@
 #include <gsl/gsl_vector_double.h>
 
 int
-gsl_extrap_richardson_step(
-  gsl_vector * Ah,
-  const gsl_vector * Aht,
+gsl_extrap_richardson_vector_step(
+  gsl_vector * const Ah,
+  const gsl_vector * const Aht,
   const double t,
   const double ki)
 {
@@ -37,6 +37,11 @@ gsl_extrap_richardson_step(
       GSL_ERROR("ki == 0 invalid as it will cause a division-by-zero",
                 GSL_EDOM);
     }
+  if (t <= 0)
+    {
+      GSL_ERROR("t <= 0 invalid as represents a nonsensical refinement",
+                GSL_EDOM);
+    }
 
   const double tki       = pow(t, ki);
   const double inv_tkim1 = 1.0 / (tki - 1.0);
@@ -50,11 +55,25 @@ gsl_extrap_richardson_step(
 }
 
 int
-gsl_extrap_richardson(
+gsl_extrap_richardson_step(
+  double * const Ah,
+  const double * const Aht,
+  const double t,
+  const double ki)
+{
+  gsl_vector_view       Ah_view  = gsl_vector_view_array(Ah, 1);
+  gsl_vector_const_view Aht_view = gsl_vector_const_view_array(Aht, 1);
+
+  return gsl_extrap_richardson_vector_step(
+           &Ah_view.vector, &Aht_view.vector, t, ki);
+}
+
+int
+gsl_extrap_richardson_vector(
   gsl_matrix * const A,
   const double t,
-  const gsl_vector * k,
-  gsl_matrix * normtable,
+  const gsl_vector * const k,
+  gsl_matrix * const normtable,
   const gsl_vector * const exact)
 {
   gsl_vector * scratch = NULL;
@@ -146,7 +165,7 @@ gsl_extrap_richardson(
           gsl_vector_view Aih  = gsl_matrix_column(A, j);
           gsl_vector_view Aiht = gsl_matrix_column(A, j + 1);
 
-          const int error = gsl_extrap_richardson_step(
+          const int error = gsl_extrap_richardson_vector_step(
                               &Aih.vector, &Aiht.vector, t, ki);
           if (error)
             {
@@ -182,3 +201,34 @@ gsl_extrap_richardson(
 
   return GSL_SUCCESS;
 }
+
+int
+gsl_extrap_richardson(
+  gsl_vector * const A,
+  const double t,
+  const gsl_vector * const k,
+  gsl_matrix * const normtable,
+  const double exact)
+{
+  gsl_matrix_view       A_view     = gsl_matrix_view_vector(A, 1, A->size);
+  gsl_vector_const_view exact_view = gsl_vector_const_view_array(&exact, 1);
+
+  if (A->stride != 1)
+    {
+      GSL_ERROR("Unable to create required matrix view for A->stride != 1",
+                GSL_EINVAL);
+    }
+
+  /* Provide exact value iff normtable is non-NULL, otherwise we
+     encounter a consistency check in gsl_extrap_richardson_vector */
+  if (normtable)
+    {
+      return gsl_extrap_richardson_vector(
+               &A_view.matrix, t, k, normtable, &exact_view.vector);
+    }
+  else
+    {
+      return gsl_extrap_richardson_vector(
+               &A_view.matrix, t, k, NULL, NULL);
+    }
+}
diff --git a/extrap/test_richardson.c b/extrap/test_richardson.c
index 903ccd7..cfc1ccc 100644
--- a/extrap/test_richardson.c
+++ b/extrap/test_richardson.c
@@ -5,470 +5,892 @@
 #include <gsl/gsl_test.h>
 
 void
-test_richardson_extrapolation_step()
+test_richardson_vector_extrapolation_step()
 {
-    const int    ki           = 1;
-    const size_t n            = 1;
-    gsl_vector *Ah   = gsl_vector_alloc(n);
-    gsl_vector *Aht  = gsl_vector_alloc(n);
-
-    {
-        const double t = 2.0;
-        gsl_vector_set(Ah,  0, 1.0);
-        gsl_vector_set(Aht, 0, 2.0);
-        gsl_test(gsl_extrap_richardson_step(Ah, Aht, t, ki),
-                "Unexpected error reported in %s for t=%f", __func__, t);
-        gsl_test_abs(gsl_vector_get(Ah, 0), 3.0, GSL_DBL_EPSILON,
-                "Extrapolation in %s for ki=%d, t=%f", __func__, ki, t);
-    }
-
-    {
-        const double t = 3.0;
-        gsl_vector_set(Ah,  0, 1.0);
-        gsl_vector_set(Aht, 0, 2.0);
-        gsl_test(gsl_extrap_richardson_step(Ah, Aht, t, ki),
-                "Unexpected error reported in %s for t=%f", __func__, t);
-        gsl_test_abs(gsl_vector_get(Ah, 0), 5.0/2.0, GSL_DBL_EPSILON,
-                "Extrapolation in %s for ki=%d, t=%f", __func__, ki, t);
-    }
-
-    gsl_vector_free(Aht);
-    gsl_vector_free(Ah);
+  const int    ki           = 1;
+  const size_t n            = 1;
+  gsl_vector *Ah   = gsl_vector_alloc(n);
+  gsl_vector *Aht  = gsl_vector_alloc(n);
+
+  {
+    const double t = 2.0;
+    gsl_vector_set(Ah,  0, 1.0);
+    gsl_vector_set(Aht, 0, 2.0);
+    gsl_test(gsl_extrap_richardson_vector_step(Ah, Aht, t, ki),
+             "Unexpected error reported in %s for t=%f", __func__, t);
+    gsl_test_abs(gsl_vector_get(Ah, 0), 3.0, GSL_DBL_EPSILON,
+                 "Extrapolation in %s for ki=%d, t=%f", __func__, ki, t);
+  }
+
+  {
+    const double t = 3.0;
+    gsl_vector_set(Ah,  0, 1.0);
+    gsl_vector_set(Aht, 0, 2.0);
+    gsl_test(gsl_extrap_richardson_vector_step(Ah, Aht, t, ki),
+             "Unexpected error reported in %s for t=%f", __func__, t);
+    gsl_test_abs(gsl_vector_get(Ah, 0), 5.0 / 2.0, GSL_DBL_EPSILON,
+                 "Extrapolation in %s for ki=%d, t=%f", __func__, ki, t);
+  }
+
+  gsl_vector_free(Aht);
+  gsl_vector_free(Ah);
 }
 
 void
-test_richardson_extrapolation_defaults()
+test_richardson_vector_extrapolation_defaults()
 {
-    gsl_matrix * data1 = gsl_matrix_alloc(1, 2);
-    gsl_matrix * data2 = gsl_matrix_alloc(data1->size1, data1->size2);
-    {
-        gsl_matrix_set(data1, 0, 0, 1.0);
-        gsl_matrix_set(data1, 0, 1, 2.0);
-
-        gsl_test(gsl_extrap_richardson(data1, 2, NULL, NULL, NULL),
-                "Unexpected error reported in %s");
-        gsl_test_abs(gsl_matrix_get(data1, 0, 0), 3.0, GSL_DBL_EPSILON,
-                "%s scalar correct result", __func__);
-    }
-
-    {
-        gsl_matrix_set(data1, 0, 0, 1.0);
-        gsl_matrix_set(data1, 0, 1, 2.0);
-
-        gsl_test(gsl_extrap_richardson(data1, 3, NULL, NULL, NULL),
-                "Unexpected error reported in %s");
-        gsl_test_abs(gsl_matrix_get(data1, 0, 0), 5.0/2.0, GSL_DBL_EPSILON,
-                "%s scalar correct result", __func__);
-    }
-
-    {
-        gsl_matrix_set(data1, 0, 0, 1.0);
-        gsl_matrix_set(data1, 0, 1, 2.0);
-        gsl_matrix_memcpy(data2, data1);
-
-        gsl_test(gsl_extrap_richardson(data1, 2, NULL, NULL, NULL),
-                "Unexpected error reported in %s");
-        const double result1 = gsl_matrix_get(data1, 0, 0);
-
-        gsl_test_abs(result1, 3.0, GSL_DBL_EPSILON,
-                "%s scalar correct result", __func__);
-
-        gsl_vector * k = gsl_vector_alloc(1);
-        gsl_vector_set(k, 0, 1);
-        gsl_test(gsl_extrap_richardson(data2, 2, k, NULL, NULL),
-                "Unexpected error reported in %s");
-        gsl_vector_free(k);
-        const double result2 = gsl_matrix_get(data2, 0, 0);
-
-        gsl_test_abs(result1, result2, GSL_DBL_EPSILON,
-                "%s with k == NULL", __func__);
-    }
-
-    {
-        int i, j;
-        gsl_matrix_set(data1, 0, 0, 1.0);
-        gsl_matrix_set(data1, 0, 1, 2.0);
-        gsl_matrix_memcpy(data2, data1);
-
-        gsl_matrix * normtable1 = gsl_matrix_alloc(data1->size2, data1->size2);
-        gsl_matrix * normtable2 = gsl_matrix_alloc(data1->size2, data1->size2);
-
-        gsl_test(gsl_extrap_richardson(
-                    data1, 2, NULL, normtable1, NULL),
-                "Unexpected error reported in %s");
-        const double result1 = gsl_matrix_get(data1, 0, 0);
-
-        gsl_test_abs(result1, 3.0, GSL_DBL_EPSILON,
-                "%s correct extrapolation answer");
-
-        for (i = 0; i < normtable1->size1 - 1; ++i) {
-            for (j = i+1; j < normtable1->size2; ++j) {
-                gsl_test(!gsl_isnan(gsl_matrix_get(normtable1, i, j)),
-                        "%s expected NAN in normtable1 at (%d, %d)",
-                        __func__, i, j);
-            }
-        }
-
-        gsl_test_abs(gsl_matrix_get(normtable1, 0, 0), 1.0, GSL_DBL_EPSILON,
-                "%s normtable1 (0,0) value", __func__);
-        gsl_test_abs(gsl_matrix_get(normtable1, 1, 0), 2.0, GSL_DBL_EPSILON,
-                "%s normtable1 (1,0) value", __func__);
-        gsl_test_abs(gsl_matrix_get(normtable1, 1, 1), 3.0, GSL_DBL_EPSILON,
-                "%s normtable1 (1,1) value", __func__);
-
-        gsl_vector * exact = gsl_vector_alloc(1);
-        gsl_vector_set_zero(exact);
-        gsl_test(gsl_extrap_richardson(
-                    data2, 2.0, NULL, normtable2, exact),
-                "Unexpected error reported in %s");
-        gsl_vector_free(exact);
-        const double result2 = gsl_matrix_get(data2, 0, 0);
-
-        for (i = 0; i < normtable2->size1 - 1; ++i) {
-            for (j = i+1; j < normtable2->size2; ++j) {
-                gsl_test(!gsl_isnan(gsl_matrix_get(normtable2, i, j)),
-                        "%s expected NAN in normtable2 at (%d, %d)",
-                        __func__, i, j);
-            }
-        }
-
-        gsl_test_abs(result1, result2, GSL_DBL_EPSILON,
-                "%s with exact == NULL", __func__);
-
-        for (i = 0; i < normtable1->size1; ++i) {
-            for (j = 0; j < i + 1; ++j) {
-                gsl_test_abs(gsl_matrix_get(normtable1, i, j),
-                             gsl_matrix_get(normtable2, i, j),
-                             GSL_DBL_EPSILON,
-                             "%s identical normtable results at (%d, %d)",
-                             __func__, i, j);
-            }
-        }
-
-        gsl_matrix_free(normtable2);
-        gsl_matrix_free(normtable1);
-    }
-
-    gsl_matrix_free(data2);
-    gsl_matrix_free(data1);
+  gsl_matrix * data1 = gsl_matrix_alloc(1, 2);
+  gsl_matrix * data2 = gsl_matrix_alloc(data1->size1, data1->size2);
+  {
+    gsl_matrix_set(data1, 0, 0, 1.0);
+    gsl_matrix_set(data1, 0, 1, 2.0);
+
+    gsl_test(gsl_extrap_richardson_vector(data1, 2, NULL, NULL, NULL),
+             "Unexpected error reported in %s");
+    gsl_test_abs(gsl_matrix_get(data1, 0, 0), 3.0, GSL_DBL_EPSILON,
+                 "%s scalar correct result", __func__);
+  }
+
+  {
+    gsl_matrix_set(data1, 0, 0, 1.0);
+    gsl_matrix_set(data1, 0, 1, 2.0);
+
+    gsl_test(gsl_extrap_richardson_vector(data1, 3, NULL, NULL, NULL),
+             "Unexpected error reported in %s");
+    gsl_test_abs(gsl_matrix_get(data1, 0, 0), 5.0 / 2.0, GSL_DBL_EPSILON,
+                 "%s scalar correct result", __func__);
+  }
+
+  {
+    gsl_matrix_set(data1, 0, 0, 1.0);
+    gsl_matrix_set(data1, 0, 1, 2.0);
+    gsl_matrix_memcpy(data2, data1);
+
+    gsl_test(gsl_extrap_richardson_vector(data1, 2, NULL, NULL, NULL),
+             "Unexpected error reported in %s");
+    const double result1 = gsl_matrix_get(data1, 0, 0);
+
+    gsl_test_abs(result1, 3.0, GSL_DBL_EPSILON,
+                 "%s scalar correct result", __func__);
+
+    gsl_vector * k = gsl_vector_alloc(1);
+    gsl_vector_set(k, 0, 1);
+    gsl_test(gsl_extrap_richardson_vector(data2, 2, k, NULL, NULL),
+             "Unexpected error reported in %s");
+    gsl_vector_free(k);
+    const double result2 = gsl_matrix_get(data2, 0, 0);
+
+    gsl_test_abs(result1, result2, GSL_DBL_EPSILON,
+                 "%s with k == NULL", __func__);
+  }
+
+  {
+    int i, j;
+    gsl_matrix_set(data1, 0, 0, 1.0);
+    gsl_matrix_set(data1, 0, 1, 2.0);
+    gsl_matrix_memcpy(data2, data1);
+
+    gsl_matrix * normtable1 = gsl_matrix_alloc(data1->size2, data1->size2);
+    gsl_matrix * normtable2 = gsl_matrix_alloc(data1->size2, data1->size2);
+
+    gsl_test(gsl_extrap_richardson_vector(
+               data1, 2, NULL, normtable1, NULL),
+             "Unexpected error reported in %s");
+    const double result1 = gsl_matrix_get(data1, 0, 0);
+
+    gsl_test_abs(result1, 3.0, GSL_DBL_EPSILON,
+                 "%s correct extrapolation answer");
+
+    for (i = 0; i < normtable1->size1 - 1; ++i)
+      {
+        for (j = i + 1; j < normtable1->size2; ++j)
+          {
+            gsl_test(!gsl_isnan(gsl_matrix_get(normtable1, i, j)),
+                     "%s expected NAN in normtable1 at (%d, %d)",
+                     __func__, i, j);
+          }
+      }
+
+    gsl_test_abs(gsl_matrix_get(normtable1, 0, 0), 1.0, GSL_DBL_EPSILON,
+                 "%s normtable1 (0,0) value", __func__);
+    gsl_test_abs(gsl_matrix_get(normtable1, 1, 0), 2.0, GSL_DBL_EPSILON,
+                 "%s normtable1 (1,0) value", __func__);
+    gsl_test_abs(gsl_matrix_get(normtable1, 1, 1), 3.0, GSL_DBL_EPSILON,
+                 "%s normtable1 (1,1) value", __func__);
+
+    gsl_vector * exact = gsl_vector_alloc(1);
+    gsl_vector_set_zero(exact);
+    gsl_test(gsl_extrap_richardson_vector(
+               data2, 2.0, NULL, normtable2, exact),
+             "Unexpected error reported in %s");
+    gsl_vector_free(exact);
+    const double result2 = gsl_matrix_get(data2, 0, 0);
+
+    for (i = 0; i < normtable2->size1 - 1; ++i)
+      {
+        for (j = i + 1; j < normtable2->size2; ++j)
+          {
+            gsl_test(!gsl_isnan(gsl_matrix_get(normtable2, i, j)),
+                     "%s expected NAN in normtable2 at (%d, %d)",
+                     __func__, i, j);
+          }
+      }
+
+    gsl_test_abs(result1, result2, GSL_DBL_EPSILON,
+                 "%s with exact == NULL", __func__);
+
+    for (i = 0; i < normtable1->size1; ++i)
+      {
+        for (j = 0; j < i + 1; ++j)
+          {
+            gsl_test_abs(gsl_matrix_get(normtable1, i, j),
+                         gsl_matrix_get(normtable2, i, j),
+                         GSL_DBL_EPSILON,
+                         "%s identical normtable results at (%d, %d)",
+                         __func__, i, j);
+          }
+      }
+
+    gsl_matrix_free(normtable2);
+    gsl_matrix_free(normtable1);
+  }
+
+  gsl_matrix_free(data2);
+  gsl_matrix_free(data1);
 }
 
 void
-test_richardson_extrapolation_twolevels()
+test_richardson_vector_extrapolation_twolevels()
 {
-    gsl_matrix * data2 = gsl_matrix_alloc(1, 2);
-    gsl_matrix * data3 = gsl_matrix_alloc(1, 3);
-    gsl_vector * k1 = gsl_vector_alloc(1);
-    gsl_vector * k2 = gsl_vector_alloc(2);
-
-    {
-        gsl_matrix_set(data2, 0, 0, 1.0);
-        gsl_matrix_set(data2, 0, 1, 2.0);
-        gsl_vector_set(k1, 0, 1);
-
-        gsl_test(gsl_extrap_richardson(data2, 2, k1, NULL, NULL),
-                "Unexpected error reported in %s");
-        gsl_test_abs(gsl_matrix_get(data2, 0, 0), 3.0, GSL_DBL_EPSILON,
-                "%s scalar correct result at %s:%d",
-                __func__, __FILE__, __LINE__);
-    }
-
-    {
-        gsl_matrix_set(data2, 0, 0, 2.0);
-        gsl_matrix_set(data2, 0, 1, 3.0);
-        gsl_vector_set(k1, 0, 1);
-
-        gsl_test(gsl_extrap_richardson(data2, 2, k1, NULL, NULL),
-                "Unexpected error reported in %s");
-        gsl_test_abs(gsl_matrix_get(data2, 0, 0), 4.0, GSL_DBL_EPSILON,
-                "%s scalar correct result at %s:%d",
-                __func__, __FILE__, __LINE__);
-    }
-
-    {
-        gsl_matrix_set(data2, 0, 0, 3.0);
-        gsl_matrix_set(data2, 0, 1, 4.0);
-        gsl_vector_set(k1, 0, 2);
-
-        gsl_test(gsl_extrap_richardson(data2, 2, k1, NULL, NULL),
-                "Unexpected error reported in %s");
-        gsl_test_abs(gsl_matrix_get(data2, 0, 0), 13.0/3.0, GSL_DBL_EPSILON,
-                "%s scalar correct result at %s:%d",
-                __func__, __FILE__, __LINE__);
-    }
-
-    {
-        gsl_matrix_set(data3, 0, 0, 1.0);
-        gsl_matrix_set(data3, 0, 1, 2.0);
-        gsl_matrix_set(data3, 0, 2, 3.0);
-        gsl_vector_set(k1, 0, 1);
-
-        gsl_test(gsl_extrap_richardson(data3, 2, k1, NULL, NULL),
-                "Unexpected error reported in %s");
-        gsl_test_abs(gsl_matrix_get(data3, 0, 0), 13.0/3.0, GSL_DBL_EPSILON,
-                "%s scalar correct result at %s:%d",
-                __func__, __FILE__, __LINE__);
-    }
-
-    {
-        int i, j;
-        gsl_matrix * normtable = gsl_matrix_alloc(data3->size2, data3->size2);
-
-        gsl_matrix_set(data3, 0, 0, 1.0);
-        gsl_matrix_set(data3, 0, 1, 2.0);
-        gsl_matrix_set(data3, 0, 2, 3.0);
-        gsl_vector_set(k2, 0, 1);
-        gsl_vector_set(k2, 1, 2);
-
-        gsl_test(gsl_extrap_richardson(data3, 2, k2, normtable, NULL),
-                "Unexpected error reported in %s");
-        gsl_test_abs(gsl_matrix_get(data3, 0, 0), 13.0/3.0, GSL_DBL_EPSILON,
-                "%s scalar correct result at %s:%d",
-                __func__, __FILE__, __LINE__);
-
-        for (i = 0; i < normtable->size1 - 1; ++i) {
-            for (j = i+1; j < normtable->size2; ++j) {
-                gsl_test(!gsl_isnan(gsl_matrix_get(normtable, i, j)),
-                        "%s expected NAN in normtable at (%d, %d)",
-                        __func__, i, j);
-            }
-        }
-
-        gsl_test_abs(gsl_matrix_get(normtable, 0, 0), 1.0, GSL_DBL_EPSILON,
-                "%s normtable (0,0) value", __func__);
-        gsl_test_abs(gsl_matrix_get(normtable, 1, 0), 2.0, GSL_DBL_EPSILON,
-                "%s normtable (1,0) value", __func__);
-        gsl_test_abs(gsl_matrix_get(normtable, 2, 0), 3.0, GSL_DBL_EPSILON,
-                "%s normtable (2,0) value", __func__);
-        gsl_test_abs(gsl_matrix_get(normtable, 1, 1), 3.0, GSL_DBL_EPSILON,
-                "%s normtable (1,1) value", __func__);
-        gsl_test_abs(gsl_matrix_get(normtable, 2, 1), 4.0, GSL_DBL_EPSILON,
-                "%s normtable (2,1) value", __func__);
-        gsl_test_abs(gsl_matrix_get(normtable, 2, 2), 13.0/3.0, GSL_DBL_EPSILON,
-                "%s normtable (2,2) value", __func__);
-
-        gsl_matrix_free(normtable);
-    }
-
-    {
-        gsl_matrix_set(data3, 0, 0, 1.0);
-        gsl_matrix_set(data3, 0, 1, 2.0);
-        gsl_matrix_set(data3, 0, 2, 3.0);
-        gsl_vector_set(k2, 0, 2);
-        gsl_vector_set(k2, 1, 3);
-
-        gsl_test(gsl_extrap_richardson(data3, 2, k2, NULL, NULL),
-                "Unexpected error reported in %s");
-        gsl_test_abs(gsl_matrix_get(data3, 0, 0), 73.0/21.0, GSL_DBL_EPSILON,
-                "%s scalar correct result at %s:%d",
-                __func__, __FILE__, __LINE__);
-    }
-
-    {
-        gsl_matrix_set(data3, 0, 0, 1.0);
-        gsl_matrix_set(data3, 0, 1, 2.0);
-        gsl_matrix_set(data3, 0, 2, 3.0);
-        gsl_vector_set(k1, 0, 2);
-
-        gsl_test(gsl_extrap_richardson(data3, 2, k1, NULL, NULL),
-                "Unexpected error reported in %s");
-        gsl_test_abs(gsl_matrix_get(data3, 0, 0), 73.0/21.0, GSL_DBL_EPSILON,
-                "%s scalar correct result at %s:%d",
-                __func__, __FILE__, __LINE__);
-    }
-
-
-    gsl_vector_free(k2);
-    gsl_vector_free(k1);
-    gsl_matrix_free(data3);
-    gsl_matrix_free(data2);
+  gsl_matrix * data2 = gsl_matrix_alloc(1, 2);
+  gsl_matrix * data3 = gsl_matrix_alloc(1, 3);
+  gsl_vector * k1 = gsl_vector_alloc(1);
+  gsl_vector * k2 = gsl_vector_alloc(2);
+
+  {
+    gsl_matrix_set(data2, 0, 0, 1.0);
+    gsl_matrix_set(data2, 0, 1, 2.0);
+    gsl_vector_set(k1, 0, 1);
+
+    gsl_test(gsl_extrap_richardson_vector(data2, 2, k1, NULL, NULL),
+             "Unexpected error reported in %s");
+    gsl_test_abs(gsl_matrix_get(data2, 0, 0), 3.0, GSL_DBL_EPSILON,
+                 "%s scalar correct result at %s:%d",
+                 __func__, __FILE__, __LINE__);
+  }
+
+  {
+    gsl_matrix_set(data2, 0, 0, 2.0);
+    gsl_matrix_set(data2, 0, 1, 3.0);
+    gsl_vector_set(k1, 0, 1);
+
+    gsl_test(gsl_extrap_richardson_vector(data2, 2, k1, NULL, NULL),
+             "Unexpected error reported in %s");
+    gsl_test_abs(gsl_matrix_get(data2, 0, 0), 4.0, GSL_DBL_EPSILON,
+                 "%s scalar correct result at %s:%d",
+                 __func__, __FILE__, __LINE__);
+  }
+
+  {
+    gsl_matrix_set(data2, 0, 0, 3.0);
+    gsl_matrix_set(data2, 0, 1, 4.0);
+    gsl_vector_set(k1, 0, 2);
+
+    gsl_test(gsl_extrap_richardson_vector(data2, 2, k1, NULL, NULL),
+             "Unexpected error reported in %s");
+    gsl_test_abs(gsl_matrix_get(data2, 0, 0), 13.0 / 3.0, GSL_DBL_EPSILON,
+                 "%s scalar correct result at %s:%d",
+                 __func__, __FILE__, __LINE__);
+  }
+
+  {
+    gsl_matrix_set(data3, 0, 0, 1.0);
+    gsl_matrix_set(data3, 0, 1, 2.0);
+    gsl_matrix_set(data3, 0, 2, 3.0);
+    gsl_vector_set(k1, 0, 1);
+
+    gsl_test(gsl_extrap_richardson_vector(data3, 2, k1, NULL, NULL),
+             "Unexpected error reported in %s");
+    gsl_test_abs(gsl_matrix_get(data3, 0, 0), 13.0 / 3.0, GSL_DBL_EPSILON,
+                 "%s scalar correct result at %s:%d",
+                 __func__, __FILE__, __LINE__);
+  }
+
+  {
+    int i, j;
+    gsl_matrix * normtable = gsl_matrix_alloc(data3->size2, data3->size2);
+
+    gsl_matrix_set(data3, 0, 0, 1.0);
+    gsl_matrix_set(data3, 0, 1, 2.0);
+    gsl_matrix_set(data3, 0, 2, 3.0);
+    gsl_vector_set(k2, 0, 1);
+    gsl_vector_set(k2, 1, 2);
+
+    gsl_test(gsl_extrap_richardson_vector(data3, 2, k2, normtable, NULL),
+             "Unexpected error reported in %s");
+    gsl_test_abs(gsl_matrix_get(data3, 0, 0), 13.0 / 3.0, GSL_DBL_EPSILON,
+                 "%s scalar correct result at %s:%d",
+                 __func__, __FILE__, __LINE__);
+
+    for (i = 0; i < normtable->size1 - 1; ++i)
+      {
+        for (j = i + 1; j < normtable->size2; ++j)
+          {
+            gsl_test(!gsl_isnan(gsl_matrix_get(normtable, i, j)),
+                     "%s expected NAN in normtable at (%d, %d)",
+                     __func__, i, j);
+          }
+      }
+
+    gsl_test_abs(gsl_matrix_get(normtable, 0, 0), 1.0, GSL_DBL_EPSILON,
+                 "%s normtable (0,0) value", __func__);
+    gsl_test_abs(gsl_matrix_get(normtable, 1, 0), 2.0, GSL_DBL_EPSILON,
+                 "%s normtable (1,0) value", __func__);
+    gsl_test_abs(gsl_matrix_get(normtable, 2, 0), 3.0, GSL_DBL_EPSILON,
+                 "%s normtable (2,0) value", __func__);
+    gsl_test_abs(gsl_matrix_get(normtable, 1, 1), 3.0, GSL_DBL_EPSILON,
+                 "%s normtable (1,1) value", __func__);
+    gsl_test_abs(gsl_matrix_get(normtable, 2, 1), 4.0, GSL_DBL_EPSILON,
+                 "%s normtable (2,1) value", __func__);
+    gsl_test_abs(gsl_matrix_get(normtable, 2, 2), 13.0 / 3.0, GSL_DBL_EPSILON,
+                 "%s normtable (2,2) value", __func__);
+
+    gsl_matrix_free(normtable);
+  }
+
+  {
+    gsl_matrix_set(data3, 0, 0, 1.0);
+    gsl_matrix_set(data3, 0, 1, 2.0);
+    gsl_matrix_set(data3, 0, 2, 3.0);
+    gsl_vector_set(k2, 0, 2);
+    gsl_vector_set(k2, 1, 3);
+
+    gsl_test(gsl_extrap_richardson_vector(data3, 2, k2, NULL, NULL),
+             "Unexpected error reported in %s");
+    gsl_test_abs(gsl_matrix_get(data3, 0, 0), 73.0 / 21.0,
+                 GSL_DBL_EPSILON*10,
+                 "%s scalar correct result at %s:%d",
+                 __func__, __FILE__, __LINE__);
+  }
+
+  {
+    gsl_matrix_set(data3, 0, 0, 1.0);
+    gsl_matrix_set(data3, 0, 1, 2.0);
+    gsl_matrix_set(data3, 0, 2, 3.0);
+    gsl_vector_set(k1, 0, 2);
+
+    gsl_test(gsl_extrap_richardson_vector(data3, 2, k1, NULL, NULL),
+             "Unexpected error reported in %s");
+    gsl_test_abs(gsl_matrix_get(data3, 0, 0), 73.0 / 21.0,
+                 GSL_DBL_EPSILON*10,
+                 "%s scalar correct result at %s:%d",
+                 __func__, __FILE__, __LINE__);
+  }
+
+
+  gsl_vector_free(k2);
+  gsl_vector_free(k1);
+  gsl_matrix_free(data3);
+  gsl_matrix_free(data2);
 
 }
 
 void
-test_richardson_extrapolation_multiplelevels()
+test_richardson_vector_extrapolation_multiplelevels()
+{
+  gsl_matrix * data2 = gsl_matrix_alloc(1, 2);
+  gsl_matrix * data4 = gsl_matrix_alloc(1, 4);
+  gsl_vector * k1 = gsl_vector_alloc(1);
+  gsl_vector * k2 = gsl_vector_alloc(2);
+  gsl_vector * k3 = gsl_vector_alloc(3);
+
+  {
+    int i, j;
+    gsl_vector_set(k1, 0, 3);
+
+    gsl_matrix_set(data2, 0, 0, 1.0);
+    gsl_matrix_set(data2, 0, 1, 2.0);
+    gsl_test(gsl_extrap_richardson_vector(data2, 2, k1, NULL, NULL),
+             "Unexpected error reported in %s");
+    const double xx = gsl_matrix_get(data2, 0, 0);
+
+    gsl_matrix_set(data2, 0, 0, 2.0);
+    gsl_matrix_set(data2, 0, 1, 3.0);
+    gsl_test(gsl_extrap_richardson_vector(data2, 2, k1, NULL, NULL),
+             "Unexpected error reported in %s");
+    const double xy = gsl_matrix_get(data2, 0, 0);
+
+    gsl_matrix_set(data2, 0, 0, 3.0);
+    gsl_matrix_set(data2, 0, 1, 4.0);
+    gsl_test(gsl_extrap_richardson_vector(data2, 2, k1, NULL, NULL),
+             "Unexpected error reported in %s");
+    const double xz = gsl_matrix_get(data2, 0, 0);
+
+    gsl_vector_set(k1, 0, 6);
+
+    gsl_matrix_set(data2, 0, 0, xx);
+    gsl_matrix_set(data2, 0, 1, xy);
+    gsl_test(gsl_extrap_richardson_vector(data2, 2, k1, NULL, NULL),
+             "Unexpected error reported in %s");
+    const double yx = gsl_matrix_get(data2, 0, 0);
+
+    gsl_matrix_set(data2, 0, 0, xy);
+    gsl_matrix_set(data2, 0, 1, xz);
+    gsl_test(gsl_extrap_richardson_vector(data2, 2, k1, NULL, NULL),
+             "Unexpected error reported in %s");
+    const double yy = gsl_matrix_get(data2, 0, 0);
+
+    gsl_vector_set(k1, 0, 9);
+
+    gsl_matrix_set(data2, 0, 0, yx);
+    gsl_matrix_set(data2, 0, 1, yy);
+    gsl_test(gsl_extrap_richardson_vector(data2, 2, k1, NULL, NULL),
+             "Unexpected error reported in %s");
+    const double zx = gsl_matrix_get(data2, 0, 0);
+
+    gsl_vector_set(k3, 0, 3);
+    gsl_vector_set(k3, 1, 6);
+    gsl_vector_set(k3, 2, 9);
+
+    gsl_matrix_set(data4, 0, 0, 1.0);
+    gsl_matrix_set(data4, 0, 1, 2.0);
+    gsl_matrix_set(data4, 0, 2, 3.0);
+    gsl_matrix_set(data4, 0, 3, 4.0);
+    gsl_test(gsl_extrap_richardson_vector(data4, 2, k3, NULL, NULL),
+             "Unexpected error reported in %s");
+    gsl_test_abs(gsl_matrix_get(data4, 0, 0), zx, GSL_DBL_EPSILON,
+                 "%s scalar correct result at %s:%d",
+                 __func__, __FILE__, __LINE__);
+
+    gsl_matrix * normtable = gsl_matrix_alloc(data4->size2, data4->size2);
+
+    gsl_vector_set(k2, 0, 3);
+    gsl_vector_set(k2, 1, 6);
+
+    gsl_matrix_set(data4, 0, 0, 1.0);
+    gsl_matrix_set(data4, 0, 1, 2.0);
+    gsl_matrix_set(data4, 0, 2, 3.0);
+    gsl_matrix_set(data4, 0, 3, 4.0);
+    gsl_test(gsl_extrap_richardson_vector(data4, 2, k2, normtable, NULL),
+             "Unexpected error reported in %s");
+    gsl_test_abs(gsl_matrix_get(data4, 0, 0), zx, GSL_DBL_EPSILON,
+                 "%s scalar correct result at %s:%d",
+                 __func__, __FILE__, __LINE__);
+
+    for (i = 0; i < normtable->size1 - 1; ++i)
+      {
+        for (j = i + 1; j < normtable->size2; ++j)
+          {
+            gsl_test(!gsl_isnan(gsl_matrix_get(normtable, i, j)),
+                     "%s expected NAN in normtable at (%d, %d)",
+                     __func__, i, j);
+          }
+      }
+
+    gsl_test_abs(gsl_matrix_get(normtable, 0, 0), 1.0, GSL_DBL_EPSILON,
+                 "%s normtable (0,0) value", __func__);
+    gsl_test_abs(gsl_matrix_get(normtable, 1, 0), 2.0, GSL_DBL_EPSILON,
+                 "%s normtable (1,0) value", __func__);
+    gsl_test_abs(gsl_matrix_get(normtable, 2, 0), 3.0, GSL_DBL_EPSILON,
+                 "%s normtable (2,0) value", __func__);
+    gsl_test_abs(gsl_matrix_get(normtable, 3, 0), 4.0, GSL_DBL_EPSILON,
+                 "%s normtable (3,0) value", __func__);
+    gsl_test_abs(gsl_matrix_get(normtable, 1, 1), xx, GSL_DBL_EPSILON,
+                 "%s normtable (1,1) value", __func__);
+    gsl_test_abs(gsl_matrix_get(normtable, 2, 1), xy, GSL_DBL_EPSILON,
+                 "%s normtable (2,1) value", __func__);
+    gsl_test_abs(gsl_matrix_get(normtable, 3, 1), xz, GSL_DBL_EPSILON,
+                 "%s normtable (3,1) value", __func__);
+    gsl_test_abs(gsl_matrix_get(normtable, 2, 2), yx, GSL_DBL_EPSILON,
+                 "%s normtable (2,2) value", __func__);
+    gsl_test_abs(gsl_matrix_get(normtable, 3, 2), yy, GSL_DBL_EPSILON,
+                 "%s normtable (3,2) value", __func__);
+    gsl_test_abs(gsl_matrix_get(normtable, 3, 3), zx, GSL_DBL_EPSILON,
+                 "%s normtable (3,3) value", __func__);
+
+    gsl_matrix_free(normtable);
+  }
+
+  gsl_vector_free(k3);
+  gsl_vector_free(k2);
+  gsl_vector_free(k1);
+  gsl_matrix_free(data4);
+  gsl_matrix_free(data2);
+}
+
+void
+test_richardson_vector_extrapolation_vectorinputs()
+{
+  {
+    gsl_matrix * data = gsl_matrix_alloc(2, 2);
+    gsl_matrix_set(data, 0, 0, 1.0);
+    gsl_matrix_set(data, 0, 1, 2.0);
+    gsl_matrix_set(data, 1, 0, 2.0);
+    gsl_matrix_set(data, 1, 1, 3.0);
+
+    gsl_matrix * normtable = gsl_matrix_alloc(data->size2, data->size2);
+
+    gsl_vector * exact = gsl_vector_alloc(data->size2);
+    gsl_vector_set(exact, 0, 3.0);
+    gsl_vector_set(exact, 1, 4.0);
+
+    gsl_test(gsl_extrap_richardson_vector(
+               data, 2, NULL, normtable, exact),
+             "Unexpected error reported in %s");
+    gsl_test_abs(gsl_matrix_get(data, 0, 0), 3.0, GSL_DBL_EPSILON,
+                 "%s data (0,0) value at %s:%d", __func__, __FILE__, __LINE__);
+    gsl_test_abs(gsl_matrix_get(data, 1, 0), 4.0, GSL_DBL_EPSILON,
+                 "%s data (1,0) value at %s:%d", __func__, __FILE__, __LINE__);
+
+    gsl_test_abs(gsl_matrix_get(normtable, 1, 1), 0.0, GSL_DBL_EPSILON,
+                 "%s normtable (1,1) value at %s:%d",
+                 __func__, __FILE__, __LINE__);
+
+    gsl_vector_free(exact);
+    gsl_matrix_free(normtable);
+    gsl_matrix_free(data);
+  }
+
+  {
+    gsl_matrix * data = gsl_matrix_alloc(2, 3);
+    gsl_matrix_set(data, 0, 0, 1.0);
+    gsl_matrix_set(data, 0, 1, 2.0);
+    gsl_matrix_set(data, 0, 2, 3.0);
+    gsl_matrix_set(data, 1, 0, 2.0);
+    gsl_matrix_set(data, 1, 1, 3.0);
+    gsl_matrix_set(data, 1, 2, 4.0);
+
+    gsl_test(gsl_extrap_richardson_vector(data, 2, NULL, NULL, NULL),
+             "Unexpected error reported in %s");
+    gsl_test_abs(gsl_matrix_get(data, 0, 0), 13.0 / 3.0, GSL_DBL_EPSILON,
+                 "%s data (0,0) value at %s:%d", __func__, __FILE__, __LINE__);
+    gsl_test_abs(gsl_matrix_get(data, 1, 0), 16.0 / 3.0, GSL_DBL_EPSILON,
+                 "%s data (1,0) value at %s:%d", __func__, __FILE__, __LINE__);
+
+    gsl_matrix_free(data);
+  }
+}
+
+void
+test_richardson_extrapolation_step()
 {
-    gsl_matrix * data2 = gsl_matrix_alloc(1, 2);
-    gsl_matrix * data4 = gsl_matrix_alloc(1, 4);
-    gsl_vector * k1 = gsl_vector_alloc(1);
-    gsl_vector * k2 = gsl_vector_alloc(2);
-    gsl_vector * k3 = gsl_vector_alloc(3);
-
-    {
-        int i, j;
-        gsl_vector_set(k1, 0, 3);
-
-        gsl_matrix_set(data2, 0, 0, 1.0);
-        gsl_matrix_set(data2, 0, 1, 2.0);
-        gsl_test(gsl_extrap_richardson(data2, 2, k1, NULL, NULL),
-                "Unexpected error reported in %s");
-        const double xx = gsl_matrix_get(data2, 0, 0);
-
-        gsl_matrix_set(data2, 0, 0, 2.0);
-        gsl_matrix_set(data2, 0, 1, 3.0);
-        gsl_test(gsl_extrap_richardson(data2, 2, k1, NULL, NULL),
-                "Unexpected error reported in %s");
-        const double xy = gsl_matrix_get(data2, 0, 0);
-
-        gsl_matrix_set(data2, 0, 0, 3.0);
-        gsl_matrix_set(data2, 0, 1, 4.0);
-        gsl_test(gsl_extrap_richardson(data2, 2, k1, NULL, NULL),
-                "Unexpected error reported in %s");
-        const double xz = gsl_matrix_get(data2, 0, 0);
-
-        gsl_vector_set(k1, 0, 6);
-
-        gsl_matrix_set(data2, 0, 0, xx);
-        gsl_matrix_set(data2, 0, 1, xy);
-        gsl_test(gsl_extrap_richardson(data2, 2, k1, NULL, NULL),
-                "Unexpected error reported in %s");
-        const double yx = gsl_matrix_get(data2, 0, 0);
-
-        gsl_matrix_set(data2, 0, 0, xy);
-        gsl_matrix_set(data2, 0, 1, xz);
-        gsl_test(gsl_extrap_richardson(data2, 2, k1, NULL, NULL),
-                "Unexpected error reported in %s");
-        const double yy = gsl_matrix_get(data2, 0, 0);
-
-        gsl_vector_set(k1, 0, 9);
-
-        gsl_matrix_set(data2, 0, 0, yx);
-        gsl_matrix_set(data2, 0, 1, yy);
-        gsl_test(gsl_extrap_richardson(data2, 2, k1, NULL, NULL),
-                "Unexpected error reported in %s");
-        const double zx = gsl_matrix_get(data2, 0, 0);
-
-        gsl_vector_set(k3, 0, 3);
-        gsl_vector_set(k3, 1, 6);
-        gsl_vector_set(k3, 2, 9);
-
-        gsl_matrix_set(data4, 0, 0, 1.0);
-        gsl_matrix_set(data4, 0, 1, 2.0);
-        gsl_matrix_set(data4, 0, 2, 3.0);
-        gsl_matrix_set(data4, 0, 3, 4.0);
-        gsl_test(gsl_extrap_richardson(data4, 2, k3, NULL, NULL),
-                "Unexpected error reported in %s");
-        gsl_test_abs(gsl_matrix_get(data4, 0, 0), zx, GSL_DBL_EPSILON,
-                "%s scalar correct result at %s:%d",
-                __func__, __FILE__, __LINE__);
-
-        gsl_matrix * normtable = gsl_matrix_alloc(data4->size2, data4->size2);
-
-        gsl_vector_set(k2, 0, 3);
-        gsl_vector_set(k2, 1, 6);
-
-        gsl_matrix_set(data4, 0, 0, 1.0);
-        gsl_matrix_set(data4, 0, 1, 2.0);
-        gsl_matrix_set(data4, 0, 2, 3.0);
-        gsl_matrix_set(data4, 0, 3, 4.0);
-        gsl_test(gsl_extrap_richardson(data4, 2, k2, normtable, NULL),
-                "Unexpected error reported in %s");
-        gsl_test_abs(gsl_matrix_get(data4, 0, 0), zx, GSL_DBL_EPSILON,
-                "%s scalar correct result at %s:%d",
-                __func__, __FILE__, __LINE__);
-
-        for (i = 0; i < normtable->size1 - 1; ++i) {
-            for (j = i+1; j < normtable->size2; ++j) {
-                gsl_test(!gsl_isnan(gsl_matrix_get(normtable, i, j)),
-                        "%s expected NAN in normtable at (%d, %d)",
-                        __func__, i, j);
-            }
-        }
-
-        gsl_test_abs(gsl_matrix_get(normtable, 0, 0), 1.0, GSL_DBL_EPSILON,
-                "%s normtable (0,0) value", __func__);
-        gsl_test_abs(gsl_matrix_get(normtable, 1, 0), 2.0, GSL_DBL_EPSILON,
-                "%s normtable (1,0) value", __func__);
-        gsl_test_abs(gsl_matrix_get(normtable, 2, 0), 3.0, GSL_DBL_EPSILON,
-                "%s normtable (2,0) value", __func__);
-        gsl_test_abs(gsl_matrix_get(normtable, 3, 0), 4.0, GSL_DBL_EPSILON,
-                "%s normtable (3,0) value", __func__);
-        gsl_test_abs(gsl_matrix_get(normtable, 1, 1), xx, GSL_DBL_EPSILON,
-                "%s normtable (1,1) value", __func__);
-        gsl_test_abs(gsl_matrix_get(normtable, 2, 1), xy, GSL_DBL_EPSILON,
-                "%s normtable (2,1) value", __func__);
-        gsl_test_abs(gsl_matrix_get(normtable, 3, 1), xz, GSL_DBL_EPSILON,
-                "%s normtable (3,1) value", __func__);
-        gsl_test_abs(gsl_matrix_get(normtable, 2, 2), yx, GSL_DBL_EPSILON,
-                "%s normtable (2,2) value", __func__);
-        gsl_test_abs(gsl_matrix_get(normtable, 3, 2), yy, GSL_DBL_EPSILON,
-                "%s normtable (3,2) value", __func__);
-        gsl_test_abs(gsl_matrix_get(normtable, 3, 3), zx, GSL_DBL_EPSILON,
-                "%s normtable (3,3) value", __func__);
-
-        gsl_matrix_free(normtable);
-    }
-
-    gsl_vector_free(k3);
-    gsl_vector_free(k2);
-    gsl_vector_free(k1);
-    gsl_matrix_free(data4);
-    gsl_matrix_free(data2);
+  const int    ki           = 1;
+
+  {
+    const double t   = 2.0;
+    double       Ah  = 1.0;
+    const double Aht = 2.0;
+
+    gsl_test(gsl_extrap_richardson_step(&Ah, &Aht, t, ki),
+             "Unexpected error reported in %s for t=%f", __func__, t);
+    gsl_test_abs(Ah, 3.0, GSL_DBL_EPSILON,
+                 "Extrapolation in %s for ki=%d, t=%f", __func__, ki, t);
+  }
+
+  {
+    const double t = 3.0;
+    double       Ah  = 1.0;
+    const double Aht = 2.0;
+    gsl_test(gsl_extrap_richardson_step(&Ah, &Aht, t, ki),
+             "Unexpected error reported in %s for t=%f", __func__, t);
+    gsl_test_abs(Ah, 5.0 / 2.0, GSL_DBL_EPSILON,
+                 "Extrapolation in %s for ki=%d, t=%f", __func__, ki, t);
+  }
 }
 
 void
-test_richardson_extrapolation_vectorinputs()
+test_richardson_extrapolation_defaults()
 {
-    {
-        gsl_matrix * data = gsl_matrix_alloc(2,2);
-        gsl_matrix_set(data, 0, 0, 1.0);
-        gsl_matrix_set(data, 0, 1, 2.0);
-        gsl_matrix_set(data, 1, 0, 2.0);
-        gsl_matrix_set(data, 1, 1, 3.0);
-
-        gsl_matrix * normtable = gsl_matrix_alloc(data->size2, data->size2);
-
-        gsl_vector * exact = gsl_vector_alloc(data->size2);
-        gsl_vector_set(exact, 0, 3.0);
-        gsl_vector_set(exact, 1, 4.0);
-
-        gsl_test(gsl_extrap_richardson(
-                    data, 2, NULL, normtable, exact),
-                "Unexpected error reported in %s");
-        gsl_test_abs(gsl_matrix_get(data, 0, 0), 3.0, GSL_DBL_EPSILON,
-                "%s data (0,0) value at %s:%d", __func__, __FILE__, __LINE__);
-        gsl_test_abs(gsl_matrix_get(data, 1, 0), 4.0, GSL_DBL_EPSILON,
-                "%s data (1,0) value at %s:%d", __func__, __FILE__, __LINE__);
-
-        gsl_test_abs(gsl_matrix_get(normtable, 1, 1), 0.0, GSL_DBL_EPSILON,
-                "%s normtable (1,1) value at %s:%d",
-                __func__, __FILE__, __LINE__);
-
-        gsl_vector_free(exact);
-        gsl_matrix_free(normtable);
-        gsl_matrix_free(data);
-    }
-
-    {
-        gsl_matrix * data = gsl_matrix_alloc(2,3);
-        gsl_matrix_set(data, 0, 0, 1.0);
-        gsl_matrix_set(data, 0, 1, 2.0);
-        gsl_matrix_set(data, 0, 2, 3.0);
-        gsl_matrix_set(data, 1, 0, 2.0);
-        gsl_matrix_set(data, 1, 1, 3.0);
-        gsl_matrix_set(data, 1, 2, 4.0);
-
-        gsl_test(gsl_extrap_richardson(data, 2, NULL, NULL, NULL),
-                "Unexpected error reported in %s");
-        gsl_test_abs(gsl_matrix_get(data, 0, 0), 13.0/3.0, GSL_DBL_EPSILON,
-                "%s data (0,0) value at %s:%d", __func__, __FILE__, __LINE__);
-        gsl_test_abs(gsl_matrix_get(data, 1, 0), 16.0/3.0, GSL_DBL_EPSILON,
-                "%s data (1,0) value at %s:%d", __func__, __FILE__, __LINE__);
-
-        gsl_matrix_free(data);
-    }
+  gsl_vector * data1 = gsl_vector_alloc(2);
+  gsl_vector * data2 = gsl_vector_alloc(data1->size);
+  {
+    gsl_vector_set(data1, 0, 1.0);
+    gsl_vector_set(data1, 1, 2.0);
+
+    gsl_test(gsl_extrap_richardson(data1, 2, NULL, NULL, 0),
+             "Unexpected error reported in %s");
+    gsl_test_abs(gsl_vector_get(data1, 0), 3.0, GSL_DBL_EPSILON,
+                 "%s scalar correct result", __func__);
+  }
+
+  {
+    gsl_vector_set(data1, 0, 1.0);
+    gsl_vector_set(data1, 1, 2.0);
+
+    gsl_test(gsl_extrap_richardson(data1, 3, NULL, NULL, 0),
+             "Unexpected error reported in %s");
+    gsl_test_abs(gsl_vector_get(data1, 0), 5.0 / 2.0, GSL_DBL_EPSILON,
+                 "%s scalar correct result", __func__);
+  }
+
+  {
+    gsl_vector_set(data1, 0, 1.0);
+    gsl_vector_set(data1, 1, 2.0);
+    gsl_vector_memcpy(data2, data1);
+
+    gsl_test(gsl_extrap_richardson(data1, 2, NULL, NULL, 0),
+             "Unexpected error reported in %s");
+    const double result1 = gsl_vector_get(data1, 0);
+
+    gsl_test_abs(result1, 3.0, GSL_DBL_EPSILON,
+                 "%s scalar correct result", __func__);
+
+    gsl_vector * k = gsl_vector_alloc(1);
+    gsl_vector_set(k, 0, 1);
+    gsl_test(gsl_extrap_richardson(data2, 2, k, NULL, 0),
+             "Unexpected error reported in %s");
+    gsl_vector_free(k);
+    const double result2 = gsl_vector_get(data2, 0);
+
+    gsl_test_abs(result1, result2, GSL_DBL_EPSILON,
+                 "%s with k == NULL", __func__);
+  }
+
+  {
+    int i, j;
+    gsl_vector_set(data1, 0, 1.0);
+    gsl_vector_set(data1, 1, 2.0);
+    gsl_vector_memcpy(data2, data1);
+
+    gsl_matrix * normtable1 = gsl_matrix_alloc(data1->size, data1->size);
+    gsl_matrix * normtable2 = gsl_matrix_alloc(data1->size, data1->size);
+
+    gsl_test(gsl_extrap_richardson(data1, 2, NULL, normtable1, 0),
+             "Unexpected error reported in %s");
+    const double result1 = gsl_vector_get(data1, 0);
+
+    gsl_test_abs(result1, 3.0, GSL_DBL_EPSILON,
+                 "%s correct extrapolation answer");
+
+    for (i = 0; i < normtable1->size1 - 1; ++i)
+      {
+        for (j = i + 1; j < normtable1->size2; ++j)
+          {
+            gsl_test(!gsl_isnan(gsl_matrix_get(normtable1, i, j)),
+                     "%s expected NAN in normtable1 at (%d, %d)",
+                     __func__, i, j);
+          }
+      }
+
+    gsl_test_abs(gsl_matrix_get(normtable1, 0, 0), 1.0, GSL_DBL_EPSILON,
+                 "%s normtable1 (0,0) value", __func__);
+    gsl_test_abs(gsl_matrix_get(normtable1, 1, 0), 2.0, GSL_DBL_EPSILON,
+                 "%s normtable1 (1,0) value", __func__);
+    gsl_test_abs(gsl_matrix_get(normtable1, 1, 1), 3.0, GSL_DBL_EPSILON,
+                 "%s normtable1 (1,1) value", __func__);
+
+    const double exact = 0.0;
+    gsl_test(gsl_extrap_richardson(data2, 2.0, NULL, normtable2, exact),
+             "Unexpected error reported in %s");
+    const double result2 = gsl_vector_get(data2, 0);
+
+    for (i = 0; i < normtable2->size1 - 1; ++i)
+      {
+        for (j = i + 1; j < normtable2->size2; ++j)
+          {
+            gsl_test(!gsl_isnan(gsl_matrix_get(normtable2, i, j)),
+                     "%s expected NAN in normtable2 at (%d, %d)",
+                     __func__, i, j);
+          }
+      }
+
+    gsl_test_abs(result1, result2, GSL_DBL_EPSILON,
+                 "%s with exact == NULL", __func__);
+
+    for (i = 0; i < normtable1->size1; ++i)
+      {
+        for (j = 0; j < i + 1; ++j)
+          {
+            gsl_test_abs(gsl_matrix_get(normtable1, i, j),
+                         gsl_matrix_get(normtable2, i, j),
+                         GSL_DBL_EPSILON,
+                         "%s identical normtable results at (%d, %d)",
+                         __func__, i, j);
+          }
+      }
+
+    gsl_matrix_free(normtable2);
+    gsl_matrix_free(normtable1);
+  }
+
+  gsl_vector_free(data2);
+  gsl_vector_free(data1);
+}
+
+void
+test_richardson_extrapolation_twolevels()
+{
+  gsl_vector * data2 = gsl_vector_alloc(2);
+  gsl_vector * data3 = gsl_vector_alloc(3);
+  gsl_vector * k1 = gsl_vector_alloc(1);
+  gsl_vector * k2 = gsl_vector_alloc(2);
+
+  {
+    gsl_vector_set(data2, 0, 1.0);
+    gsl_vector_set(data2, 1, 2.0);
+    gsl_vector_set(k1, 0, 1);
+
+    gsl_test(gsl_extrap_richardson(data2, 2, k1, NULL, 0),
+             "Unexpected error reported in %s");
+    gsl_test_abs(gsl_vector_get(data2, 0), 3.0, GSL_DBL_EPSILON,
+                 "%s scalar correct result at %s:%d",
+                 __func__, __FILE__, __LINE__);
+  }
+
+  {
+    gsl_vector_set(data2, 0, 2.0);
+    gsl_vector_set(data2, 1, 3.0);
+    gsl_vector_set(k1, 0, 1);
+
+    gsl_test(gsl_extrap_richardson(data2, 2, k1, NULL, 0),
+             "Unexpected error reported in %s");
+    gsl_test_abs(gsl_vector_get(data2, 0), 4.0, GSL_DBL_EPSILON,
+                 "%s scalar correct result at %s:%d",
+                 __func__, __FILE__, __LINE__);
+  }
+
+  {
+    gsl_vector_set(data2, 0, 3.0);
+    gsl_vector_set(data2, 1, 4.0);
+    gsl_vector_set(k1, 0, 2);
+
+    gsl_test(gsl_extrap_richardson(data2, 2, k1, NULL, 0),
+             "Unexpected error reported in %s");
+    gsl_test_abs(gsl_vector_get(data2, 0), 13.0 / 3.0, GSL_DBL_EPSILON,
+                 "%s scalar correct result at %s:%d",
+                 __func__, __FILE__, __LINE__);
+  }
+
+  {
+    gsl_vector_set(data3, 0, 1.0);
+    gsl_vector_set(data3, 1, 2.0);
+    gsl_vector_set(data3, 2, 3.0);
+    gsl_vector_set(k1, 0, 1);
+
+    gsl_test(gsl_extrap_richardson(data3, 2, k1, NULL, 0),
+             "Unexpected error reported in %s");
+    gsl_test_abs(gsl_vector_get(data3, 0), 13.0 / 3.0, GSL_DBL_EPSILON,
+                 "%s scalar correct result at %s:%d",
+                 __func__, __FILE__, __LINE__);
+  }
+
+  {
+    int i, j;
+    gsl_matrix * normtable = gsl_matrix_alloc(data3->size, data3->size);
+
+    gsl_vector_set(data3, 0, 1.0);
+    gsl_vector_set(data3, 1, 2.0);
+    gsl_vector_set(data3, 2, 3.0);
+    gsl_vector_set(k2, 0, 1);
+    gsl_vector_set(k2, 1, 2);
+
+    gsl_test(gsl_extrap_richardson(data3, 2, k2, normtable, 0),
+             "Unexpected error reported in %s");
+    gsl_test_abs(gsl_vector_get(data3, 0), 13.0 / 3.0, GSL_DBL_EPSILON,
+                 "%s scalar correct result at %s:%d",
+                 __func__, __FILE__, __LINE__);
+
+    for (i = 0; i < normtable->size1 - 1; ++i)
+      {
+        for (j = i + 1; j < normtable->size2; ++j)
+          {
+            gsl_test(!gsl_isnan(gsl_matrix_get(normtable, i, j)),
+                     "%s expected NAN in normtable at (%d, %d)",
+                     __func__, i, j);
+          }
+      }
+
+    gsl_test_abs(gsl_matrix_get(normtable, 0, 0), 1.0, GSL_DBL_EPSILON,
+                 "%s normtable (0,0) value", __func__);
+    gsl_test_abs(gsl_matrix_get(normtable, 1, 0), 2.0, GSL_DBL_EPSILON,
+                 "%s normtable (1,0) value", __func__);
+    gsl_test_abs(gsl_matrix_get(normtable, 2, 0), 3.0, GSL_DBL_EPSILON,
+                 "%s normtable (2,0) value", __func__);
+    gsl_test_abs(gsl_matrix_get(normtable, 1, 1), 3.0, GSL_DBL_EPSILON,
+                 "%s normtable (1,1) value", __func__);
+    gsl_test_abs(gsl_matrix_get(normtable, 2, 1), 4.0, GSL_DBL_EPSILON,
+                 "%s normtable (2,1) value", __func__);
+    gsl_test_abs(gsl_matrix_get(normtable, 2, 2), 13.0 / 3.0, GSL_DBL_EPSILON,
+                 "%s normtable (2,2) value", __func__);
+
+    gsl_matrix_free(normtable);
+  }
+
+  {
+    gsl_vector_set(data3, 0, 1.0);
+    gsl_vector_set(data3, 1, 2.0);
+    gsl_vector_set(data3, 2, 3.0);
+    gsl_vector_set(k2, 0, 2);
+    gsl_vector_set(k2, 1, 3);
+
+    gsl_test(gsl_extrap_richardson(data3, 2, k2, NULL, 0),
+             "Unexpected error reported in %s");
+    gsl_test_abs(gsl_vector_get(data3, 0), 73.0 / 21.0,
+                 GSL_DBL_EPSILON*10,
+                 "%s scalar correct result at %s:%d",
+                 __func__, __FILE__, __LINE__);
+  }
+
+  {
+    gsl_vector_set(data3, 0, 1.0);
+    gsl_vector_set(data3, 1, 2.0);
+    gsl_vector_set(data3, 2, 3.0);
+    gsl_vector_set(k1, 0, 2);
+
+    gsl_test(gsl_extrap_richardson(data3, 2, k1, NULL, 0),
+             "Unexpected error reported in %s");
+    gsl_test_abs(gsl_vector_get(data3, 0), 73.0 / 21.0,
+                 GSL_DBL_EPSILON*10,
+                 "%s scalar correct result at %s:%d",
+                 __func__, __FILE__, __LINE__);
+  }
+
+
+  gsl_vector_free(k2);
+  gsl_vector_free(k1);
+  gsl_vector_free(data3);
+  gsl_vector_free(data2);
+
+}
+
+void
+test_richardson_extrapolation_multiplelevels()
+{
+  gsl_vector * data2 = gsl_vector_alloc(2);
+  gsl_vector * data4 = gsl_vector_alloc(4);
+  gsl_vector * k1 = gsl_vector_alloc(1);
+  gsl_vector * k2 = gsl_vector_alloc(2);
+  gsl_vector * k3 = gsl_vector_alloc(3);
+
+  {
+    int i, j;
+    gsl_vector_set(k1, 0, 3);
+
+    gsl_vector_set(data2, 0, 1.0);
+    gsl_vector_set(data2, 1, 2.0);
+    gsl_test(gsl_extrap_richardson(data2, 2, k1, NULL, 0),
+             "Unexpected error reported in %s");
+    const double xx = gsl_vector_get(data2, 0);
+
+    gsl_vector_set(data2, 0, 2.0);
+    gsl_vector_set(data2, 1, 3.0);
+    gsl_test(gsl_extrap_richardson(data2, 2, k1, NULL, 0),
+             "Unexpected error reported in %s");
+    const double xy = gsl_vector_get(data2, 0);
+
+    gsl_vector_set(data2, 0, 3.0);
+    gsl_vector_set(data2, 1, 4.0);
+    gsl_test(gsl_extrap_richardson(data2, 2, k1, NULL, 0),
+             "Unexpected error reported in %s");
+    const double xz = gsl_vector_get(data2, 0);
+
+    gsl_vector_set(k1, 0, 6);
+
+    gsl_vector_set(data2, 0, xx);
+    gsl_vector_set(data2, 1, xy);
+    gsl_test(gsl_extrap_richardson(data2, 2, k1, NULL, 0),
+             "Unexpected error reported in %s");
+    const double yx = gsl_vector_get(data2, 0);
+
+    gsl_vector_set(data2, 0, xy);
+    gsl_vector_set(data2, 1, xz);
+    gsl_test(gsl_extrap_richardson(data2, 2, k1, NULL, 0),
+             "Unexpected error reported in %s");
+    const double yy = gsl_vector_get(data2, 0);
+
+    gsl_vector_set(k1, 0, 9);
+
+    gsl_vector_set(data2, 0, yx);
+    gsl_vector_set(data2, 1, yy);
+    gsl_test(gsl_extrap_richardson(data2, 2, k1, NULL, 0),
+             "Unexpected error reported in %s");
+    const double zx = gsl_vector_get(data2, 0);
+
+    gsl_vector_set(k3, 0, 3);
+    gsl_vector_set(k3, 1, 6);
+    gsl_vector_set(k3, 2, 9);
+
+    gsl_vector_set(data4, 0, 1.0);
+    gsl_vector_set(data4, 1, 2.0);
+    gsl_vector_set(data4, 2, 3.0);
+    gsl_vector_set(data4, 3, 4.0);
+    gsl_test(gsl_extrap_richardson(data4, 2, k3, NULL, 0),
+             "Unexpected error reported in %s");
+    gsl_test_abs(gsl_vector_get(data4, 0), zx, GSL_DBL_EPSILON,
+                 "%s scalar correct result at %s:%d",
+                 __func__, __FILE__, __LINE__);
+
+    gsl_matrix * normtable = gsl_matrix_alloc(data4->size, data4->size);
+
+    gsl_vector_set(k2, 0, 3);
+    gsl_vector_set(k2, 1, 6);
+
+    gsl_vector_set(data4, 0, 1.0);
+    gsl_vector_set(data4, 1, 2.0);
+    gsl_vector_set(data4, 2, 3.0);
+    gsl_vector_set(data4, 3, 4.0);
+    gsl_test(gsl_extrap_richardson(data4, 2, k2, normtable, 0),
+             "Unexpected error reported in %s");
+    gsl_test_abs(gsl_vector_get(data4, 0), zx, GSL_DBL_EPSILON,
+                 "%s scalar correct result at %s:%d",
+                 __func__, __FILE__, __LINE__);
+
+    for (i = 0; i < normtable->size1 - 1; ++i)
+      {
+        for (j = i + 1; j < normtable->size2; ++j)
+          {
+            gsl_test(!gsl_isnan(gsl_matrix_get(normtable, i, j)),
+                     "%s expected NAN in normtable at (%d, %d)",
+                     __func__, i, j);
+          }
+      }
+
+    gsl_test_abs(gsl_matrix_get(normtable, 0, 0), 1.0, GSL_DBL_EPSILON,
+                 "%s normtable (0,0) value", __func__);
+    gsl_test_abs(gsl_matrix_get(normtable, 1, 0), 2.0, GSL_DBL_EPSILON,
+                 "%s normtable (1,0) value", __func__);
+    gsl_test_abs(gsl_matrix_get(normtable, 2, 0), 3.0, GSL_DBL_EPSILON,
+                 "%s normtable (2,0) value", __func__);
+    gsl_test_abs(gsl_matrix_get(normtable, 3, 0), 4.0, GSL_DBL_EPSILON,
+                 "%s normtable (3,0) value", __func__);
+    gsl_test_abs(gsl_matrix_get(normtable, 1, 1), xx, GSL_DBL_EPSILON,
+                 "%s normtable (1,1) value", __func__);
+    gsl_test_abs(gsl_matrix_get(normtable, 2, 1), xy, GSL_DBL_EPSILON,
+                 "%s normtable (2,1) value", __func__);
+    gsl_test_abs(gsl_matrix_get(normtable, 3, 1), xz, GSL_DBL_EPSILON,
+                 "%s normtable (3,1) value", __func__);
+    gsl_test_abs(gsl_matrix_get(normtable, 2, 2), yx, GSL_DBL_EPSILON,
+                 "%s normtable (2,2) value", __func__);
+    gsl_test_abs(gsl_matrix_get(normtable, 3, 2), yy, GSL_DBL_EPSILON,
+                 "%s normtable (3,2) value", __func__);
+    gsl_test_abs(gsl_matrix_get(normtable, 3, 3), zx, GSL_DBL_EPSILON,
+                 "%s normtable (3,3) value", __func__);
+
+    gsl_matrix_free(normtable);
+  }
+
+  gsl_vector_free(k3);
+  gsl_vector_free(k2);
+  gsl_vector_free(k1);
+  gsl_vector_free(data4);
+  gsl_vector_free(data2);
 }
 
 int
 main(int argc, char **argv)
 {
-    gsl_ieee_env_setup();
+  gsl_ieee_env_setup();
+
+  test_richardson_vector_extrapolation_step();
+  test_richardson_vector_extrapolation_defaults();
+  test_richardson_vector_extrapolation_twolevels();
+  test_richardson_vector_extrapolation_multiplelevels();
+  test_richardson_vector_extrapolation_vectorinputs();
 
-    test_richardson_extrapolation_step();
-    test_richardson_extrapolation_defaults();
-    test_richardson_extrapolation_twolevels();
-    test_richardson_extrapolation_multiplelevels();
-    test_richardson_extrapolation_vectorinputs();
+  test_richardson_extrapolation_step();
+  test_richardson_extrapolation_defaults();
+  test_richardson_extrapolation_twolevels();
+  test_richardson_extrapolation_multiplelevels();
 
-    exit(gsl_test_summary());
+  exit(gsl_test_summary());
 }
-- 
1.6.0.4

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

[Brian Gough]

At Tue, 15 Dec 2009 11:31:10 -0600,
Rhys Ulerich wrote:
> You're right that each component could be extrapolated independently.
> I wrote it using vectors because doing so allows using BLAS calls for
> the linear algebra and speeds up many component use cases.  I'd prefer
> to keep the code vector-capable under the covers under the theory that
> people extrapolating only a scalar at a time aren't all that worried
> about speed.

Ideally I'd like to start with a scalar version using normal C arrays,
similar to the gsl_sum functions, for implementation simplicity.
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.

-- 
Brian Gough

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

[Rhys Ulerich]

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

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

[Brian Gough]

At Wed, 16 Dec 2009 18:43:44 -0600,
Rhys Ulerich wrote:
> > 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?

Where possible it is good to have an API that does not depend on
gsl_vector and gsl_matrix as it makes it simpler to use (or reuse) the
code from C programs with a minimal dependence on GSL.  Implementing
each routine to be minimal (without extra features) also makes life
simpler.

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

I will send that once I've looked at it a bit more.

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

[Brian Gough]

At Thu, 17 Dec 2009 22:50:22 +0000,
Brian Gough wrote:
> > 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?
> 
> I will send that once I've looked at it a bit more.

I think the best interface would be the analog of the existing levin
sum routines:

int gsl_sum_levin_u_accel (const double *array,
                           const size_t n,
                           gsl_sum_levin_u_workspace * w,
                           double *sum_accel, double *abserr);

e.g.

int gsl_extrap_richardson (const double *array,
                           const size_t n,
                           const double t,
                           const size_t k
                           gsl_extrap_richardson_workspace * w,
                           double *result, double *abserr);

It would be useful to have both a result and abserr, to take into
account cancellation error.

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

[Rhys Ulerich]

Sorry for the long hiatus on this Richardson extrapolation thread:

> I think the best interface would be the analog of the existing levin
> sum routines:
>
> ...
>
> int gsl_extrap_richardson (const double *array,
>                           const size_t n,
>                           const double t,
>                           const size_t k
>                           gsl_extrap_richardson_workspace * w,
>                           double *result, double *abserr);
>
> It would be useful to have both a result and abserr, to take into
> account cancellation error.

I'm good on moving the interface away from gsl_matrix and gsl_vector.
No problem there.  But, from looking a bit at the levin sum
documentation, I'm not sure how to apply the abserr concept to
Richardson extrapolation.  Any guidance would be much appreciated.

Brian, if this is an easier thing to talk through on the phone and
you're willing to do so, please contact me directly.

- Rhys

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

[Brian Gough]

At Sat, 6 Feb 2010 09:25:47 -0600,
Rhys Ulerich wrote:
> I'm good on moving the interface away from gsl_matrix and gsl_vector.
> No problem there.  But, from looking a bit at the levin sum
> documentation, I'm not sure how to apply the abserr concept to
> Richardson extrapolation.  Any guidance would be much appreciated.

The error term in the Levin sum is calculated by simple error
propagation assuming the individual terms are accurate to double
precision.  Therefore it is just a measure of the cumulative
cancellation error. Since there are a lot of subtractions, it's worth
tracking this.

-- 
Brian Gough