From 0e6d3460d5cac6777a940b07a5f29288f2addab1 Mon Sep 17 00:00:00 2001 From: Rhys Ulerich 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 + + * added scalar interface and tests by wrapping vector routines + 2009-12-13 Rhys Ulerich * 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 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 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