This is the mail archive of the gsl-discuss@sourceware.org mailing list for the GSL project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]
Other format: [Raw text]

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


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


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]