/* ode-initval2/msbdf.c
*
* Copyright (C) 2009, 2010 Tuomo Keskitalo
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 3 of the License, or (at
* your option) any later version.
*
* This program is distributed in the hope that it will be useful, but
* WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/
/* A variable-coefficient linear multistep backward differentiation
formula (BDF) method in Nordsieck form. This stepper uses the
explicit BDF formula as predictor and implicit BDF formula as
corrector. A modified Newton iteration method is used to
solve the system of non-linear equations. Method order varies
dynamically between 1 and 5.
References:
Byrne, G. D., and Hindmarsh, A. C., A Polyalgorithm for the
Numerical Solution of Ordinary Differential Equations,
ACM Trans. Math. Software, 1 (1975), pp. 71-96.
Brown, P. N., Byrne, G. D., and Hindmarsh, A. C., VODE: A
Variable-coefficient ODE Solver, SIAM J. Sci. Stat. Comput. 10,
(1989), pp. 1038-1051.
Hindmarsh, A. C., Brown, P. N., Grant, K. E., Lee, S. L., Serban,
R., Shumaker, D. E., and Woodward, C. S., SUNDIALS: Suite of
Nonlinear and Differential/Algebraic Equation Solvers, ACM
Trans. Math. Software 31 (2005), pp. 363-396.
Note: The algorithms have been adapted for GSL ode-initval2
framework.
*/
#include <config.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_odeiv2.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include "odeiv_util.h"
/* Maximum order of BDF methods */
#define MSBDF_MAX_ORD 5
/* Steps until Jacobian evaluation is forced */
#define MSBDF_JAC_WAIT 50
/* Steps until iteration matrix M evaluation is forced */
#define MSBDF_M_WAIT 20
typedef struct
{
/* Nordsieck history matrix. Includes concatenated
Nordsieck vectors [y_n, h*y_n', (h^2/2!)*y_n'', ...,
(h^ord/ord!)*d^(ord)(y_n)]. Nordsieck vector number i is located
at z[i*dim] (i=0..ord).
*/
double *z;
double *zbackup; /* backup of Nordsieck matrix */
double *ytmp; /* work area */
double *ytmp2; /* work area */
double *l; /* polynomial coefficients */
double *hprev; /* previous step sizes */
double *hprevbackup; /* backup of hprev */
size_t *ordprev; /* orders of previous calls */
size_t *ordprevbackup; /* backup of ordprev */
double *errlev; /* desired error level of y */
gsl_vector *abscor; /* absolute y values for correction */
gsl_vector *abscorscaled; /* scaled abscor for order evaluation */
gsl_vector *relcor; /* relative y values for correction */
gsl_vector *svec; /* saved abscor & work area */
gsl_vector *tempvec; /* work area */
const gsl_odeiv2_driver *driver; /* pointer to gsl_odeiv2_driver object */
gsl_matrix *dfdy; /* Jacobian */
double *dfdt; /* storage for time derivative of f */
gsl_matrix *M; /* Newton iteration matrix */
gsl_permutation *p; /* permutation for LU decomposition of M */
gsl_vector *rhs; /* right hand side equations (-G) */
long int ni; /* stepper call counter */
size_t ord; /* current order of method */
double tprev; /* t point of previous call */
size_t ordwait; /* counter for order change */
size_t ordwaitbackup; /* backup of ordwait */
size_t failord; /* order of convergence failure */
double failt; /* t point of convergence failure */
double ordp1coeffprev; /* saved order coefficient */
size_t nJ; /* step counter for Jacobian evaluation */
size_t nM; /* step counter for update of M */
double gammaprev; /* gamma of previous call */
double gammaprevbackup; /* backup of gammaprev */
size_t failcount; /* counter for rejected steps */
}
msbdf_state_t;
/* Introduce msbdf_reset for use in msbdf_alloc and _apply */
static int msbdf_reset (void *, size_t);
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);
}
state->z = (double *) malloc ((MSBDF_MAX_ORD + 1) * dim * sizeof (double));
if (state->z == 0)
{
free (state);
GSL_ERROR_NULL ("failed to allocate space for z", GSL_ENOMEM);
}
state->zbackup =
(double *) malloc ((MSBDF_MAX_ORD + 1) * dim * sizeof (double));
if (state->zbackup == 0)
{
free (state->z);
free (state);
GSL_ERROR_NULL ("failed to allocate space for zbackup", GSL_ENOMEM);
}
state->ytmp = (double *) malloc (dim * sizeof (double));
if (state->ytmp == 0)
{
free (state->zbackup);
free (state->z);
free (state);
GSL_ERROR_NULL ("failed to allocate space for ytmp", GSL_ENOMEM);
}
state->ytmp2 = (double *) malloc (dim * sizeof (double));
if (state->ytmp2 == 0)
{
free (state->ytmp);
free (state->zbackup);
free (state->z);
free (state);
GSL_ERROR_NULL ("failed to allocate space for ytmp2", GSL_ENOMEM);
}
state->l = (double *) malloc ((MSBDF_MAX_ORD + 1) * sizeof (double));
if (state->l == 0)
{
free (state->ytmp);
free (state->zbackup);
free (state->z);
free (state);
GSL_ERROR_NULL ("failed to allocate space for l", GSL_ENOMEM);
}
state->hprev = (double *) malloc (MSBDF_MAX_ORD * sizeof (double));
if (state->hprev == 0)
{
free (state->l);
free (state->ytmp2);
free (state->ytmp);
free (state->zbackup);
free (state->z);
free (state);
GSL_ERROR_NULL ("failed to allocate space for hprev", GSL_ENOMEM);
}
state->hprevbackup = (double *) malloc (MSBDF_MAX_ORD * sizeof (double));
if (state->hprevbackup == 0)
{
free (state->hprev);
free (state->l);
free (state->ytmp2);
free (state->ytmp);
free (state->zbackup);
free (state->z);
free (state);
GSL_ERROR_NULL ("failed to allocate space for hprevbackup", GSL_ENOMEM);
}
state->ordprev = (size_t *) malloc (MSBDF_MAX_ORD * sizeof (size_t));
if (state->ordprev == 0)
{
free (state->hprevbackup);
free (state->hprev);
free (state->l);
free (state->ytmp2);
free (state->ytmp);
free (state->zbackup);
free (state->z);
free (state);
GSL_ERROR_NULL ("failed to allocate space for ordprev", GSL_ENOMEM);
}
state->ordprevbackup = (size_t *) malloc (MSBDF_MAX_ORD * sizeof (size_t));
if (state->ordprevbackup == 0)
{
free (state->ordprev);
free (state->hprevbackup);
free (state->hprev);
free (state->l);
free (state->ytmp2);
free (state->ytmp);
free (state->zbackup);
free (state->z);
free (state);
GSL_ERROR_NULL ("failed to allocate space for ordprevbackup",
GSL_ENOMEM);
}
state->errlev = (double *) malloc (dim * sizeof (double));
if (state->errlev == 0)
{
free (state->ordprevbackup);
free (state->ordprev);
free (state->hprevbackup);
free (state->hprev);
free (state->l);
free (state->ytmp2);
free (state->ytmp);
free (state->zbackup);
free (state->z);
free (state);
GSL_ERROR_NULL ("failed to allocate space for errlev", GSL_ENOMEM);
}
state->abscor = gsl_vector_alloc (dim);
if (state->abscor == 0)
{
free (state->errlev);
free (state->ordprevbackup);
free (state->ordprev);
free (state->hprevbackup);
free (state->hprev);
free (state->l);
free (state->ytmp2);
free (state->ytmp);
free (state->zbackup);
free (state->z);
free (state);
GSL_ERROR_NULL ("failed to allocate space for abscor", GSL_ENOMEM);
}
state->relcor = gsl_vector_alloc (dim);
if (state->relcor == 0)
{
gsl_vector_free (state->abscor);
free (state->errlev);
free (state->ordprevbackup);
free (state->ordprev);
free (state->hprevbackup);
free (state->hprev);
free (state->l);
free (state->ytmp2);
free (state->ytmp);
free (state->zbackup);
free (state->z);
free (state);
GSL_ERROR_NULL ("failed to allocate space for relcor", GSL_ENOMEM);
}
state->svec = gsl_vector_alloc (dim);
if (state->svec == 0)
{
gsl_vector_free (state->relcor);
gsl_vector_free (state->abscor);
free (state->errlev);
free (state->ordprevbackup);
free (state->ordprev);
free (state->hprevbackup);
free (state->hprev);
free (state->l);
free (state->ytmp2);
free (state->ytmp);
free (state->zbackup);
free (state->z);
free (state);
GSL_ERROR_NULL ("failed to allocate space for svec", GSL_ENOMEM);
}
state->tempvec = gsl_vector_alloc (dim);
if (state->tempvec == 0)
{
gsl_vector_free (state->svec);
gsl_vector_free (state->relcor);
gsl_vector_free (state->abscor);
free (state->errlev);
free (state->ordprevbackup);
free (state->ordprev);
free (state->hprevbackup);
free (state->hprev);
free (state->l);
free (state->ytmp2);
free (state->ytmp);
free (state->zbackup);
free (state->z);
free (state);
GSL_ERROR_NULL ("failed to allocate space for tempvec", GSL_ENOMEM);
}
state->dfdy = gsl_matrix_alloc (dim, dim);
if (state->dfdy == 0)
{
gsl_vector_free (state->tempvec);
gsl_vector_free (state->svec);
gsl_vector_free (state->relcor);
gsl_vector_free (state->abscor);
free (state->errlev);
free (state->ordprevbackup);
free (state->ordprev);
free (state->hprevbackup);
free (state->hprev);
free (state->l);
free (state->ytmp2);
free (state->ytmp);
free (state->zbackup);
free (state->z);
free (state);
GSL_ERROR_NULL ("failed to allocate space for dfdy", GSL_ENOMEM);
}
state->dfdt = (double *) malloc (dim * sizeof (double));
if (state->dfdt == 0)
{
gsl_matrix_free (state->dfdy);
gsl_vector_free (state->tempvec);
gsl_vector_free (state->svec);
gsl_vector_free (state->relcor);
gsl_vector_free (state->abscor);
free (state->errlev);
free (state->ordprevbackup);
free (state->ordprev);
free (state->hprevbackup);
free (state->hprev);
free (state->l);
free (state->ytmp2);
free (state->ytmp);
free (state->zbackup);
free (state->z);
free (state);
GSL_ERROR_NULL ("failed to allocate space for dfdt", GSL_ENOMEM);
}
state->M = gsl_matrix_alloc (dim, dim);
if (state->M == 0)
{
free (state->dfdt);
gsl_matrix_free (state->dfdy);
gsl_vector_free (state->tempvec);
gsl_vector_free (state->svec);
gsl_vector_free (state->relcor);
gsl_vector_free (state->abscor);
free (state->errlev);
free (state->ordprevbackup);
free (state->ordprev);
free (state->hprevbackup);
free (state->hprev);
free (state->l);
free (state->ytmp2);
free (state->ytmp);
free (state->zbackup);
free (state->z);
free (state);
GSL_ERROR_NULL ("failed to allocate space for M", GSL_ENOMEM);
}
state->p = gsl_permutation_alloc (dim);
if (state->p == 0)
{
gsl_matrix_free (state->M);
free (state->dfdt);
gsl_matrix_free (state->dfdy);
gsl_vector_free (state->tempvec);
gsl_vector_free (state->svec);
gsl_vector_free (state->relcor);
gsl_vector_free (state->abscor);
free (state->errlev);
free (state->ordprevbackup);
free (state->ordprev);
free (state->hprevbackup);
free (state->hprev);
free (state->l);
free (state->ytmp2);
free (state->ytmp);
free (state->zbackup);
free (state->z);
free (state);
GSL_ERROR_NULL ("failed to allocate space for p", GSL_ENOMEM);
}
state->rhs = gsl_vector_alloc (dim);
if (state->rhs == 0)
{
gsl_permutation_free (state->p);
gsl_matrix_free (state->M);
free (state->dfdt);
gsl_matrix_free (state->dfdy);
gsl_vector_free (state->tempvec);
gsl_vector_free (state->svec);
gsl_vector_free (state->relcor);
gsl_vector_free (state->abscor);
free (state->errlev);
free (state->ordprevbackup);
free (state->ordprev);
free (state->hprevbackup);
free (state->hprev);
free (state->l);
free (state->ytmp2);
free (state->ytmp);
free (state->zbackup);
free (state->z);
free (state);
GSL_ERROR_NULL ("failed to allocate space for rhs", GSL_ENOMEM);
}
state->abscorscaled = gsl_vector_alloc (dim);
if (state->abscorscaled == 0)
{
gsl_vector_free (state->rhs);
gsl_permutation_free (state->p);
gsl_matrix_free (state->M);
free (state->dfdt);
gsl_matrix_free (state->dfdy);
gsl_vector_free (state->tempvec);
gsl_vector_free (state->svec);
gsl_vector_free (state->relcor);
gsl_vector_free (state->abscor);
free (state->errlev);
free (state->ordprevbackup);
free (state->ordprev);
free (state->hprevbackup);
free (state->hprev);
free (state->l);
free (state->ytmp2);
free (state->ytmp);
free (state->zbackup);
free (state->z);
free (state);
GSL_ERROR_NULL ("failed to allocate space for abscorscaled", GSL_ENOMEM);
}
msbdf_reset ((void *) state, dim);
state->driver = NULL;
return state;
}
static int
msbdf_failurehandler (void *vstate, const size_t dim, const double t)
{
/* Internal failure handler routine for msbdf. Adjusted strategy
for GSL: Decrease order if this is the second time a failure
has occurred at this order and point.
*/
msbdf_state_t *state = (msbdf_state_t *) vstate;
const size_t ord = state->ord;
if (ord > 1 && (ord - state->ordprev[0] == 0) &&
ord == state->failord && t == state->failt)
{
state->ord--;
}
/* Save information about failure */
state->failord = ord;
state->failt = t;
state->ni++;
/* Force reinitialization if failure took place at lowest
order
*/
if (ord == 1)
{
msbdf_reset (vstate, dim);
}
return GSL_SUCCESS;
}
static int
msbdf_calccoeffs (const size_t ord, const size_t ordwait,
const double h, const double hprev[],
double l[],
double *errcoeff, double *ordm1coeff,
double *ordp1coeff, double *ordp2coeff, double *gamma)
{
/* Calculates coefficients (l) of polynomial Lambda, error and
auxiliary order change evaluation coefficients.
*/
if (ord == 1)
{
l[0] = 1.0;
l[1] = 1.0;
*errcoeff = 0.5;
*ordp1coeff = 2.0;
{
const double hsum = h + hprev[0];
const double a5 = -1.5;
const double a6 = -1.0 - h / hsum;
const double c2 = 2.0 / (1.0 - a6 + a5);
*ordp2coeff = fabs (c2 * (h / hsum) * 3.0 * a5);
}
}
else
{
size_t i, j;
double hsum = h;
double coeff1 = -1.0;
double x;
/* Calculate the actual polynomial coefficients (l) */
DBL_ZERO_MEMSET (l, MSBDF_MAX_ORD + 1);
l[0] = 1.0;
l[1] = 1.0;
for (i = 2; i < ord; i++)
{
hsum += hprev[i - 2];
coeff1 += -1.0 / i;
for (j = i; j > 0; j--)
{
l[j] += h / hsum * l[j - 1];
}
}
coeff1 += -1.0 / ord;
x = -l[1] - coeff1;
for (i = ord; i > 0; i--)
{
l[i] += l[i - 1] * x;
}
#ifdef DEBUG
{
size_t di;
printf ("-- calccoeffs l: ");
for (di = 0; di < ord + 1; di++)
{
printf ("%.5e ", l[di]);
}
printf ("\n");
}
#endif
hsum += hprev[ord - 2];
{
const double coeff2 = -l[1] - h / hsum;
const double a1 = 1.0 - coeff2 + coeff1;
const double a2 = 1.0 + ord * a1;
/* Calculate error coefficient */
*errcoeff = fabs (a1 / (coeff1 * a2));
/* Calculate auxiliary coefficients used in evaluation of change
of order
*/
if (ordwait < 2)
{
const double a3 = coeff1 + 1.0 / ord;
const double a4 = coeff2 + h / hsum;
const double c1 = a3 / (1.0 - a4 + a3);
*ordm1coeff = fabs (c1 / (x / l[ord]));
*ordp1coeff = fabs (a2 / (l[ord] * (h / hsum) / x));
hsum += hprev[ord - 1];
{
const double a5 = coeff1 - 1.0 / (ord + 1.0);
const double a6 = coeff2 - h / hsum;
const double c2 = a2 / (1.0 - a6 + a5);
*ordp2coeff = fabs (c2 * (h / hsum) * (ord + 2) * a5);
}
}
}
}
*gamma = h / l[1];
#ifdef DEBUG
printf ("-- calccoeffs ordm1coeff=%.5e ", *ordm1coeff);
printf ("ordp1coeff=%.5e ", *ordp1coeff);
printf ("ordp2coeff=%.5e ", *ordp2coeff);
printf ("errcoeff=%.5e\n", *errcoeff);
#endif
return GSL_SUCCESS;
}
static int
msbdf_update (void *vstate, const size_t dim, gsl_matrix * dfdy, double *dfdt,
const double t, const double *y, const gsl_odeiv2_system * sys,
gsl_matrix * M, gsl_permutation * p,
const size_t iter, size_t * nJ, size_t * nM,
const double tprev, const double failt,
const double gamma, const double gammaprev, const double hratio)
{
/* Evaluates Jacobian dfdy and updates iteration matrix M
if criteria for update is met.
*/
/* Jacobian is evaluated
- at first step
- if MSBDF_JAC_WAIT steps have been made without re-evaluation
- in case of a convergence failure if
--- change in gamma is small, or
--- convergence failure resulted in step size decrease
*/
const double c = 0.2;
const double gammarel = fabs (gamma / gammaprev - 1.0);
if (*nJ == 0 || *nJ > MSBDF_JAC_WAIT ||
(t == failt && (gammarel < c || hratio < 1.0)))
{
#ifdef DEBUG
printf ("-- evaluate jacobian\n");
#endif
int s = GSL_ODEIV_JA_EVAL (sys, t, y, dfdy->data, dfdt);
if (s == GSL_EBADFUNC)
{
return s;
}
if (s != GSL_SUCCESS)
{
msbdf_failurehandler (vstate, dim, t);
#ifdef DEBUG
printf ("-- FAIL at jacobian function evaluation\n");
#endif
return s;
}
/* Reset counter */
*nJ = 0;
}
/* Iteration matrix M (and it's LU decomposition) is generated
- at first step
- if MSBDF_M_WAIT steps have been made without an update
- if change in gamma is significant (e.g. change in step size)
- if previous step was rejected
*/
if (*nM == 0 || *nM > MSBDF_M_WAIT || gammarel >= c ||
t == tprev || t == failt)
{
#ifdef DEBUG
printf ("-- update M, gamma=%.5e\n", gamma);
#endif
size_t i;
gsl_matrix_memcpy (M, dfdy);
gsl_matrix_scale (M, -gamma);
for (i = 0; i < dim; i++)
{
gsl_matrix_set (M, i, i, gsl_matrix_get (M, i, i) + 1.0);
}
{
int signum;
int s = gsl_linalg_LU_decomp (M, p, &signum);
if (s != GSL_SUCCESS)
{
return GSL_FAILURE;
}
}
/* Reset counter */
*nM = 0;
}
return GSL_SUCCESS;
}
static int
msbdf_corrector (void *vstate, const gsl_odeiv2_system * sys,
const double t, const double h, const size_t dim,
const double z[], const double errlev[],
const double l[], const double errcoeff,
gsl_vector * abscor, gsl_vector * relcor,
double ytmp[], double ytmp2[],
gsl_matrix * dfdy, double dfdt[], gsl_matrix * M,
gsl_permutation * p, gsl_vector * rhs,
size_t * nJ, size_t * nM,
const double tprev, const double failt,
const double gamma, const double gammaprev,
const double hprev0)
{
/* Calculates the correction step (abscor). Equation
system M = I - gamma * dfdy = -G is solved by Newton iteration.
*/
size_t mi, i;
const size_t max_iter = 3; /* Maximum number of iterations */
double convrate = 1.0; /* convergence rate */
double stepnorm = 0.0; /* norm of correction step */
double stepnormprev = 0.0; /* previous norm value */
/* Evaluate at predicted values */
{
int s = GSL_ODEIV_FN_EVAL (sys, t + h, z, ytmp);
if (s == GSL_EBADFUNC)
{
return s;
}
if (s != GSL_SUCCESS)
{
msbdf_failurehandler (vstate, dim, t);
#ifdef DEBUG
printf ("-- FAIL at user function evaluation\n");
#endif
return s;
}
}
/* Calculate correction step (abscor) */
gsl_vector_set_zero (abscor);
for (mi = 0; mi < max_iter; mi++)
{
const double safety = 0.3;
const double safety2 = 0.1;
/* Generate or update Jacobian and/or iteration matrix M if needed */
if (mi == 0)
{
int s = msbdf_update (vstate, dim, dfdy, dfdt, t + h, z,
sys, M, p, mi,
nJ, nM, tprev, failt,
gamma, gammaprev,
h / hprev0);
if (s != GSL_SUCCESS)
{
return s;
}
}
/* Evaluate the right hand side (-G) */
for (i = 0; i < dim; i++)
{
const double r = -1.0 * gsl_vector_get (abscor, i) -
z[1 * dim + i] / l[1] + gamma * ytmp[i];
gsl_vector_set (rhs, i, r);
}
/* Solve system of equations */
{
int s = gsl_linalg_LU_solve (M, p, rhs, relcor);
if (s != GSL_SUCCESS)
{
msbdf_failurehandler (vstate, dim, t);
#ifdef DEBUG
printf ("-- FAIL at LU_solve\n");
#endif
return GSL_FAILURE;
}
}
#ifdef DEBUG
{
size_t di;
printf ("-- dstep: ");
for (di = 0; di < dim; di++)
{
printf ("%.5e ", gsl_vector_get (relcor, di));
}
printf ("\n");
}
#endif
/* Add iteration results */
for (i = 0; i < dim; i++)
{
const double r =
gsl_vector_get (abscor, i) + gsl_vector_get (relcor, i);
gsl_vector_set (abscor, i, r);
ytmp2[i] = z[i] + r;
gsl_vector_set (relcor, i, gsl_vector_get (relcor, i) / errlev[i]);
}
#ifdef DEBUG
{
size_t di;
printf ("-- abscor: ");
for (di = 0; di < dim; di++)
{
printf ("%.5e ", gsl_vector_get (abscor, di));
}
printf ("\n");
}
#endif
/* Convergence test. Norms used are root-mean-square norms. */
stepnorm = gsl_blas_dnrm2 (relcor) / sqrt (dim);
if (mi > 0)
{
convrate = GSL_MAX_DBL (safety * convrate, stepnorm / stepnormprev);
}
else
{
convrate = 1.0;
}
{
const double convtest =
GSL_MIN_DBL (convrate, 1.0) * stepnorm * errcoeff / safety2;
#ifdef DEBUG
printf
("-- newt iter loop %d, errcoeff=%.5e, stepnorm =%.5e, convrate = %.5e, convtest = %.5e\n",
(int) mi, errcoeff, stepnorm, convrate, convtest);
#endif
if (convtest <= 1.0)
{
break;
}
}
/* Check for divergence during iteration */
{
const double div_const = 2.0;
if (mi > 1 && stepnorm > div_const * stepnormprev)
{
msbdf_failurehandler (vstate, dim, t);
#ifdef DEBUG
printf ("-- FAIL, diverging Newton iteration\n");
#endif
return GSL_FAILURE;
}
}
/* Evaluate at new y */
{
int s = GSL_ODEIV_FN_EVAL (sys, t + h, ytmp2, ytmp);
if (s == GSL_EBADFUNC)
{
return s;
}
if (s != GSL_SUCCESS)
{
msbdf_failurehandler (vstate, dim, t);
#ifdef DEBUG
printf ("-- FAIL at user function evaluation\n");
#endif
return s;
}
}
stepnormprev = stepnorm;
}
#ifdef DEBUG
printf ("-- Newton iteration exit at mi=%d\n", (int) mi);
#endif
/* Handle convergence failure */
if (mi == max_iter)
{
msbdf_failurehandler (vstate, dim, t);
#ifdef DEBUG
printf ("-- FAIL, max_iter reached\n");
#endif
return GSL_FAILURE;
}
return GSL_SUCCESS;
}
static int
msbdf_eval_order (gsl_vector * abscorscaled, gsl_vector * tempvec,
gsl_vector * svec, const double errcoeff,
const size_t dim, const double errlev[],
const double ordm1coeff, const double ordp1coeff,
const double ordp1coeffprev, const double ordp2coeff,
const double hprev[],
const double h, const double z[],
size_t * ord)
{
/* Evaluates and executes change in method order (current, current-1
or current+1). Order which maximizes the step length is selected.
*/
size_t i;
/* step size estimates at current order, order-1 and order+1 */
double ordest = 0.0;
double ordm1est = 0.0;
double ordp1est = 0.0;
const double safety = 1e-6;
const double bias = 6.0;
const double bias2 = 10.0;
const double min_incr = 1.5;
/* Relative step length estimate for current order */
ordest = 1.0 / (pow (bias * gsl_blas_dnrm2 (abscorscaled) / sqrt (dim)
* errcoeff, 1.0 / (*ord + 1)) + safety);
/* Relative step length estimate for order ord - 1 */
if (*ord > 1)
{
for (i = 0; i < dim; i++)
{
gsl_vector_set (tempvec, i, z[*ord * dim + i] / errlev[i]);
}
ordm1est = 1.0 / (pow (bias * gsl_blas_dnrm2 (tempvec) / sqrt (dim)
/ ordm1coeff, 1.0 / (*ord)) + safety);
}
else
{
ordm1est = 0.0;
}
/* Relative step length estimate for order ord + 1 */
if (*ord < MSBDF_MAX_ORD)
{
const double c = -ordp1coeff / ordp1coeffprev *
pow (h / hprev[1], *ord + 1);
for (i = 0; i < dim; i++)
{
gsl_vector_set (svec, i, gsl_vector_get (svec, i) * c +
gsl_vector_get (abscorscaled, i));
}
ordp1est = 1.0 / (pow (bias2 * gsl_blas_dnrm2 (svec) / sqrt (dim)
/ ordp2coeff, 1.0 / (*ord + 2)) + safety);
}
else
{
ordp1est = 0.0;
}
#ifdef DEBUG
printf
("-- eval_order ord=%d, ordest=%.5e, ordm1est=%.5e, ordp1est=%.5e\n",
(int) *ord, ordest, ordm1est, ordp1est);
#endif
/* Choose order that maximises step size and increases step
size markedly compared to current step
*/
if (ordm1est > ordest && ordm1est > ordp1est && ordm1est > min_incr)
{
*ord -= 1;
#ifdef DEBUG
printf ("-- eval_order order DECREASED to %d\n", (int) *ord);
#endif
}
else if (ordp1est > ordest && ordp1est > ordm1est && ordp1est > min_incr)
{
*ord += 1;
#ifdef DEBUG
printf ("-- eval_order order INCREASED to %d\n", (int) *ord);
#endif
}
return GSL_SUCCESS;
}
static int
msbdf_check_no_order_decrease (size_t const ordprev[])
{
/* Checks if order has not been decreased according to order history
array. Used in stability enhancement.
*/
size_t i;
for (i = 0; i < MSBDF_MAX_ORD - 1; i++)
{
if (ordprev[i + 1] > ordprev[i])
{
return 0;
}
}
return 1;
}
static int
msbdf_check_step_size_decrease (double const hprev[])
{
/* Checks if step size has decreased markedly according to
step size history array. Used in stability enhancement.
*/
size_t i;
double max = fabs (hprev[0]);
const double min = fabs (hprev[0]);
const double decrease_limit = 0.5;
for (i = 1; i < MSBDF_MAX_ORD; i++)
{
const double h = fabs (hprev[i]);
if (h > min && h > max)
{
max = h;
}
}
if (min / max < decrease_limit)
{
return 1;
}
return 0;
}
static int
msbdf_apply (void *vstate, size_t dim, double t, double h,
double y[], double yerr[],
const double dydt_in[], double dydt_out[],
const gsl_odeiv2_system * sys)
{
/* Carries out a step by BDF linear multistep methods. */
msbdf_state_t *state = (msbdf_state_t *) vstate;
double *const z = state->z;
double *const zbackup = state->zbackup;
double *const ytmp = state->ytmp;
double *const ytmp2 = state->ytmp2;
double *const l = state->l;
double *const hprev = state->hprev;
double *const hprevbackup = state->hprevbackup;
size_t *const ordprev = state->ordprev;
size_t *const ordprevbackup = state->ordprevbackup;
double *const errlev = state->errlev;
gsl_vector *const abscor = state->abscor;
gsl_vector *const abscorscaled = state->abscorscaled;
gsl_vector *const relcor = state->relcor;
gsl_vector *const svec = state->svec;
gsl_vector *const tempvec = state->tempvec;
size_t ord = state->ord; /* order for this step */
double ordm1coeff = 0.0;
double ordp1coeff = 0.0;
double ordp2coeff = 0.0;
double errcoeff = 0.0; /* error coefficient */
double gamma = 0.0; /* gamma coefficient */
const size_t max_failcount = 3;
size_t i;
#ifdef DEBUG
{
size_t di;
printf ("msbdf_apply: t=%.5e, ord=%d, h=%.5e, y:", t, (int) ord, h);
for (di = 0; di < dim; di++)
{
printf ("%.5e ", y[di]);
}
printf ("\n");
}
#endif
/* Check if t is the same as on previous stepper call (or last
failed call). This means that calculation of previous step failed
or the step was rejected, and therefore previous state will be
restored or the method will be reset.
*/
if (state->ni > 0 && (t == state->tprev || t == state->failt))
{
if (state->ni == 1)
{
/* No step has been accepted yet, reset method */
msbdf_reset (vstate, dim);
#ifdef DEBUG
printf ("-- first step was REJECTED, msbdf_reset called\n");
#endif
}
else
{
/* A succesful step has been saved, restore previous state. */
/* If previous step suggests order increase, but the step was
rejected, then do not increase order.
*/
if (ord > ordprev[0])
{
state->ord = ordprev[0];
ord = state->ord;
}
/* Restore previous state */
DBL_MEMCPY (z, zbackup, (MSBDF_MAX_ORD + 1) * dim);
DBL_MEMCPY (hprev, hprevbackup, MSBDF_MAX_ORD);
for (i = 0; i < MSBDF_MAX_ORD; i++)
{
ordprev[i] = ordprevbackup[i];
}
state->ordwait = state->ordwaitbackup;
state->gammaprev = state->gammaprevbackup;
#ifdef DEBUG
printf ("-- previous step was REJECTED, state restored\n");
#endif
}
/* If step is repeatedly rejected, then reset method */
state->failcount++;
if (state->failcount > max_failcount && state->ni > 1)
{
msbdf_reset (vstate, dim);
ord = state->ord;
#ifdef DEBUG
printf ("-- max_failcount reached, msbdf_reset called\n");
#endif
}
}
else
{
/* The previous step was accepted. Backup current state. */
DBL_MEMCPY (zbackup, z, (MSBDF_MAX_ORD + 1) * dim);
DBL_MEMCPY (hprevbackup, hprev, MSBDF_MAX_ORD);
for (i = 0; i < MSBDF_MAX_ORD; i++)
{
ordprevbackup[i] = ordprev[i];
}
state->ordwaitbackup = state->ordwait;
state->gammaprevbackup = state->gammaprev;
state->failcount = 0;
#ifdef DEBUG
if (state->ni > 0)
{
printf ("-- previous step was ACCEPTED, state saved\n");
}
#endif
}
#ifdef DEBUG
printf ("-- ord=%d, ni=%ld, ordwait=%d\n", (int) ord, state->ni,
(int) state->ordwait);
size_t di;
printf ("-- ordprev: ");
for (di = 0; di < MSBDF_MAX_ORD; di++)
{
printf ("%d ", (int) ordprev[di]);
}
printf ("\n");
#endif
/* Get desired error levels via gsl_odeiv2_control object through driver
object, which is a requirement for this stepper.
*/
if (state->driver == NULL)
{
return GSL_EFAULT;
}
else
{
size_t i;
for (i = 0; i < dim; i++)
{
if (dydt_in != NULL)
{
gsl_odeiv2_control_errlevel (state->driver->c, y[i],
dydt_in[i], h, i, &errlev[i]);
}
else
{
gsl_odeiv2_control_errlevel (state->driver->c, y[i],
0.0, h, i, &errlev[i]);
}
}
}
#ifdef DEBUG
{
size_t di;
printf ("-- errlev: ");
for (di = 0; di < dim; di++)
{
printf ("%.5e ", errlev[di]);
}
printf ("\n");
}
#endif
/* On first call initialize Nordsieck matrix */
if (state->ni == 0)
{
size_t i;
DBL_ZERO_MEMSET (z, (MSBDF_MAX_ORD + 1) * dim);
if (dydt_in != NULL)
{
DBL_MEMCPY (ytmp, dydt_in, dim);
}
else
{
int s = GSL_ODEIV_FN_EVAL (sys, t, y, ytmp);
if (s != GSL_SUCCESS)
{
return s;
}
}
DBL_MEMCPY (&z[0 * dim], y, dim);
DBL_MEMCPY (&z[1 * dim], ytmp, dim);
for (i = 0; i < dim; i++)
{
z[1 * dim + i] *= h;
}
}
/* Stability enhancement heuristic for msbdf: If order > 1 and order
has not been changed, check for decrease in step size, that is
not accompanied by a decrease in method order. This condition may
be indication of BDF method stability problems, a change in ODE
system, or convergence problems in Newton iteration. In all
cases, the strategy is to decrease method order.
*/
#ifdef DEBUG
printf ("-- check_no_order_decrease %d, check_step_size_decrease %d\n",
msbdf_check_no_order_decrease (ordprev),
msbdf_check_step_size_decrease (hprev));
#endif
if (ord > 1 &&
ord - ordprev[0] == 0 &&
msbdf_check_no_order_decrease (ordprev) &&
msbdf_check_step_size_decrease (hprev))
{
state->ord--;
state->ordwait = ord + 2;
ord = state->ord;
#ifdef DEBUG
printf ("-- stability enhancement decreased order to %d\n", (int) ord);
#endif
}
/* Sanity check */
{
const int deltaord = ord - ordprev[0];
if (deltaord > 1 || deltaord < -1)
{
printf ("-- order change %d\n", deltaord);
GSL_ERROR_NULL ("msbdf_apply too large order change", GSL_ESANITY);
}
/* Modify Nordsieck matrix if order or step length has been changed */
/* If order increased by 1, adjust Nordsieck matrix */
if (deltaord == 1)
{
if (ord > 2)
{
size_t i, j;
double hsum = h;
double coeff1 = -1.0;
double coeff2 = 1.0;
double hrelprev = 1.0;
double hrelprod = 1.0;
double hrel = 0.0;
/* Calculate coefficients used in adjustment to l */
DBL_ZERO_MEMSET (l, MSBDF_MAX_ORD + 1);
l[2] = 1.0;
for (i = 1; i < ord - 1; i++)
{
hsum += hprev[i];
hrel = hsum / h;
hrelprod *= hrel;
coeff1 -= 1.0 / (i + 1);
coeff2 += 1.0 / hrel;
for (j = i + 2; j > 1; j--)
{
l[j] *= hrelprev;
l[j] += l[j - 1];
}
hrelprev = hrel;
}
/* Scale Nordsieck matrix */
{
const double c = (-coeff1 - coeff2) / hrelprod;
for (i = 0; i < dim; i++)
{
z[ord * dim + i] = c * gsl_vector_get (abscor, i);
}
}
for (i = 2; i < ord; i++)
for (j = 0; j < dim; j++)
{
z[i * dim + j] += l[i] * z[ord * dim + j];
}
}
else
{
/* zero new vector for order incease from 1 to 2 */
DBL_ZERO_MEMSET (&z[ord * dim], dim);
}
#ifdef DEBUG
printf ("-- order increase detected, Nordsieck modified\n");
#endif
}
/* If order decreased by 1, adjust Nordsieck matrix */
if (deltaord == -1)
{
size_t i, j;
double hsum = 0.0;
/* Calculate coefficients used in adjustment to l */
DBL_ZERO_MEMSET (l, MSBDF_MAX_ORD + 1);
l[2] = 1.0;
for (i = 1; i < ord; i++)
{
hsum += hprev[i - 1];
for (j = i + 2; j > 1; j--)
{
l[j] *= hsum / h;
l[j] += l[j - 1];
}
}
/* Scale Nordsieck matrix */
for (i = 2; i < ord + 1; i++)
for (j = 0; j < dim; j++)
{
z[i * dim + j] += -l[i] * z[(ord + 1) * dim + j];
}
#ifdef DEBUG
printf ("-- order decrease detected, Nordsieck modified\n");
#endif
}
/* Scale Nordsieck vectors if step size has been changed */
if (state->ni > 0 && h != hprev[0])
{
size_t i, j;
const double hrel = h / hprev[0];
double coeff = hrel;
for (i = 1; i < ord + 1; i++)
{
for (j = 0; j < dim; j++)
{
z[i * dim + j] *= coeff;
}
coeff *= hrel;
}
#ifdef DEBUG
printf ("-- h != hprev, Nordsieck modified\n");
#endif
}
/* Calculate polynomial coefficients (l), error coefficient and
auxiliary coefficients
*/
msbdf_calccoeffs (ord, state->ordwait, h, hprev, l, &errcoeff,
º1coeff, &ordp1coeff, &ordp2coeff, &gamma);
/* Carry out the prediction step */
{
size_t i, j, k;
for (i = 1; i < ord + 1; i++)
for (j = ord; j > i - 1; j--)
for (k = 0; k < dim; k++)
{
z[(j - 1) * dim + k] += z[j * dim + k];
}
#ifdef DEBUG
{
size_t di;
printf ("-- predicted y: ");
for (di = 0; di < dim; di++)
{
printf ("%.5e ", z[di]);
}
printf ("\n");
}
#endif
}
/* Calculate correction step to abscor */
{
int s;
s = msbdf_corrector (vstate, sys, t, h, dim, z, errlev, l, errcoeff,
abscor, relcor, ytmp, ytmp2,
state->dfdy, state->dfdt, state->M,
state->p, state->rhs,
&(state->nJ), &(state->nM),
state->tprev, state->failt, gamma,
state->gammaprev, hprev[0]);
if (s != GSL_SUCCESS)
{
return s;
}
}
{
/* Add accepted final correction step to Nordsieck matrix */
size_t i, j;
for (i = 0; i < ord + 1; i++)
for (j = 0; j < dim; j++)
{
z[i * dim + j] += l[i] * gsl_vector_get (abscor, j);
}
#ifdef DEBUG
{
size_t di;
printf ("---- l: ");
for (di = 0; di < ord + 1; di++)
{
printf ("%.5e ", l[di]);
}
printf ("\n");
printf ("-- corrected y: ");
for (di = 0; di < dim; di++)
{
printf ("%.5e ", z[di]);
}
printf ("\n");
}
#endif
/* Derivatives at output */
if (dydt_out != NULL)
{
int s = GSL_ODEIV_FN_EVAL (sys, t + h, z, dydt_out);
if (s == GSL_EBADFUNC)
{
return s;
}
if (s != GSL_SUCCESS)
{
msbdf_failurehandler (vstate, dim, t);
#ifdef DEBUG
printf ("-- FAIL at user function evaluation\n");
#endif
return s;
}
}
/* Calculate error estimate */
for (i = 0; i < dim; i++)
{
yerr[i] = fabs (gsl_vector_get (abscor, i)) * errcoeff;
}
#ifdef DEBUG
{
size_t di;
printf ("-- yerr: ");
for (di = 0; di < dim; di++)
{
printf ("%.5e ", yerr[di]);
}
printf ("\n");
}
#endif
/* Save y values */
for (i = 0; i < dim; i++)
{
y[i] = z[0 * dim + i];
}
}
/* Scale abscor with errlev for later use in norm calculations
in order evaluation in msbdf_eval_order
*/
{
size_t i;
for (i = 0; i < dim; i++)
{
gsl_vector_set (abscorscaled, i,
gsl_vector_get (abscor, i) / errlev[i]);
}
}
/* Save items needed for evaluation of order increase on next
call, if needed
*/
if (state->ordwait == 1 && ord < MSBDF_MAX_ORD)
{
size_t i;
state->ordp1coeffprev = ordp1coeff;
for (i = 0; i < dim; i++)
{
gsl_vector_set (svec, i, gsl_vector_get (abscorscaled, i));
}
}
/* Consider and execute order change for next step if order is unchanged. */
if (state->ordwait == 0)
{
if (ord - ordprev[0] == 0)
{
msbdf_eval_order (abscorscaled, tempvec, svec, errcoeff, dim, errlev,
ordm1coeff, ordp1coeff,
state->ordp1coeffprev, ordp2coeff,
hprev, h, z, &(state->ord));
state->ordwait = state->ord + 2;
}
else
{
/* Postpone order evaluation if order has been modified elsewhere */
state->ordwait = 2;
}
}
/* Save information about current step in state and update counters */
{
size_t i;
for (i = MSBDF_MAX_ORD - 1; i > 0; i--)
{
hprev[i] = hprev[i - 1];
ordprev[i] = ordprev[i - 1];
}
}
hprev[0] = h;
ordprev[0] = ord;
#ifdef DEBUG
{
size_t di;
printf ("-- hprev: ");
for (di = 0; di < MSBDF_MAX_ORD; di++)
{
printf ("%.5e ", hprev[di]);
}
printf ("\n");
}
#endif
state->tprev = t;
state->ordwait--;
state->ni++;
state->gammaprev = gamma;
state->nJ++;
state->nM++;
#ifdef DEBUG
printf ("-- nJ=%d, nM=%d\n", (int) state->nJ, (int) state->nM);
#endif
}
return GSL_SUCCESS;
}
static int
msbdf_set_driver (void *vstate, const gsl_odeiv2_driver * d)
{
msbdf_state_t *state = (msbdf_state_t *) vstate;
state->driver = d;
return GSL_SUCCESS;
}
static int
msbdf_reset (void *vstate, size_t dim)
{
msbdf_state_t *state = (msbdf_state_t *) vstate;
size_t i;
state->ni = 0;
state->ord = 1;
state->ordwait = 2;
state->ordwaitbackup = 2;
state->failord = 0;
state->failt = GSL_NAN;
state->tprev = GSL_NAN;
state->gammaprev = 1.0;
state->gammaprevbackup = 1.0;
state->nJ = 0;
state->nM = 0;
state->failcount = 0;
state->ordp1coeffprev = 0.0;
DBL_ZERO_MEMSET (state->hprev, MSBDF_MAX_ORD);
DBL_ZERO_MEMSET (state->hprevbackup, MSBDF_MAX_ORD);
DBL_ZERO_MEMSET (state->z, (MSBDF_MAX_ORD + 1) * dim);
DBL_ZERO_MEMSET (state->zbackup, (MSBDF_MAX_ORD + 1) * dim);
for (i = 0; i < MSBDF_MAX_ORD; i++)
{
state->ordprev[i] = 1;
state->ordprevbackup[i] = 1;
}
#ifdef DEBUG
printf ("-- msbdf_reset called\n");
#endif
return GSL_SUCCESS;
}
static unsigned int
msbdf_order (void *vstate)
{
msbdf_state_t *state = (msbdf_state_t *) vstate;
return state->ord;
}
static void
msbdf_free (void *vstate)
{
msbdf_state_t *state = (msbdf_state_t *) vstate;
gsl_vector_free (state->rhs);
gsl_permutation_free (state->p);
gsl_matrix_free (state->M);
free (state->dfdt);
gsl_matrix_free (state->dfdy);
gsl_vector_free (state->tempvec);
gsl_vector_free (state->svec);
gsl_vector_free (state->relcor);
gsl_vector_free (state->abscor);
gsl_vector_free (state->abscorscaled);
free (state->errlev);
free (state->ordprevbackup);
free (state->ordprev);
free (state->hprevbackup);
free (state->hprev);
free (state->l);
free (state->ytmp2);
free (state->ytmp);
free (state->zbackup);
free (state->z);
free (state);
}
static const gsl_odeiv2_step_type msbdf_type = {
"msbdf", /* name */
1, /* can use dydt_in? */
1, /* gives exact dydt_out? */
&msbdf_alloc,
&msbdf_apply,
&msbdf_set_driver,
&msbdf_reset,
&msbdf_order,
&msbdf_free
};
const gsl_odeiv2_step_type *gsl_odeiv2_step_msbdf = &msbdf_type;