rpms/gsl
3
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

temporary diagnostic

Browse Source
This commit is contained in:
Frantisek Kluknavsky 2013-01-29 13:45:30 +01:00
parent 9747dd763d
commit 38920d77a6
1 changed files with 176 additions and 150 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

326
wrk.patch
View File

@ -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@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
--- 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