htswanalysis/MACS/lib/gsl/gsl-1.11/ode-initval/evolve.c

   1 /* ode-initval/evolve.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 /* Author:  G. Jungman
  21  */
  22 #include <config.h>
  23 #include <string.h>
  24 #include <stdlib.h>
  25 #include <gsl/gsl_math.h>
  26 #include <gsl/gsl_errno.h>
  27 #include <gsl/gsl_odeiv.h>
  28
  29 #include "odeiv_util.h"
  30
  31 gsl_odeiv_evolve *
  32 gsl_odeiv_evolve_alloc (size_t dim)
  33 {
  34   gsl_odeiv_evolve *e =
  35     (gsl_odeiv_evolve *) malloc (sizeof (gsl_odeiv_evolve));
  36
  37   if (e == 0)
  38     {
  39       GSL_ERROR_NULL ("failed to allocate space for evolve struct",
  40                       GSL_ENOMEM);
  41     }
  42
  43   e->y0 = (double *) malloc (dim * sizeof (double));
  44
  45   if (e->y0 == 0)
  46     {
  47       free (e);
  48       GSL_ERROR_NULL ("failed to allocate space for y0", GSL_ENOMEM);
  49     }
  50
  51   e->yerr = (double *) malloc (dim * sizeof (double));
  52
  53   if (e->yerr == 0)
  54     {
  55       free (e->y0);
  56       free (e);
  57       GSL_ERROR_NULL ("failed to allocate space for yerr", GSL_ENOMEM);
  58     }
  59
  60   e->dydt_in = (double *) malloc (dim * sizeof (double));
  61
  62   if (e->dydt_in == 0)
  63     {
  64       free (e->yerr);
  65       free (e->y0);
  66       free (e);
  67       GSL_ERROR_NULL ("failed to allocate space for dydt_in", GSL_ENOMEM);
  68     }
  69
  70   e->dydt_out = (double *) malloc (dim * sizeof (double));
  71
  72   if (e->dydt_out == 0)
  73     {
  74       free (e->dydt_in);
  75       free (e->yerr);
  76       free (e->y0);
  77       free (e);
  78       GSL_ERROR_NULL ("failed to allocate space for dydt_out", GSL_ENOMEM);
  79     }
  80
  81   e->dimension = dim;
  82   e->count = 0;
  83   e->failed_steps = 0;
  84   e->last_step = 0.0;
  85
  86   return e;
  87 }
  88
  89 int
  90 gsl_odeiv_evolve_reset (gsl_odeiv_evolve * e)
  91 {
  92   e->count = 0;
  93   e->failed_steps = 0;
  94   e->last_step = 0.0;
  95   return GSL_SUCCESS;
  96 }
  97
  98 void
  99 gsl_odeiv_evolve_free (gsl_odeiv_evolve * e)
 100 {
 101   free (e->dydt_out);
 102   free (e->dydt_in);
 103   free (e->yerr);
 104   free (e->y0);
 105   free (e);
 106 }
 107
 108 /* Evolution framework method.
 109  *
 110  * Uses an adaptive step control object
 111  */
 112 int
 113 gsl_odeiv_evolve_apply (gsl_odeiv_evolve * e,
 114                         gsl_odeiv_control * con,
 115                         gsl_odeiv_step * step,
 116                         const gsl_odeiv_system * dydt,
 117                         double *t, double t1, double *h, double y[])
 118 {
 119   const double t0 = *t;
 120   double h0 = *h;
 121   int step_status;
 122   int final_step = 0;
 123   double dt = t1 - t0;  /* remaining time, possibly less than h */
 124
 125   if (e->dimension != step->dimension)
 126     {
 127       GSL_ERROR ("step dimension must match evolution size", GSL_EINVAL);
 128     }
 129
 130   if ((dt < 0.0 && h0 > 0.0) || (dt > 0.0 && h0 < 0.0))
 131     {
 132       GSL_ERROR ("step direction must match interval direction", GSL_EINVAL);
 133     }
 134
 135   /* No need to copy if we cannot control the step size. */
 136
 137   if (con != NULL)
 138     {
 139       DBL_MEMCPY (e->y0, y, e->dimension);
 140     }
 141
 142   /* Calculate initial dydt once if the method can benefit. */
 143
 144   if (step->type->can_use_dydt_in)
 145     {
 146       int status = GSL_ODEIV_FN_EVAL (dydt, t0, y, e->dydt_in);
 147
 148       if (status)
 149         {
 150           return status;
 151         }
 152     }
 153
 154 try_step:
 155
 156   if ((dt >= 0.0 && h0 > dt) || (dt < 0.0 && h0 < dt))
 157     {
 158       h0 = dt;
 159       final_step = 1;
 160     }
 161   else
 162     {
 163       final_step = 0;
 164     }
 165
 166   if (step->type->can_use_dydt_in)
 167     {
 168       step_status =
 169         gsl_odeiv_step_apply (step, t0, h0, y, e->yerr, e->dydt_in,
 170                               e->dydt_out, dydt);
 171     }
 172   else
 173     {
 174       step_status =
 175         gsl_odeiv_step_apply (step, t0, h0, y, e->yerr, NULL, e->dydt_out,
 176                               dydt);
 177     }
 178
 179   /* Check for stepper internal failure */
 180
 181   if (step_status != GSL_SUCCESS)
 182     {
 183       *h = h0;  /* notify user of step-size which caused the failure */
 184       return step_status;
 185     }
 186
 187   e->count++;
 188   e->last_step = h0;
 189
 190   if (final_step)
 191     {
 192       *t = t1;
 193     }
 194   else
 195     {
 196       *t = t0 + h0;
 197     }
 198
 199   if (con != NULL)
 200     {
 201       /* Check error and attempt to adjust the step. */
 202
 203       double h_old = h0;
 204
 205       const int hadjust_status
 206         = gsl_odeiv_control_hadjust (con, step, y, e->yerr, e->dydt_out, &h0);
 207
 208       if (hadjust_status == GSL_ODEIV_HADJ_DEC)
 209         {
 210           /* Check that the reported status is correct (i.e. an actual
 211              decrease in h0 occured) and the suggested h0 will change
 212              the time by at least 1 ulp */
 213
 214           double t_curr = GSL_COERCE_DBL(*t);
 215           double t_next = GSL_COERCE_DBL((*t) + h0);
 216
 217           if (fabs(h0) < fabs(h_old) && t_next != t_curr)
 218             {
 219               /* Step was decreased. Undo step, and try again with new h0. */
 220               DBL_MEMCPY (y, e->y0, dydt->dimension);
 221               e->failed_steps++;
 222               goto try_step;
 223             }
 224           else
 225             {
 226               h0 = h_old; /* keep current step size */
 227             }
 228         }
 229     }
 230
 231   *h = h0;  /* suggest step size for next time-step */
 232
 233   return step_status;
 234 }