diff -up wrk/ode-initval2/bsimp.c.wrk wrk/ode-initval2/bsimp.c diff -U0 wrk/ode-initval2/ChangeLog.wrk wrk/ode-initval2/ChangeLog --- wrk/ode-initval2/ChangeLog.wrk 2013-01-29 13:40:22.470035141 +0100 +++ wrk/ode-initval2/ChangeLog 2013-01-29 12:58:21.000000000 +0100 @@ -0,0 +1,52 @@ +2013-01-27 Tuomo Keskitalo + + * msbdf.c: Corrected bug which enabled order to be changed + by two (first by stability enhancement, then by order evaluation + after a rejected step). Thanks for Frantisek Kluknavsky and + Andreas Schwab for bug reports! + + * msbdf.c: Added more state variables explicitly to be reset in + msbdf_reset. + + *msbdf.c: Added abscorscaled to remove division of abscor for + use in msbdf_eval_order and subsequent backwards multiplication. + + *test.c: test_extreme_problems: Increased ringmod case tolerances + from 1e-7 to 1e-12 to increase precision. Because that + increases computational burden, I also decreased end time + from 1e-3 to 1e-5. That decreased the acidity of the test + significantly, but the test case is now more appropriate + for normal "make check". Thanks for Frantisek Kluknavsky + for pointing out bad design of the test case. + +2012-09-10 Rhys Ulerich + + * test.c: Correct two out-of-order declarations. + Thanks to Brian Gladman for spotting the problem. + +2012-03-10 Tuomo Keskitalo + + * driver.c: Added function gsl_odeiv2_driver_reset_hstart to + allow resetting step size, possibly change its sign, too. + Check for non-zero hstart in initializations. + + * test.c: Modified test_driver to carry out tests with all + steppers. Small printout changes. + +2012-01-21 Tuomo Keskitalo + + * rkf45.c: Bug correction: Set gives_exact_dydt_out to 1 in + rkf45_type + +2012-01-14 Tuomo Keskitalo + + * evolve.c: Modified initial derivative evaluation in evolve_apply + to reuse previously calculated values. This saves calls to the + user function. Thanks for Illes Farkas for pointing out redundant + function calls! + +2011-06-28 Brian Gough + + * rk4imp.c (rk4imp_apply): use M_SQRT3 instead of sqrt(3) in array + initialiser. + diff -up wrk/ode-initval2/control.c.wrk wrk/ode-initval2/control.c diff -up wrk/ode-initval2/control_utils.c.wrk wrk/ode-initval2/control_utils.c diff -up wrk/ode-initval2/cscal.c.wrk wrk/ode-initval2/cscal.c diff -up wrk/ode-initval2/cstd.c.wrk wrk/ode-initval2/cstd.c diff -up wrk/ode-initval2/driver.c.wrk wrk/ode-initval2/driver.c --- wrk/ode-initval2/driver.c.wrk 2013-01-29 13:40:22.508035371 +0100 +++ wrk/ode-initval2/driver.c 2013-01-29 13:07:29.000000000 +0100 @@ -81,7 +81,7 @@ driver_alloc (const gsl_odeiv2_system * GSL_ERROR_NULL ("failed to allocate evolve object", GSL_ENOMEM); } - if (hstart >= 0.0 || hstart < 0.0) + if (hstart > 0.0 || hstart < 0.0) { state->h = hstart; } @@ -459,6 +459,29 @@ gsl_odeiv2_driver_reset (gsl_odeiv2_driv return GSL_SUCCESS; } +int +gsl_odeiv2_driver_reset_hstart (gsl_odeiv2_driver * d, const double hstart) +{ + /* Resets current driver and sets initial step size to hstart */ + + gsl_odeiv2_driver_reset (d); + + if ((d->hmin > fabs (hstart)) || (fabs (hstart) > d->hmax)) + { + GSL_ERROR_NULL ("hmin <= fabs(h) <= hmax required", GSL_EINVAL); + } + + if (hstart > 0.0 || hstart < 0.0) + { + d->h = hstart; + } + else + { + GSL_ERROR_NULL ("invalid hstart", GSL_EINVAL); + } + + return GSL_SUCCESS; +} void gsl_odeiv2_driver_free (gsl_odeiv2_driver * state) diff -up wrk/ode-initval2/evolve.c.wrk wrk/ode-initval2/evolve.c --- wrk/ode-initval2/evolve.c.wrk 2013-01-29 13:40:22.482035214 +0100 +++ wrk/ode-initval2/evolve.c 2013-01-29 13:07:32.000000000 +0100 @@ -138,16 +138,24 @@ gsl_odeiv2_evolve_apply (gsl_odeiv2_evol DBL_MEMCPY (e->y0, y, e->dimension); - /* Calculate initial dydt once if the method can benefit. */ - + /* Calculate initial dydt once or reuse previous value if the method + can benefit. */ + if (step->type->can_use_dydt_in) { - int status = GSL_ODEIV_FN_EVAL (dydt, t0, y, e->dydt_in); + if (e->count == 0) + { + int status = GSL_ODEIV_FN_EVAL (dydt, t0, y, e->dydt_in); - if (status) - { - return status; - } + if (status) + { + return status; + } + } + else + { + DBL_MEMCPY (e->dydt_in, e->dydt_out, e->dimension); + } } try_step: diff -up wrk/ode-initval2/gsl_odeiv2.h.wrk wrk/ode-initval2/gsl_odeiv2.h --- wrk/ode-initval2/gsl_odeiv2.h.wrk 2013-01-29 13:40:22.468035129 +0100 +++ wrk/ode-initval2/gsl_odeiv2.h 2013-01-29 13:07:52.000000000 +0100 @@ -326,6 +326,7 @@ int gsl_odeiv2_driver_apply_fixed_step ( const unsigned long int n, double y[]); int gsl_odeiv2_driver_reset (gsl_odeiv2_driver * d); +int gsl_odeiv2_driver_reset_hstart (gsl_odeiv2_driver * d, const double hstart); void gsl_odeiv2_driver_free (gsl_odeiv2_driver * state); __END_DECLS diff -up wrk/ode-initval2/Makefile.am.wrk wrk/ode-initval2/Makefile.am diff -up wrk/ode-initval2/Makefile.in.wrk wrk/ode-initval2/Makefile.in diff -up wrk/ode-initval2/modnewton1.c.wrk wrk/ode-initval2/modnewton1.c diff -up wrk/ode-initval2/msadams.c.wrk wrk/ode-initval2/msadams.c diff -up wrk/ode-initval2/msbdf.c.wrk wrk/ode-initval2/msbdf.c --- wrk/ode-initval2/msbdf.c.wrk 2013-01-29 13:40:22.473035159 +0100 +++ 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); } + state->abscorscaled = gsl_vector_alloc (dim); + + 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); state->driver = NULL; @@ -926,14 +953,14 @@ msbdf_corrector (void *vstate, const gsl } 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) { /* 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 /* Relative step length estimate for current order */ - 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); /* 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)); } ordp1est = 1.0 / (pow (bias2 * gsl_blas_dnrm2 (svec) / sqrt (dim) @@ -1020,8 +1047,6 @@ msbdf_eval_order (gsl_vector * abscor, g #endif } - *ordwait = *ord + 2; - return GSL_SUCCESS; } @@ -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 } } - /* 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; 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]); } } @@ -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/odeiv_util.h.wrk wrk/ode-initval2/odeiv_util.h diff -up wrk/ode-initval2/rk1imp.c.wrk wrk/ode-initval2/rk1imp.c diff -up wrk/ode-initval2/rk2.c.wrk wrk/ode-initval2/rk2.c diff -up wrk/ode-initval2/rk2imp.c.wrk wrk/ode-initval2/rk2imp.c diff -up wrk/ode-initval2/rk4.c.wrk wrk/ode-initval2/rk4.c diff -up wrk/ode-initval2/rk4imp.c.wrk wrk/ode-initval2/rk4imp.c --- wrk/ode-initval2/rk4imp.c.wrk 2013-01-29 13:40:22.479035196 +0100 +++ wrk/ode-initval2/rk4imp.c 2013-01-29 13:08:32.000000000 +0100 @@ -249,7 +249,7 @@ rk4imp_apply (void *vstate, size_t dim, gsl_matrix *A = state->A; const double b[] = { 0.5, 0.5 }; - const double c[] = { (3 - sqrt (3)) / 6, (3 + sqrt (3)) / 6 }; + const double c[] = { (3.0 - M_SQRT3) / 6.0, (3.0 + M_SQRT3) / 6.0 }; gsl_matrix_set (A, 0, 0, 1.0 / 4); gsl_matrix_set (A, 0, 1, (3 - 2 * sqrt (3)) / 12); diff -up wrk/ode-initval2/rk8pd.c.wrk wrk/ode-initval2/rk8pd.c diff -up wrk/ode-initval2/rkck.c.wrk wrk/ode-initval2/rkck.c diff -up wrk/ode-initval2/rkf45.c.wrk wrk/ode-initval2/rkf45.c --- wrk/ode-initval2/rkf45.c.wrk 2013-01-29 13:40:22.464035105 +0100 +++ wrk/ode-initval2/rkf45.c 2013-01-29 13:08:45.000000000 +0100 @@ -369,7 +369,7 @@ rkf45_free (void *vstate) static const gsl_odeiv2_step_type rkf45_type = { "rkf45", /* name */ 1, /* can use dydt_in */ - 0, /* gives exact dydt_out */ + 1, /* gives exact dydt_out */ &rkf45_alloc, &rkf45_apply, &stepper_set_driver_null, diff -up wrk/ode-initval2/rksubs.c.wrk wrk/ode-initval2/rksubs.c diff -up wrk/ode-initval2/step.c.wrk wrk/ode-initval2/step.c diff -up wrk/ode-initval2/step_utils.c.wrk wrk/ode-initval2/step_utils.c diff -up wrk/ode-initval2/test.c.wrk wrk/ode-initval2/test.c --- wrk/ode-initval2/test.c.wrk 2013-01-29 13:40:22.487035244 +0100 +++ 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; - + if (rhs_xsin_reset) { rhs_xsin_reset = 0; @@ -1073,6 +1073,7 @@ test_odeiv_stepper (const gsl_odeiv2_ste { 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)); 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); - + extern int nfe, nje; nfe = 0; nje = 0; - + while (t < t1) { 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 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 }; /* 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 */ steppers[0] = gsl_odeiv2_step_bsimp; @@ -2024,7 +2024,7 @@ test_extreme_problems (void) const double val1 = y[k]; 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) } 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); s = gsl_odeiv2_driver_apply (d, &t, t1, y); if (s != GSL_SUCCESS) { - 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); } 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 */ @@ -2493,14 +2525,12 @@ main (void) } /* Special tests */ - + test_nonstiff_problems (); test_stiff_problems (); test_extreme_problems (); - test_driver (); - exit (gsl_test_summary ()); } diff -up wrk/ode-initval2/TODO.wrk wrk/ode-initval2/TODO