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

   1 /* ode-initval/cscal.c
   2  *
   3  * Copyright (C) 2002, 2007 Brian Gough
   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 #include <config.h>
  21 #include <stdlib.h>
  22 #include <string.h>
  23 #include <math.h>
  24 #include <gsl/gsl_errno.h>
  25 #include <gsl/gsl_math.h>
  26 #include <gsl/gsl_odeiv.h>
  27
  28 typedef struct
  29 {
  30   double eps_abs;
  31   double eps_rel;
  32   double a_y;
  33   double a_dydt;
  34   double * scale_abs;
  35 }
  36 sc_control_state_t;
  37
  38 static void *
  39 sc_control_alloc (void)
  40 {
  41   sc_control_state_t * s =
  42     (sc_control_state_t *) malloc (sizeof(sc_control_state_t));
  43
  44   if (s == 0)
  45     {
  46       GSL_ERROR_NULL ("failed to allocate space for sc_control_state",
  47                       GSL_ENOMEM);
  48     }
  49
  50   return s;
  51 }
  52
  53 static int
  54 sc_control_init (void * vstate,
  55                   double eps_abs, double eps_rel, double a_y, double a_dydt)
  56 {
  57   sc_control_state_t * s = (sc_control_state_t *) vstate;
  58
  59   if (eps_abs < 0)
  60     {
  61       GSL_ERROR ("eps_abs is negative", GSL_EINVAL);
  62     }
  63   else if (eps_rel < 0)
  64     {
  65       GSL_ERROR ("eps_rel is negative", GSL_EINVAL);
  66     }
  67   else if (a_y < 0)
  68     {
  69       GSL_ERROR ("a_y is negative", GSL_EINVAL);
  70     }
  71   else if (a_dydt < 0)
  72     {
  73       GSL_ERROR ("a_dydt is negative", GSL_EINVAL);
  74     }
  75
  76   s->eps_rel = eps_rel;
  77   s->eps_abs = eps_abs;
  78   s->a_y = a_y;
  79   s->a_dydt = a_dydt;
  80
  81   return GSL_SUCCESS;
  82 }
  83
  84 static int
  85 sc_control_hadjust(void * vstate, size_t dim, unsigned int ord, const double y[], const double yerr[], const double yp[], double * h)
  86 {
  87   sc_control_state_t *state = (sc_control_state_t *) vstate;
  88
  89   const double eps_abs = state->eps_abs;
  90   const double eps_rel = state->eps_rel;
  91   const double a_y     = state->a_y;
  92   const double a_dydt  = state->a_dydt;
  93   const double * scale_abs = state->scale_abs;
  94
  95   const double S = 0.9;
  96   const double h_old = *h;
  97
  98   double rmax = DBL_MIN;
  99   size_t i;
 100
 101   for(i=0; i<dim; i++) {
 102     const double D0 =
 103       eps_rel * (a_y * fabs(y[i]) + a_dydt * fabs(h_old * yp[i]))
 104       + eps_abs * scale_abs[i];
 105     const double r  = fabs(yerr[i]) / fabs(D0);
 106     rmax = GSL_MAX_DBL(r, rmax);
 107   }
 108
 109   if(rmax > 1.1) {
 110     /* decrease step, no more than factor of 5, but a fraction S more
 111        than scaling suggests (for better accuracy) */
 112     double r =  S / pow(rmax, 1.0/ord);
 113
 114     if (r < 0.2)
 115       r = 0.2;
 116
 117     *h = r * h_old;
 118
 119     return GSL_ODEIV_HADJ_DEC;
 120   }
 121   else if(rmax < 0.5) {
 122     /* increase step, no more than factor of 5 */
 123     double r = S / pow(rmax, 1.0/(ord+1.0));
 124
 125     if (r > 5.0)
 126       r = 5.0;
 127
 128     if (r < 1.0)  /* don't allow any decrease caused by S<1 */
 129       r = 1.0;
 130
 131     *h = r * h_old;
 132
 133     return GSL_ODEIV_HADJ_INC;
 134   }
 135   else {
 136     /* no change */
 137     return GSL_ODEIV_HADJ_NIL;
 138   }
 139 }
 140
 141
 142 static void
 143 sc_control_free (void * vstate)
 144 {
 145   sc_control_state_t *state = (sc_control_state_t *) vstate;
 146   free (state->scale_abs);
 147   free (state);
 148 }
 149
 150 static const gsl_odeiv_control_type sc_control_type =
 151 {"scaled",                      /* name */
 152  &sc_control_alloc,
 153  &sc_control_init,
 154  &sc_control_hadjust,
 155  &sc_control_free};
 156
 157 const gsl_odeiv_control_type *gsl_odeiv_control_scaled = &sc_control_type;
 158
 159
 160 gsl_odeiv_control *
 161 gsl_odeiv_control_scaled_new(double eps_abs, double eps_rel,
 162                              double a_y, double a_dydt,
 163                              const double scale_abs[],
 164                              size_t dim)
 165 {
 166   gsl_odeiv_control * c =
 167     gsl_odeiv_control_alloc (gsl_odeiv_control_scaled);
 168
 169   int status = gsl_odeiv_control_init (c, eps_abs, eps_rel, a_y, a_dydt);
 170
 171   if (status != GSL_SUCCESS)
 172     {
 173       gsl_odeiv_control_free (c);
 174       GSL_ERROR_NULL ("error trying to initialize control", status);
 175     }
 176
 177   {
 178     sc_control_state_t * s = (sc_control_state_t *) c->state;
 179
 180     s->scale_abs = (double *)malloc(dim * sizeof(double));
 181
 182     if (s->scale_abs == 0)
 183       {
 184         free (s);
 185         GSL_ERROR_NULL ("failed to allocate space for scale_abs",
 186                         GSL_ENOMEM);
 187       }
 188
 189     memcpy(s->scale_abs, scale_abs, dim * sizeof(double));
 190   }
 191
 192   return c;
 193 }