htswanalysis/MACS/lib/gsl/gsl-1.11/multiroots/fsolver.c

   1 /* multiroots/fsolver.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 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 <gsl/gsl_errno.h>
  24 #include <gsl/gsl_multiroots.h>
  25
  26 gsl_multiroot_fsolver *
  27 gsl_multiroot_fsolver_alloc (const gsl_multiroot_fsolver_type * T,
  28                                      size_t n)
  29 {
  30   int status;
  31
  32   gsl_multiroot_fsolver * s;
  33
  34   s = (gsl_multiroot_fsolver *) malloc (sizeof (gsl_multiroot_fsolver));
  35
  36   if (s == 0)
  37     {
  38       GSL_ERROR_VAL ("failed to allocate space for multiroot solver struct",
  39                         GSL_ENOMEM, 0);
  40     }
  41
  42   s->x = gsl_vector_calloc (n);
  43
  44   if (s->x == 0)
  45     {
  46       free (s);
  47       GSL_ERROR_VAL ("failed to allocate space for x", GSL_ENOMEM, 0);
  48     }
  49
  50   s->f = gsl_vector_calloc (n);
  51
  52   if (s->f == 0)
  53     {
  54       gsl_vector_free (s->x);
  55       free (s);
  56       GSL_ERROR_VAL ("failed to allocate space for f", GSL_ENOMEM, 0);
  57     }
  58
  59   s->dx = gsl_vector_calloc (n);
  60
  61   if (s->dx == 0)
  62     {
  63       gsl_vector_free (s->x);
  64       gsl_vector_free (s->f);
  65       free (s);
  66       GSL_ERROR_VAL ("failed to allocate space for dx", GSL_ENOMEM, 0);
  67     }
  68
  69   s->state = malloc (T->size);
  70
  71   if (s->state == 0)
  72     {
  73       gsl_vector_free (s->dx);
  74       gsl_vector_free (s->x);
  75       gsl_vector_free (s->f);
  76       free (s);         /* exception in constructor, avoid memory leak */
  77
  78       GSL_ERROR_VAL ("failed to allocate space for multiroot solver state",
  79                         GSL_ENOMEM, 0);
  80     }
  81
  82   s->type = T ;
  83
  84   status = (s->type->alloc)(s->state, n);
  85
  86   if (status != GSL_SUCCESS)
  87     {
  88       (s->type->free)(s->state);
  89       free (s->state);
  90       gsl_vector_free (s->dx);
  91       gsl_vector_free (s->x);
  92       gsl_vector_free (s->f);
  93       free (s);         /* exception in constructor, avoid memory leak */
  94
  95       GSL_ERROR_VAL ("failed to set solver", status, 0);
  96     }
  97
  98   s->function = NULL;
  99
 100   return s;
 101 }
 102
 103 int
 104 gsl_multiroot_fsolver_set (gsl_multiroot_fsolver * s,
 105                            gsl_multiroot_function * f,
 106                            const gsl_vector * x)
 107 {
 108   if (s->x->size != f->n)
 109     {
 110       GSL_ERROR ("function incompatible with solver size", GSL_EBADLEN);
 111     }
 112
 113   if (x->size != f->n)
 114     {
 115       GSL_ERROR ("vector length not compatible with function", GSL_EBADLEN);
 116     }
 117
 118   s->function = f;
 119   gsl_vector_memcpy(s->x,x);
 120
 121   return (s->type->set) (s->state, s->function, s->x, s->f, s->dx);
 122 }
 123
 124 int
 125 gsl_multiroot_fsolver_iterate (gsl_multiroot_fsolver * s)
 126 {
 127   return (s->type->iterate) (s->state, s->function, s->x, s->f, s->dx);
 128 }
 129
 130 void
 131 gsl_multiroot_fsolver_free (gsl_multiroot_fsolver * s)
 132 {
 133   (s->type->free) (s->state);
 134   free (s->state);
 135   gsl_vector_free (s->dx);
 136   gsl_vector_free (s->x);
 137   gsl_vector_free (s->f);
 138   free (s);
 139 }
 140
 141 const char *
 142 gsl_multiroot_fsolver_name (const gsl_multiroot_fsolver * s)
 143 {
 144   return s->type->name;
 145 }
 146
 147 gsl_vector *
 148 gsl_multiroot_fsolver_root (const gsl_multiroot_fsolver * s)
 149 {
 150   return s->x;
 151 }
 152
 153 gsl_vector *
 154 gsl_multiroot_fsolver_dx (const gsl_multiroot_fsolver * s)
 155 {
 156   return s->dx;
 157 }
 158
 159 gsl_vector *
 160 gsl_multiroot_fsolver_f (const gsl_multiroot_fsolver * s)
 161 {
 162   return s->f;
 163 }