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

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