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  <tuomo.keskitalo@iki.fi>

+

+	* 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  <rhys.ulerich@gmail.com>

+

+	* test.c: Correct two out-of-order declarations.

+	Thanks to Brian Gladman for spotting the problem.

+

+2012-03-10  Tuomo Keskitalo  <tuomo.keskitalo@iki.fi>

+

+	* 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  <tuomo.keskitalo@iki.fi>

+

+	* rkf45.c: Bug correction: Set gives_exact_dydt_out to 1 in

+	rkf45_type

+

+2012-01-14  Tuomo Keskitalo  <tuomo.keskitalo@iki.fi>

+

+	* 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  <bjg@network-theory.co.uk>

+

+	* 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 */