[gsl/temp_odecrash] temporary diagnostic
Frantisek Kluknavsky
fkluknav at fedoraproject.org
Tue Jan 29 15:25:24 UTC 2013
commit 38920d77a606cddb8b6e247d1389f52082b1ddd4
Author: Frantisek Kluknavsky <fkluknav at redhat.com>
Date: Tue Jan 29 13:45:30 2013 +0100
temporary diagnostic
wrk.patch | 324 +++++++++++++++++++++++++++++++++----------------------------
1 files changed, 175 insertions(+), 149 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 <tuomo.keskitalo at 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 at 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 at 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 at iki.fi>
++
++ * rkf45.c: Bug correction: Set gives_exact_dydt_out to 1 in
++ rkf45_type
++
++2012-01-14 Tuomo Keskitalo <tuomo.keskitalo at 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 at 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
---- 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;
-
- 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));
++ 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;
++}
-@@ -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
- }
-
- try_step:
+- /* Calculate initial dydt once if the method can benefit. */
-
-+//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:
-
++ /* Calculate initial dydt once or reuse previous value if the method
++ can benefit. */
++
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;
+- 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);
++ }
}
-+// printf("vystup z apply\n");
-+// PD(*h);
-+
- return step_status;
- }
-
+ 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-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
More information about the scm-commits
mailing list