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

[Rhys Ulerich <rhys.ulerich@gmail.com>]

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