328 lines
8.6 KiB
Diff
328 lines
8.6 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 15:26:46.105532490 +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;
|
|
@@ -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,
|
|
º1coeff, &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\n");
|
|
+//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 ());
|
|
}
|