[gsl/temp_odecrash] temporary diagnostic
Frantisek Kluknavsky
fkluknav at fedoraproject.org
Tue Jan 29 15:25:19 UTC 2013
commit 9747dd763d020f85fdc60f259f6823afa99f9faf
Author: Frantisek Kluknavsky <fkluknav at redhat.com>
Date: Tue Jan 29 12:51:40 2013 +0100
temporary diagnostic
wrk.patch | 685 ++++++++++++++++++++++++++++++++++++++++++++-----------------
1 files changed, 493 insertions(+), 192 deletions(-)
---
diff --git a/wrk.patch b/wrk.patch
index 4c755f7..445c1cc 100644
--- a/wrk.patch
+++ b/wrk.patch
@@ -162,244 +162,545 @@ diff -up wrk/ode-initval2/evolve.c.wrk wrk/ode-initval2/evolve.c
diff -up wrk/ode-initval2/msbdf.c.wrk wrk/ode-initval2/msbdf.c
--- wrk/ode-initval2/msbdf.c.wrk 2013-01-09 10:35:45.259960403 +0100
-+++ wrk/ode-initval2/msbdf.c 2013-01-23 18:17:00.641056195 +0100
-@@ -43,6 +43,8 @@
- framework.
- */
++++ wrk/ode-initval2/msbdf.c 2013-01-28 10:42:35.000000000 +0100
+@@ -82,6 +82,7 @@ typedef struct
+ size_t *ordprevbackup; /* backup of ordprev */
+ double *errlev; /* desired error level of y */
+ gsl_vector *abscor; /* absolute y values for correction */
++ gsl_vector *abscorscaled; /* scaled abscor for order evaluation */
+ gsl_vector *relcor; /* relative y values for correction */
+ gsl_vector *svec; /* saved abscor & work area */
+ gsl_vector *tempvec; /* work area */
+@@ -91,7 +92,6 @@ typedef struct
+ gsl_matrix *M; /* Newton iteration matrix */
+ gsl_permutation *p; /* permutation for LU decomposition of M */
+ gsl_vector *rhs; /* right hand side equations (-G) */
+-
+ long int ni; /* stepper call counter */
+ size_t ord; /* current order of method */
+ double tprev; /* t point of previous call */
+@@ -446,6 +446,33 @@ msbdf_alloc (size_t dim)
+ GSL_ERROR_NULL ("failed to allocate space for rhs", GSL_ENOMEM);
+ }
-+#define PD(fl) printf("%lu\n", *((unsigned long *)&(fl)))
++ state->abscorscaled = gsl_vector_alloc (dim);
+
- #include <config.h>
- #include <stdlib.h>
- #include <string.h>
-@@ -116,7 +118,6 @@ static void *
- msbdf_alloc (size_t dim)
- {
- msbdf_state_t *state = (msbdf_state_t *) malloc (sizeof (msbdf_state_t));
--
- if (state == 0)
- {
- GSL_ERROR_NULL ("failed to allocate space for msbdf_state", GSL_ENOMEM);
-@@ -611,6 +612,15 @@ msbdf_calccoeffs (const size_t ord, cons
- printf ("errcoeff=%.5e\n", *errcoeff);
- #endif
++ if (state->abscorscaled == 0)
++ {
++ gsl_vector_free (state->rhs);
++ gsl_permutation_free (state->p);
++ gsl_matrix_free (state->M);
++ free (state->dfdt);
++ gsl_matrix_free (state->dfdy);
++ gsl_vector_free (state->tempvec);
++ gsl_vector_free (state->svec);
++ gsl_vector_free (state->relcor);
++ gsl_vector_free (state->abscor);
++ free (state->errlev);
++ free (state->ordprevbackup);
++ free (state->ordprev);
++ free (state->hprevbackup);
++ free (state->hprev);
++ free (state->l);
++ free (state->ytmp2);
++ free (state->ytmp);
++ free (state->zbackup);
++ free (state->z);
++ free (state);
++ GSL_ERROR_NULL ("failed to allocate space for abscorscaled", GSL_ENOMEM);
++ }
++
+ msbdf_reset ((void *) state, dim);
-+/*PD(l[0]);
-+PD(l[1]);
-+PD(*errcoeff);
-+PD(*ordm1coeff);
-+PD(*ordp1coeff);
-+PD(*ordp2coeff);
-+PD(*gamma);
-+printf("\n");
-+*/
- return GSL_SUCCESS;
+ state->driver = NULL;
+@@ -926,14 +953,14 @@ msbdf_corrector (void *vstate, const gsl
}
-@@ -1007,9 +1017,9 @@ msbdf_eval_order (gsl_vector * abscor, g
- if (ordm1est > ordest && ordm1est > ordp1est && ordm1est > min_incr)
- {
- *ord -= 1;
--#ifdef DEBUG
-+//#ifdef DEBUG
- printf ("-- eval_order order DECREASED to %d\n", (int) *ord);
--#endif
-+//#endif
- }
-
- else if (ordp1est > ordest && ordp1est > ordm1est && ordp1est > min_incr)
-@@ -1082,7 +1092,8 @@ msbdf_apply (void *vstate, size_t dim, d
- const gsl_odeiv2_system * sys)
+ static int
+-msbdf_eval_order (gsl_vector * abscor, gsl_vector * tempvec,
++msbdf_eval_order (gsl_vector * abscorscaled, gsl_vector * tempvec,
+ gsl_vector * svec, const double errcoeff,
+ const size_t dim, const double errlev[],
+ const double ordm1coeff, const double ordp1coeff,
+ const double ordp1coeffprev, const double ordp2coeff,
+ const double hprev[],
+ const double h, const double z[],
+- size_t * ord, size_t * ordwait)
++ size_t * ord)
{
- /* Carries out a step by BDF linear multistep methods. */
--
-+//printf("zaciatok msbdf\n");
-+//PD(h);
- msbdf_state_t *state = (msbdf_state_t *) vstate;
+ /* Evaluates and executes change in method order (current, current-1
+ or current+1). Order which maximizes the step length is selected.
+@@ -953,7 +980,7 @@ msbdf_eval_order (gsl_vector * abscor, g
- double *const z = state->z;
-@@ -1109,7 +1120,7 @@ msbdf_apply (void *vstate, size_t dim, d
+ /* Relative step length estimate for current order */
- const size_t max_failcount = 3;
- size_t i;
--
-+printf("ord at msbdf start\n%d\n", ord);
- #ifdef DEBUG
- {
- size_t di;
-@@ -1317,6 +1328,7 @@ msbdf_apply (void *vstate, size_t dim, d
- msbdf_check_step_size_decrease (hprev))
- {
- state->ord--;
-+ printf("stability order decrease\n");
- state->ordwait = ord + 2;
- ord = state->ord;
+- ordest = 1.0 / (pow (bias * gsl_blas_dnrm2 (abscor) / sqrt (dim)
++ ordest = 1.0 / (pow (bias * gsl_blas_dnrm2 (abscorscaled) / sqrt (dim)
+ * errcoeff, 1.0 / (*ord + 1)) + safety);
-@@ -1329,7 +1341,7 @@ msbdf_apply (void *vstate, size_t dim, d
+ /* Relative step length estimate for order ord - 1 */
+@@ -983,7 +1010,7 @@ msbdf_eval_order (gsl_vector * abscor, g
+ for (i = 0; i < dim; i++)
+ {
+ gsl_vector_set (svec, i, gsl_vector_get (svec, i) * c +
+- gsl_vector_get (abscor, i));
++ gsl_vector_get (abscorscaled, i));
+ }
- {
- const int deltaord = ord - ordprev[0];
--
-+printf("ord sanity check\n%d\n%d\n", ord, ordprev[0]);
- if (deltaord > 1 || deltaord < -1)
- {
- printf ("-- order change %d\n", deltaord);
-@@ -1466,7 +1478,10 @@ msbdf_apply (void *vstate, size_t dim, d
- /* Calculate polynomial coefficients (l), error coefficient and
- auxiliary coefficients
- */
--
-+/*printf("%d\n", ord);
-+printf("%lu\n", state->ordwait);
-+PD(h);
-+PD(hprev[0]);*/
- msbdf_calccoeffs (ord, state->ordwait, h, hprev, l, &errcoeff,
- º1coeff, &ordp1coeff, &ordp2coeff, &gamma);
+ ordp1est = 1.0 / (pow (bias2 * gsl_blas_dnrm2 (svec) / sqrt (dim)
+@@ -1020,8 +1047,6 @@ msbdf_eval_order (gsl_vector * abscor, g
+ #endif
+ }
-@@ -1619,6 +1634,7 @@ msbdf_apply (void *vstate, size_t dim, d
+- *ordwait = *ord + 2;
+-
+ return GSL_SUCCESS;
+ }
- /* Consider and execute order change for next step */
+@@ -1096,6 +1121,7 @@ msbdf_apply (void *vstate, size_t dim, d
+ size_t *const ordprevbackup = state->ordprevbackup;
+ double *const errlev = state->errlev;
+ gsl_vector *const abscor = state->abscor;
++ gsl_vector *const abscorscaled = state->abscorscaled;
+ gsl_vector *const relcor = state->relcor;
+ gsl_vector *const svec = state->svec;
+ gsl_vector *const tempvec = state->tempvec;
+@@ -1591,13 +1617,16 @@ msbdf_apply (void *vstate, size_t dim, d
+ }
+ }
-+printf("ord before eval\n%d\n", state->ord);
- if (state->ordwait == 0)
- {
- msbdf_eval_order (abscor, tempvec, svec, errcoeff, dim, errlev,
-@@ -1626,7 +1642,7 @@ msbdf_apply (void *vstate, size_t dim, d
- state->ordp1coeffprev, ordp2coeff,
- hprev, h, z, &(state->ord), &(state->ordwait));
- }
--
-+printf("ord after eval\n%d\n", state->ord);
- /* Undo scaling of abscor for possible order increase on next step */
+- /* Scale abscor with errlev for later use in norm calculations */
++ /* Scale abscor with errlev for later use in norm calculations
++ in order evaluation in msbdf_eval_order
++ */
{
size_t i;
-@@ -1675,7 +1691,9 @@ msbdf_apply (void *vstate, size_t dim, d
- printf ("-- nJ=%d, nM=%d\n", (int) state->nJ, (int) state->nM);
- #endif
+
+ for (i = 0; i < dim; i++)
+ {
+- gsl_vector_set (abscor, i, gsl_vector_get (abscor, i) / errlev[i]);
++ gsl_vector_set (abscorscaled, i,
++ gsl_vector_get (abscor, i) / errlev[i]);
+ }
}
--
-+printf("koniec msbdf, ord:\n%d\n", ord);
-+printf("ord at msbdf end\n%d\n", state->ord);
-+//PD(h);
- return GSL_SUCCESS;
- }
+@@ -1613,27 +1642,28 @@ msbdf_apply (void *vstate, size_t dim, d
+
+ for (i = 0; i < dim; i++)
+ {
+- gsl_vector_set (svec, i, gsl_vector_get (abscor, i));
++ gsl_vector_set (svec, i, gsl_vector_get (abscorscaled, i));
+ }
+ }
+
+- /* Consider and execute order change for next step */
++ /* Consider and execute order change for next step if order is unchanged. */
+
+- if (state->ordwait == 0)
+- {
+- msbdf_eval_order (abscor, tempvec, svec, errcoeff, dim, errlev,
+- ordm1coeff, ordp1coeff,
+- state->ordp1coeffprev, ordp2coeff,
+- hprev, h, z, &(state->ord), &(state->ordwait));
+- }
++ if (state->ordwait == 0)
++ {
++ if (ord - ordprev[0] == 0)
++ {
++ msbdf_eval_order (abscorscaled, tempvec, svec, errcoeff, dim, errlev,
++ ordm1coeff, ordp1coeff,
++ state->ordp1coeffprev, ordp2coeff,
++ hprev, h, z, &(state->ord));
+
+- /* Undo scaling of abscor for possible order increase on next step */
+- {
+- size_t i;
+-
+- for (i = 0; i < dim; i++)
++ state->ordwait = state->ord + 2;
++ }
++ else
+ {
+- gsl_vector_set (abscor, i, gsl_vector_get (abscor, i) * errlev[i]);
++ /* Postpone order evaluation if order has been modified elsewhere */
++
++ state->ordwait = 2;
+ }
+ }
+
+@@ -1701,10 +1731,13 @@ msbdf_reset (void *vstate, size_t dim)
+ state->ordwaitbackup = 2;
+ state->failord = 0;
+ state->failt = GSL_NAN;
++ state->tprev = GSL_NAN;
+ state->gammaprev = 1.0;
++ state->gammaprevbackup = 1.0;
+ state->nJ = 0;
+ state->nM = 0;
+ state->failcount = 0;
++ state->ordp1coeffprev = 0.0;
+
+ DBL_ZERO_MEMSET (state->hprev, MSBDF_MAX_ORD);
+ DBL_ZERO_MEMSET (state->hprevbackup, MSBDF_MAX_ORD);
diff -up wrk/ode-initval2/test.c.wrk wrk/ode-initval2/test.c
--- wrk/ode-initval2/test.c.wrk 2013-01-09 10:48:22.051415928 +0100
-+++ wrk/ode-initval2/test.c 2013-01-23 16:15:54.368720700 +0100
-@@ -923,6 +923,18 @@ rhs_ringmod (double t, const double y[],
- f[13] = (-y[0] + uin1 - (ri + rg1) * y[13]) / ls1;
- f[14] = (-y[1] - (rc + rg1) * y[14]) / ls1;
-
-+ /*int i;
-+ printf("temp states\n");
-+ for (i=0; i<15; ++i) {
-+
-+ printf("%lu\n", *((unsigned long *)(&(y[i]))));
-+ }
++++ wrk/ode-initval2/test.c 2013-01-28 10:42:50.000000000 +0100
+@@ -163,7 +163,7 @@ rhs_xsin (double t, const double y[], do
+ static int n = 0, m = 0;
+ extern int nfe;
+ nfe += 1;
+-
+
-+ printf("%lu\n", *((unsigned long *)(&(rg2))));
-+ printf("%lu\n", *((unsigned long *)(&(ls2))));
-+ printf("%lu\n", *((unsigned long *)(&(rg3))));
-+ printf("%lu\n", *((unsigned long *)(&(ls3))));
-+*/
- return GSL_SUCCESS;
- }
+ if (rhs_xsin_reset)
+ {
+ rhs_xsin_reset = 0;
+@@ -1073,6 +1073,7 @@ test_odeiv_stepper (const gsl_odeiv2_ste
-@@ -1239,6 +1251,8 @@ test_evolve_system (const gsl_odeiv2_ste
- gsl_odeiv2_driver_free (d);
- }
+ {
+ int s = gsl_odeiv2_step_apply (d->s, t, h, y, yerr, 0, 0, sys);
++
+ if (s != GSL_SUCCESS)
+ {
+ gsl_test (s, "test_odeiv_stepper: %s step_apply returned %d", desc,
+@@ -1191,13 +1192,13 @@ test_evolve_system (const gsl_odeiv2_ste
+ if (s != GSL_SUCCESS)
+ {
+ /* check that t and y are unchanged */
+- gsl_test_abs (t, t_orig, 0.0, "%s, t must be restored on failure",
++ gsl_test_abs (t, t_orig, 0.0, "%s t must be restored on failure",
+ gsl_odeiv2_step_name (d->s));
-+#define PD(fl) printf("%lu\n", *((unsigned long *)&(fl)))
+ for (i = 0; i < sys->dimension; i++)
+ {
+ gsl_test_abs (y[i], y_orig[i], 0.0,
+- "%s, y must be restored on failure",
++ "%s y must be restored on failure",
+ gsl_odeiv2_step_name (d->s), desc, i);
+ }
+
+@@ -1259,11 +1260,11 @@ sys_driver (const gsl_odeiv2_step_type *
+ gsl_odeiv2_driver *d =
+ gsl_odeiv2_driver_alloc_standard_new (sys, T, h, epsabs, epsrel,
+ 1.0, 0.0);
+-
+
- int
- sys_driver (const gsl_odeiv2_step_type * T,
- const gsl_odeiv2_system * sys,
-@@ -1264,10 +1278,23 @@ sys_driver (const gsl_odeiv2_step_type *
+ extern int nfe, nje;
nfe = 0;
nje = 0;
-
-+ int poc=0;
+-
++
while (t < t1)
{
-+ // printf("pred evolve apply\n");
-+// PD(h);
s = gsl_odeiv2_evolve_apply (d->e, d->c, d->s, sys, &t, t1, &h, y);
+@@ -1570,7 +1571,8 @@ test_user_break (const gsl_odeiv2_step_t
-+//printf("po evolve apply\n");
-+ printf("\nkrok %d\n", ++poc);
-+//PD(h);
-+/* int k;
-+ for (k=0; k<15; ++k) {
-+ printf(" %lu\n",(*(unsigned long int *)(&(y[k]))));
-+ }
-+ printf(" %lu\n",(*(unsigned long int *)(&t)));
-+*/
-+ //printf("\n");
- #ifdef DEBUG
- printf ("%.5e %.5e %.5e %d\n", t, y[0], y[1],
- gsl_odeiv2_step_order (d->s));
-@@ -2003,10 +2030,11 @@ test_extreme_problems (void)
-
- for (i = 0; steppers[i] != 0; i++)
- {
-+ // printf("spustam driver p=%ld i=%ld\n", p, i);
- int s = sys_driver (steppers[i], prob[p], start[p], end[p],
- initstepsize[p], &y[sd[p] * i],
- epsabs[p], epsrel[p], probname[p]);
--
-+ //printf("koniec drivera\n");
- if (s != GSL_SUCCESS)
- {
- printf ("start=%.5e, initstepsize=%.5e\n", start[p],
-@@ -2018,7 +2046,7 @@ test_extreme_problems (void)
+ int s = gsl_odeiv2_driver_apply (d, &t, t1, y);
+
+- gsl_test ((s - GSL_EBADFUNC), "test_user_break returned %d", s);
++ gsl_test ((s - GSL_EBADFUNC), "%s test_user_break returned %d",
++ gsl_odeiv2_step_name (d->s), s);
+
+ gsl_odeiv2_driver_free (d);
+ }
+@@ -1938,17 +1940,15 @@ test_extreme_problems (void)
+ const size_t sd[] = { 4, 3, NRINGMOD };
- /* Compare results */
+ /* Integration interval for problems */
++
+ const double start[CONST_EXTREME_NPROB] = { 0.0 };
+- const double end[CONST_EXTREME_NPROB] = { 1e11, 1e11, 1e-3 };
++ const double end[CONST_EXTREME_NPROB] = { 1e11, 1e11, 1e-5 };
+
+ const double epsabs[CONST_EXTREME_NPROB] =
+- { 1e1 * GSL_DBL_MIN, 1e1 * GSL_DBL_MIN, 1e-7 };
+- const double epsrel[CONST_EXTREME_NPROB] = { 1e-12, 1e-12, 1e-7 };
++ { 1e1 * GSL_DBL_MIN, 1e1 * GSL_DBL_MIN, 1e-12 };
++ const double epsrel[CONST_EXTREME_NPROB] = { 1e-12, 1e-12, 1e-12 };
+ const double initstepsize[CONST_EXTREME_NPROB] = { 1e-5, 1e-5, 1e-10 };
+
+- /* Problem specific tolerance modification coefficients */
+- const double pec[CONST_EXTREME_NPROB] = { 1.0, 1.0, 1e1 };
+-
+ /* Steppers */
-- for (i = 1; steppers[i] != 0; i++)
-+ for (i = 0; steppers[i] != 0; i++)
- for (k = 0; k < sd[p]; k++)
- {
+ steppers[0] = gsl_odeiv2_step_bsimp;
+@@ -2024,7 +2024,7 @@ test_extreme_problems (void)
const double val1 = y[k];
-@@ -2450,7 +2478,7 @@ main (void)
+ const double val2 = y[sd[p] * i + k];
+ gsl_test_rel (val1, val2,
+- (pec[p] * GSL_MAX (err_target[0], err_target[i])),
++ (GSL_MAX (err_target[0], err_target[i])),
+ "%s/%s %s [%d]",
+ steppers[0]->name, steppers[i]->name,
+ probname[p], k);
+@@ -2033,57 +2033,81 @@ test_extreme_problems (void)
+ }
- /* Basic tests for all steppers */
+ void
+-test_driver (void)
++test_driver (const gsl_odeiv2_step_type * T)
+ {
+ /* Tests for gsl_odeiv2_driver object */
++
+ int s;
+ const double tol = 1e-8;
+ const double hmax = 1e-2;
+
+ double y[] = { 1.0, 0.0 };
+ double t = 0.0;
++ const double tinit = 0.0;
+ const double t1 = 8.25;
+ const double t2 = 100;
+ const double t3 = -2.5;
+ const size_t minsteps = ceil (t1 / hmax);
+ const size_t maxsteps = 20;
+ const double hmin = 1e-10;
+- const size_t nfsteps = 100;
+-
+- gsl_odeiv2_driver *d =
+- gsl_odeiv2_driver_alloc_y_new (&rhs_func_sin, gsl_odeiv2_step_rkf45,
+- 1e-3, tol, tol);
++
++ const unsigned long int nfsteps = 100000;
++ const double hfixed = 0.000025;
++
++ gsl_odeiv2_driver *d = gsl_odeiv2_driver_alloc_y_new (&rhs_func_sin, T,
++ 1e-3, tol, tol);
+ gsl_odeiv2_driver_set_hmax (d, hmax);
-- for (i = 0; p[i].type != 0; i++)
-+ /*for (i = 0; p[i].type != 0; i++)
+ s = gsl_odeiv2_driver_apply (d, &t, t1, y);
+
+ if (s != GSL_SUCCESS)
{
- test_stepper (p[i].type);
+- gsl_test (s, "test_driver apply returned %d", s);
++ gsl_test (s, "%s test_driver apply returned %d",
++ gsl_odeiv2_step_name (d->s), s);
+ }
+
+ /* Test that result is correct */
+
+- gsl_test_rel (y[0], cos (t1), d->n * tol, "test_driver y0 ");
+- gsl_test_rel (y[1], sin (t1), d->n * tol, "test_driver y1 ");
++ gsl_test_rel (y[0], cos (t1), d->n * tol, "%s test_driver y0",
++ gsl_odeiv2_step_name (d->s));
++ gsl_test_rel (y[1], sin (t1), d->n * tol, "%s test_driver y1",
++ gsl_odeiv2_step_name (d->s));
+
+ /* Test that maximum step length is obeyed */
+
+ if (d->n < minsteps)
+ {
+- gsl_test (1, "test_driver steps %d < minsteps %d \n", d->n, minsteps);
++ gsl_test (1, "%s test_driver steps %d < minsteps %d \n",
++ gsl_odeiv2_step_name (d->s), d->n, minsteps);
}
-@@ -2469,9 +2497,9 @@ main (void)
+ else
+ {
+- gsl_test (0, "test_driver max step length test");
++ gsl_test (0, "%s test_driver max step length test",
++ gsl_odeiv2_step_name (d->s));
+ }
+
++ /* Test changing integration direction from forward to backward */
++
++ gsl_odeiv2_driver_reset_hstart (d, -1e-3);
++
++ s = gsl_odeiv2_driver_apply (d, &t, tinit, y);
++
++ if (s != GSL_SUCCESS)
++ {
++ gsl_test (s, "%s test_driver apply returned %d",
++ gsl_odeiv2_step_name (d->s), s);
++ }
++
++ gsl_test_rel (y[0], cos (tinit), d->n * tol,
++ "%s test_driver y0", gsl_odeiv2_step_name (d->s));
++ gsl_test_rel (y[1], sin (tinit), d->n * tol,
++ "%s test_driver y1", gsl_odeiv2_step_name (d->s));
++
+ gsl_odeiv2_driver_free (d);
+
+ /* Test that maximum number of steps is obeyed */
+
+- d = gsl_odeiv2_driver_alloc_y_new (&rhs_func_sin, gsl_odeiv2_step_rkf45,
+- 1e-3, tol, tol);
++ d = gsl_odeiv2_driver_alloc_y_new (&rhs_func_sin, T, 1e-3, tol, tol);
+ gsl_odeiv2_driver_set_hmax (d, hmax);
+ gsl_odeiv2_driver_set_nmax (d, maxsteps);
+
+@@ -2091,19 +2115,20 @@ test_driver (void)
+
+ if (d->n != maxsteps + 1)
+ {
+- gsl_test (1, "test_driver steps %d, expected %d", d->n, maxsteps + 1);
++ gsl_test (1, "%s test_driver steps %d, expected %d",
++ gsl_odeiv2_step_name (d->s), d->n, maxsteps + 1);
+ }
+ else
+ {
+- gsl_test (0, "test_driver max steps test");
++ gsl_test (0, "%s test_driver max steps test",
++ gsl_odeiv2_step_name (d->s));
+ }
+
+ gsl_odeiv2_driver_free (d);
+
+ /* Test that minimum step length is obeyed */
+
+- d = gsl_odeiv2_driver_alloc_y_new (&rhs_func_broken, gsl_odeiv2_step_rk8pd,
+- 1e-3, tol, tol);
++ d = gsl_odeiv2_driver_alloc_y_new (&rhs_func_broken, T, 1e-3, tol, tol);
+
+ gsl_odeiv2_driver_set_hmin (d, hmin);
+ y[0] = 0.0;
+@@ -2113,77 +2138,87 @@ test_driver (void)
+
+ if (s != GSL_ENOPROG)
+ {
+- gsl_test (1, "test_driver min step test returned %d", s);
++ gsl_test (1, "%s test_driver min step test returned %d",
++ gsl_odeiv2_step_name (d->s), s);
+ }
+ else
+ {
+- gsl_test (0, "test_driver min step test");
++ gsl_test (0, "%s test_driver min step test",
++ gsl_odeiv2_step_name (d->s));
+ }
+
+ gsl_odeiv2_driver_free (d);
+
+ /* Test negative integration direction */
+
+- d = gsl_odeiv2_driver_alloc_y_new (&rhs_func_sin, gsl_odeiv2_step_rkf45,
+- -1e-3, tol, tol);
++ d = gsl_odeiv2_driver_alloc_y_new (&rhs_func_sin, T, -1e-3, tol, tol);
+
+ y[0] = 1.0;
+ y[1] = 0.0;
+ t = 2.5;
+
+- gsl_odeiv2_driver_set_nmax (d, 1000);
+ s = gsl_odeiv2_driver_apply (d, &t, t3, y);
+
+ {
+- const double tol = 1e-6;
++ const double tol = 1e-3;
+ const double test = fabs (t - t3);
+ const double val = fabs (sin (-5.0) - y[1]);
+
+ if (test > GSL_DBL_EPSILON)
+ {
+ gsl_test (1,
+- "test_driver negative dir diff %e, expected less than %e",
+- test, GSL_DBL_EPSILON);
++ "%s test_driver negative dir diff %e, expected less than %e",
++ gsl_odeiv2_step_name (d->s), test, GSL_DBL_EPSILON);
+ }
+ else if (val > tol)
+ {
+- gsl_test (1, "test_driver negative dir val %e, expected less than %e",
+- val, tol);
++ gsl_test (1,
++ "%s test_driver negative dir val %e, expected less than %e",
++ gsl_odeiv2_step_name (d->s), val, tol);
+ }
+ else
+ {
+- gsl_test (s, "test_driver negative direction test");
++ gsl_test (s, "%s test_driver negative direction test",
++ gsl_odeiv2_step_name (d->s));
+ }
+ }
+
+ /* Test driver_apply_fixed_step */
+
+- s = gsl_odeiv2_driver_apply_fixed_step (d, &t, 0.025, nfsteps, y);
++ gsl_odeiv2_driver_reset_hstart (d, 1e-3);
++ s = gsl_odeiv2_driver_apply_fixed_step (d, &t, hfixed, nfsteps, y);
++
++ if (s != GSL_SUCCESS)
++ {
++ gsl_test (s, "%s test_driver apply_fixed_step returned %d",
++ gsl_odeiv2_step_name (d->s), s);
++ }
+
+ {
+- const double tol = 1e-6;
++ const double tol = 1e-3;
+ const double val = fabs (sin (-2.5) - y[1]);
+
+ if (fabs (t) > nfsteps * GSL_DBL_EPSILON)
+ {
+ gsl_test (1,
+- "test_driver apply_fixed_step diff %e, expected less than %e",
+- fabs (t), nfsteps * GSL_DBL_EPSILON);
++ "%s test_driver apply_fixed_step t %e, expected less than %e",
++ gsl_odeiv2_step_name (d->s), fabs (t),
++ nfsteps * GSL_DBL_EPSILON);
+ }
+ else if (val > tol)
+ {
+ gsl_test (1,
+- "test_driver apply_fixed_step val %e, expected less than %e",
+- val, tol);
++ "%s test_driver apply_fixed_step val %e, expected less than %e",
++ gsl_odeiv2_step_name (d->s), val, tol);
+ }
+ else
+ {
+- gsl_test (s, "test_driver apply_fixed_step test");
++ gsl_test (s, "%s test_driver apply_fixed_step test",
++ gsl_odeiv2_step_name (d->s));
+ }
+ }
+
+ gsl_odeiv2_driver_free (d);
+-
+ }
+
+ void
+@@ -2407,7 +2442,7 @@ main (void)
+ {
+
+ /* Benchmark routine to compare stepper performance */
+- /* benchmark_precision(); return 0;*/
++ /* benchmark_precision(); return 0; */
+
+ /* Test single problem, for debugging purposes */
+ /* test_evolve_temp (gsl_odeiv2_step_msadams, 1e-3, 1e-8); return 0; */
+@@ -2453,10 +2488,6 @@ main (void)
+ for (i = 0; p[i].type != 0; i++)
+ {
+ test_stepper (p[i].type);
+- }
+-
+- for (i = 0; p[i].type != 0; i++)
+- {
+ test_evolve_linear (p[i].type, p[i].h, 1e-10);
+ test_evolve_exp (p[i].type, p[i].h, 1e-5);
+ test_evolve_sin (p[i].type, p[i].h, 1e-8);
+@@ -2468,6 +2499,7 @@ main (void)
+ test_broken (p[i].type, p[i].h, 1e-8);
test_stepsize_fail (p[i].type, p[i].h);
test_user_break (p[i].type, p[i].h);
++ test_driver (p[i].type);
}
--
-+*/
+
/* Derivative test for explicit steppers */
--
-+/*
- explicit_stepper[0].type = gsl_odeiv2_step_rk4;
- explicit_stepper[0].h = 1.0e-3;
- explicit_stepper[1].type = gsl_odeiv2_step_rk2;
-@@ -2491,16 +2519,16 @@ main (void)
- test_stepfn (explicit_stepper[i].type);
- test_stepfn2 (explicit_stepper[i].type);
+@@ -2493,14 +2525,12 @@ main (void)
}
--
-+*/
+
/* Special tests */
-
-- test_nonstiff_problems ();
-+ // test_nonstiff_problems ();
+-
++
+ test_nonstiff_problems ();
-- test_stiff_problems ();
-+ // test_stiff_problems ();
+ test_stiff_problems ();
test_extreme_problems ();
- test_driver ();
-+ //test_driver ();
-
+-
exit (gsl_test_summary ());
}
More information about the scm-commits
mailing list