gsl/wrk.patch
Frantisek Kluknavsky 9bc8e33e2c temporary diagnostic
2013-01-23 16:00:58 +01:00

337 lines
8.8 KiB
Diff

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-22 12:26:29.266471308 +0100
@@ -80,11 +80,14 @@ std_control_init (void *vstate,
return GSL_SUCCESS;
}
+#define PD(fl) printf("%lu\n", *((unsigned long *)&(fl)))
+
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));
@@ -137,6 +150,7 @@ std_control_hadjust (void *vstate, size_
}
else
{
+// printf("nemenim\n");
/* no change */
return GSL_ODEIV_HADJ_NIL;
}
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-22 10:52:40.895672072 +0100
@@ -28,6 +28,8 @@
#include "odeiv_util.h"
+#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:
-
+//printf("po navasti\n");
+//PD(h0);
if ((dt >= 0.0 && h0 > dt) || (dt < 0.0 && h0 < dt))
{
h0 = dt;
@@ -206,6 +211,7 @@ try_step:
DBL_MEMCPY (y, e->y0, dydt->dimension);
e->failed_steps++;
+// printf("skacem prvy\n");
goto try_step;
}
else
@@ -228,7 +234,8 @@ try_step:
{
*t = t0 + h0;
}
-
+//printf("pred con\n");
+//PD(h0);
if (con != NULL)
{
/* Check error and attempt to adjust the step. */
@@ -238,6 +245,8 @@ try_step:
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 +263,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 +290,9 @@ try_step:
*h = h0;
}
+// printf("vystup z apply\n");
+// PD(*h);
+
return step_status;
}
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 2013-01-23 16:00:47.836028674 +0100
@@ -43,6 +43,8 @@
framework.
*/
+#define PD(fl) printf("%lu\n", *((unsigned long *)&(fl)))
+
#include <config.h>
#include <stdlib.h>
#include <string.h>
@@ -116,7 +118,6 @@ static void *
msbdf_alloc (size_t dim)
{
msbdf_state_t *state = (msbdf_state_t *) malloc (sizeof (msbdf_state_t));
-
if (state == 0)
{
GSL_ERROR_NULL ("failed to allocate space for msbdf_state", GSL_ENOMEM);
@@ -611,6 +612,15 @@ msbdf_calccoeffs (const size_t ord, cons
printf ("errcoeff=%.5e\n", *errcoeff);
#endif
+/*PD(l[0]);
+PD(l[1]);
+PD(*errcoeff);
+PD(*ordm1coeff);
+PD(*ordp1coeff);
+PD(*ordp2coeff);
+PD(*gamma);
+printf("\n");
+*/
return GSL_SUCCESS;
}
@@ -1082,7 +1092,8 @@ msbdf_apply (void *vstate, size_t dim, d
const gsl_odeiv2_system * sys)
{
/* Carries out a step by BDF linear multistep methods. */
-
+//printf("zaciatok msbdf\n");
+//PD(h);
msbdf_state_t *state = (msbdf_state_t *) vstate;
double *const z = state->z;
@@ -1109,7 +1120,7 @@ msbdf_apply (void *vstate, size_t dim, d
const size_t max_failcount = 3;
size_t i;
-
+printf("ord at msbdf start\n%d\n", ord);
#ifdef DEBUG
{
size_t di;
@@ -1329,7 +1340,7 @@ msbdf_apply (void *vstate, size_t dim, d
{
const int deltaord = ord - ordprev[0];
-
+printf("ord\n%d\n%d\n", ord, ordprev[0]);
if (deltaord > 1 || deltaord < -1)
{
printf ("-- order change %d\n", deltaord);
@@ -1466,7 +1477,10 @@ msbdf_apply (void *vstate, size_t dim, d
/* Calculate polynomial coefficients (l), error coefficient and
auxiliary coefficients
*/
-
+/*printf("%d\n", ord);
+printf("%lu\n", state->ordwait);
+PD(h);
+PD(hprev[0]);*/
msbdf_calccoeffs (ord, state->ordwait, h, hprev, l, &errcoeff,
&ordm1coeff, &ordp1coeff, &ordp2coeff, &gamma);
@@ -1675,7 +1689,8 @@ msbdf_apply (void *vstate, size_t dim, d
printf ("-- nJ=%d, nM=%d\n", (int) state->nJ, (int) state->nM);
#endif
}
-
+printf("koniec msbdf, ord:\n%d\n", ord);
+//PD(h);
return GSL_SUCCESS;
}
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 2013-01-23 15:26:48.797549396 +0100
@@ -923,6 +923,18 @@ rhs_ringmod (double t, const double y[],
f[13] = (-y[0] + uin1 - (ri + rg1) * y[13]) / ls1;
f[14] = (-y[1] - (rc + rg1) * y[14]) / ls1;
+ /*int i;
+ printf("temp states\n");
+ for (i=0; i<15; ++i) {
+
+ printf("%lu\n", *((unsigned long *)(&(y[i]))));
+ }
+
+ printf("%lu\n", *((unsigned long *)(&(rg2))));
+ printf("%lu\n", *((unsigned long *)(&(ls2))));
+ printf("%lu\n", *((unsigned long *)(&(rg3))));
+ printf("%lu\n", *((unsigned long *)(&(ls3))));
+*/
return GSL_SUCCESS;
}
@@ -1239,6 +1251,8 @@ test_evolve_system (const gsl_odeiv2_ste
gsl_odeiv2_driver_free (d);
}
+#define PD(fl) printf("%lu\n", *((unsigned long *)&(fl)))
+
int
sys_driver (const gsl_odeiv2_step_type * T,
const gsl_odeiv2_system * sys,
@@ -1264,10 +1278,23 @@ sys_driver (const gsl_odeiv2_step_type *
nfe = 0;
nje = 0;
+ int poc=0;
while (t < t1)
{
+// printf("pred\n");
+// PD(h);
s = gsl_odeiv2_evolve_apply (d->e, d->c, d->s, sys, &t, t1, &h, y);
+ printf("\nkrok %d\n", ++poc);
+//printf("po\n");
+//PD(h);
+/* int k;
+ for (k=0; k<15; ++k) {
+ printf(" %lu\n",(*(unsigned long int *)(&(y[k]))));
+ }
+ printf(" %lu\n",(*(unsigned long int *)(&t)));
+*/
+ //printf("\n");
#ifdef DEBUG
printf ("%.5e %.5e %.5e %d\n", t, y[0], y[1],
gsl_odeiv2_step_order (d->s));
@@ -2003,10 +2030,11 @@ test_extreme_problems (void)
for (i = 0; steppers[i] != 0; i++)
{
+ // printf("spustam driver p=%ld i=%ld\n", p, i);
int s = sys_driver (steppers[i], prob[p], start[p], end[p],
initstepsize[p], &y[sd[p] * i],
epsabs[p], epsrel[p], probname[p]);
-
+ //printf("koniec drivera\n");
if (s != GSL_SUCCESS)
{
printf ("start=%.5e, initstepsize=%.5e\n", start[p],
@@ -2018,7 +2046,7 @@ test_extreme_problems (void)
/* Compare results */
- for (i = 1; steppers[i] != 0; i++)
+ for (i = 0; steppers[i] != 0; i++)
for (k = 0; k < sd[p]; k++)
{
const double val1 = y[k];
@@ -2450,7 +2478,7 @@ main (void)
/* Basic tests for all steppers */
- for (i = 0; p[i].type != 0; i++)
+ /*for (i = 0; p[i].type != 0; i++)
{
test_stepper (p[i].type);
}
@@ -2469,9 +2497,9 @@ main (void)
test_stepsize_fail (p[i].type, p[i].h);
test_user_break (p[i].type, p[i].h);
}
-
+*/
/* Derivative test for explicit steppers */
-
+/*
explicit_stepper[0].type = gsl_odeiv2_step_rk4;
explicit_stepper[0].h = 1.0e-3;
explicit_stepper[1].type = gsl_odeiv2_step_rk2;
@@ -2491,16 +2519,16 @@ main (void)
test_stepfn (explicit_stepper[i].type);
test_stepfn2 (explicit_stepper[i].type);
}
-
+*/
/* Special tests */
- test_nonstiff_problems ();
+ // test_nonstiff_problems ();
- test_stiff_problems ();
+ // test_stiff_problems ();
test_extreme_problems ();
- test_driver ();
+ //test_driver ();
exit (gsl_test_summary ());
}