From 38920d77a606cddb8b6e247d1389f52082b1ddd4 Mon Sep 17 00:00:00 2001 From: Frantisek Kluknavsky Date: Tue, 29 Jan 2013 13:45:30 +0100 Subject: [PATCH] temporary diagnostic --- wrk.patch | 326 +++++++++++++++++++++++++++++------------------------- 1 file changed, 176 insertions(+), 150 deletions(-) diff --git a/wrk.patch b/wrk.patch index 445c1cc..9c4dd3d 100644 --- a/wrk.patch +++ b/wrk.patch @@ -1,167 +1,158 @@ +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 ---- wrk/ode-initval2/control.c.wrk 2013-01-23 17:06:22.402902397 +0100 -+++ wrk/ode-initval2/control.c 2013-01-23 17:36:43.945159375 +0100 -@@ -78,8 +78,11 @@ gsl_odeiv2_control_hadjust (gsl_odeiv2_c - const double y[], const double yerr[], - const double dydt[], double *h) - { -- return c->type->hadjust (c->state, s->dimension, s->type->order (s->state), -+// printf("pred hadjust\n%d\n", s->type->order (s->state)); -+ int ret = c->type->hadjust (c->state, s->dimension, s->type->order (s->state), - y, yerr, dydt, h); -+ // printf("po hadjust\n%d\n", s->type->order (s->state)); -+return ret; - } - - int +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 ---- wrk/ode-initval2/cstd.c.wrk 2013-01-18 16:02:19.566833381 +0100 -+++ wrk/ode-initval2/cstd.c 2013-01-23 17:37:01.010263153 +0100 -@@ -80,11 +80,14 @@ std_control_init (void *vstate, +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; } -+#define PD(fl) printf("%lu\n", *((unsigned long *)&(fl))) ++int ++gsl_odeiv2_driver_reset_hstart (gsl_odeiv2_driver * d, const double hstart) ++{ ++ /* Resets current driver and sets initial step size to hstart */ + - static int - std_control_hadjust (void *vstate, size_t dim, unsigned int ord, - const double y[], const double yerr[], const double yp[], - double *h) - { -+//printf("hadjustim\n"); - std_control_state_t *state = (std_control_state_t *) vstate; ++ 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; ++} - const double eps_abs = state->eps_abs; -@@ -109,9 +112,18 @@ std_control_hadjust (void *vstate, size_ - - if (rmax > 1.1) - { -+// printf("skracujem\n"); - /* decrease step, no more than factor of 5, but a fraction S more - than scaling suggests (for better accuracy) */ -- double r = S / pow (rmax, 1.0 / ord); -+ double pom = pow (rmax, 1.0 / ord); -+ double r = S / pom; -+// printf("ord %d", ord); -+//printf ("S\n"); -+//PD(S); -+//printf("rmax\n"); -+//PD(rmax); -+//printf("pom\n"); -+//PD(pom); - - if (r < 0.2) - r = 0.2; -@@ -122,6 +134,7 @@ std_control_hadjust (void *vstate, size_ - } - else if (rmax < 0.5) - { -+// printf("predlzujem\n"); - /* increase step, no more than factor of 5 */ - double r = S / pow (rmax, 1.0 / (ord + 1.0)); - -@@ -137,6 +150,7 @@ std_control_hadjust (void *vstate, size_ - } - else - { -+// printf("nemenim\n"); - /* no change */ - return GSL_ODEIV_HADJ_NIL; - } + 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-17 18:48:06.473062671 +0100 -+++ wrk/ode-initval2/evolve.c 2013-01-23 17:36:26.634053595 +0100 -@@ -28,6 +28,8 @@ +--- 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 - #include "odeiv_util.h" + DBL_MEMCPY (e->y0, y, e->dimension); -+#define PD(fl) printf("%lu\n", *((unsigned long *)&(fl))) -+ - gsl_odeiv2_evolve * - gsl_odeiv2_evolve_alloc (size_t dim) - { -@@ -118,6 +120,8 @@ gsl_odeiv2_evolve_apply (gsl_odeiv2_evol - const gsl_odeiv2_system * dydt, - double *t, double t1, double *h, double y[]) - { -+// printf("vstup do apply\n"); -+// PD(*h); - const double t0 = *t; - double h0 = *h; - int step_status; -@@ -151,7 +155,8 @@ gsl_odeiv2_evolve_apply (gsl_odeiv2_evol +- /* 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: -- -+//printf("po navasti\n"); -+//PD(h0); - if ((dt >= 0.0 && h0 > dt) || (dt < 0.0 && h0 < dt)) - { - h0 = dt; -@@ -164,9 +169,11 @@ try_step: - - if (step->type->can_use_dydt_in) - { -+ printf("pred stepapply vonku\n%d\n", step->type->order (step->state)); - step_status = - gsl_odeiv2_step_apply (step, t0, h0, y, e->yerr, e->dydt_in, - e->dydt_out, dydt); -+ printf("po stepapply vonku\n%d\n", step->type->order (step->state)); - } - else - { -@@ -206,6 +213,7 @@ try_step: - - DBL_MEMCPY (y, e->y0, dydt->dimension); - e->failed_steps++; -+ //printf("skacem prvy\n"); - goto try_step; - } - else -@@ -228,16 +236,19 @@ try_step: - { - *t = t0 + h0; - } -- -+//printf("pred con\n"); -+//PD(h0); - if (con != NULL) - { - /* Check error and attempt to adjust the step. */ - - double h_old = h0; -- -+// printf("pred hadjust vonku\n%d\n", step->type->order (step->state)); - const int hadjust_status - = - gsl_odeiv2_control_hadjust (con, step, y, e->yerr, e->dydt_out, &h0); -+//printf("po hadjust\n"); -+//PD(h0); - - if (hadjust_status == GSL_ODEIV_HADJ_DEC) - { -@@ -254,6 +265,8 @@ try_step: - - DBL_MEMCPY (y, e->y0, dydt->dimension); - e->failed_steps++; -+//printf("skacem druhy\n"); -+//PD(h0); - goto try_step; - } - else -@@ -279,6 +292,9 @@ try_step: - *h = h0; - } - -+// printf("vystup z apply\n"); -+// PD(*h); -+ - return step_status; - } +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-09 10:35:45.259960403 +0100 +--- 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 */ @@ -342,8 +333,42 @@ diff -up wrk/ode-initval2/msbdf.c.wrk wrk/ode-initval2/msbdf.c 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-09 10:48:22.051415928 +0100 +--- 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; @@ -704,3 +729,4 @@ diff -up wrk/ode-initval2/test.c.wrk wrk/ode-initval2/test.c - exit (gsl_test_summary ()); } +diff -up wrk/ode-initval2/TODO.wrk wrk/ode-initval2/TODO