This is the mail archive of the
gsl-discuss@sourceware.org
mailing list for the GSL project.
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