text: /* ode-initval/evolve.c
   * 
   * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
   * 
   * 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 2 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.
   */
  
  /* Author:  G. Jungman
   */
  #include <config.h>
  #include <string.h>
  #include <stdlib.h>
  #include <gsl/gsl_math.h>
  #include <gsl/gsl_errno.h>
  #include <gsl/gsl_odeiv.h>
  
  #include "odeiv_util.h"
  
  gsl_odeiv_evolve *
  gsl_odeiv_evolve_alloc (size_t dim)
  {
    gsl_odeiv_evolve *e =
      (gsl_odeiv_evolve *) malloc (sizeof (gsl_odeiv_evolve));
  
    if (e == 0)
      {
        GSL_ERROR_NULL ("failed to allocate space for evolve struct",
                        GSL_ENOMEM);
      }
  
    e->y0 = (double *) malloc (dim * sizeof (double));
  
    if (e->y0 == 0)
      {
        free (e);
        GSL_ERROR_NULL ("failed to allocate space for y0", GSL_ENOMEM);
      }
  
    e->yerr = (double *) malloc (dim * sizeof (double));
  
    if (e->yerr == 0)
      {
        free (e->y0);
        free (e);
        GSL_ERROR_NULL ("failed to allocate space for yerr", GSL_ENOMEM);
      }
  
    e->dydt_in = (double *) malloc (dim * sizeof (double));
  
    if (e->dydt_in == 0)
      {
        free (e->yerr);
        free (e->y0);
        free (e);
        GSL_ERROR_NULL ("failed to allocate space for dydt_in", GSL_ENOMEM);
      }
  
    e->dydt_out = (double *) malloc (dim * sizeof (double));
  
    if (e->dydt_out == 0)
      {
        free (e->dydt_in);
        free (e->yerr);
        free (e->y0);
        free (e);
        GSL_ERROR_NULL ("failed to allocate space for dydt_out", GSL_ENOMEM);
      }
  
    e->dimension = dim;
    e->count = 0;
    e->failed_steps = 0;
    e->last_step = 0.0;
  
    return e;
  }
  
  int
  gsl_odeiv_evolve_reset (gsl_odeiv_evolve * e)
  {
    e->count = 0;
    e->failed_steps = 0;
    e->last_step = 0.0;
    return GSL_SUCCESS;
  }
  
  void
  gsl_odeiv_evolve_free (gsl_odeiv_evolve * e)
  {
    free (e->dydt_out);
    free (e->dydt_in);
    free (e->yerr);
    free (e->y0);
    free (e);
  }
  
  /* Evolution framework method.
   *
   * Uses an adaptive step control object
   */
  int
  gsl_odeiv_evolve_apply (gsl_odeiv_evolve * e,
                          gsl_odeiv_control * con,
                          gsl_odeiv_step * step,
                          const gsl_odeiv_system * dydt,
                          double *t, double t1, double *h, double y[])
  {
    const double t0 = *t;
    double h0 = *h;
    int step_status;
    int final_step = 0;
    double dt = t1 - t0;  /* remaining time, possibly less than h */
  
    if (e->dimension != step->dimension)
      {
        GSL_ERROR ("step dimension must match evolution size", GSL_EINVAL);
      }
  
    if ((dt < 0.0 && h0 > 0.0) || (dt > 0.0 && h0 < 0.0))
      {
        GSL_ERROR ("step direction must match interval direction", GSL_EINVAL);
      }
  
    /* No need to copy if we cannot control the step size. */
  
    if (con != NULL)
      {
        DBL_MEMCPY (e->y0, y, e->dimension);
      }
  
    /* Calculate initial dydt once 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 (status) 
          {
            return status;
          }
      }
  
  try_step:
      
    if ((dt >= 0.0 && h0 > dt) || (dt < 0.0 && h0 < dt))
      {
        h0 = dt;
        final_step = 1;
      }
    else
      {
        final_step = 0;
      }
  
    if (step->type->can_use_dydt_in)
      {
        step_status =
          gsl_odeiv_step_apply (step, t0, h0, y, e->yerr, e->dydt_in,
                                e->dydt_out, dydt);
      }
    else
      {
        step_status =
          gsl_odeiv_step_apply (step, t0, h0, y, e->yerr, NULL, e->dydt_out,
                                dydt);
      }
  
    /* Check for stepper internal failure */
  
    if (step_status != GSL_SUCCESS) 
      {
        *h = h0;  /* notify user of step-size which caused the failure */
        return step_status;
      }
  
    e->count++;
    e->last_step = h0;
  
    if (final_step)
      {
        *t = t1;
      }
    else
      {
        *t = t0 + h0;
      }
  
    if (con != NULL)
      {
        /* Check error and attempt to adjust the step. */
        const int hadjust_status 
          = gsl_odeiv_control_hadjust (con, step, y, e->yerr, e->dydt_out, &h0);
  
        if (hadjust_status == GSL_ODEIV_HADJ_DEC)
          {
            /* Step was decreased. Undo and go back to try again. */
            DBL_MEMCPY (y, e->y0, dydt->dimension);
            e->failed_steps++;
            goto try_step;
          }
      }
  
    *h = h0;  /* suggest step size for next time-step */
  
    return step_status;
  }
  