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

   1 /* multiroots/fdfsolver.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_fdfsolver *
  27 gsl_multiroot_fdfsolver_alloc (const gsl_multiroot_fdfsolver_type * T,
  28                                        size_t n)
  29 {
  30   int status;
  31
  32   gsl_multiroot_fdfsolver * s;
  33
  34   s = (gsl_multiroot_fdfsolver *) malloc (sizeof (gsl_multiroot_fdfsolver));
  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->J = gsl_matrix_calloc (n,n);
  60
  61   if (s->J == 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 g", GSL_ENOMEM, 0);
  67     }
  68
  69   s->dx = gsl_vector_calloc (n);
  70
  71   if (s->dx == 0)
  72     {
  73       gsl_matrix_free (s->J);
  74       gsl_vector_free (s->x);
  75       gsl_vector_free (s->f);
  76       free (s);
  77       GSL_ERROR_VAL ("failed to allocate space for dx", GSL_ENOMEM, 0);
  78     }
  79
  80   s->state = malloc (T->size);
  81
  82   if (s->state == 0)
  83     {
  84       gsl_vector_free (s->dx);
  85       gsl_vector_free (s->x);
  86       gsl_vector_free (s->f);
  87       gsl_matrix_free (s->J);
  88       free (s);         /* exception in constructor, avoid memory leak */
  89
  90       GSL_ERROR_VAL ("failed to allocate space for multiroot solver state",
  91                         GSL_ENOMEM, 0);
  92     }
  93
  94   s->type = T ;
  95
  96   status = (s->type->alloc)(s->state, n);
  97
  98   if (status != GSL_SUCCESS)
  99     {
 100       free (s->state);
 101       gsl_vector_free (s->dx);
 102       gsl_vector_free (s->x);
 103       gsl_vector_free (s->f);
 104       gsl_matrix_free (s->J);
 105       free (s);         /* exception in constructor, avoid memory leak */
 106
 107       GSL_ERROR_VAL ("failed to set solver", status, 0);
 108     }
 109
 110   s->fdf = NULL;
 111
 112   return s;
 113 }
 114
 115 int
 116 gsl_multiroot_fdfsolver_set (gsl_multiroot_fdfsolver * s,
 117                              gsl_multiroot_function_fdf * f,
 118                              const gsl_vector * x)
 119 {
 120   if (s->x->size != f->n)
 121     {
 122       GSL_ERROR ("function incompatible with solver size", GSL_EBADLEN);
 123     }
 124
 125   if (x->size != f->n)
 126     {
 127       GSL_ERROR ("vector length not compatible with function", GSL_EBADLEN);
 128     }
 129
 130   s->fdf = f;
 131   gsl_vector_memcpy(s->x,x);
 132
 133   return (s->type->set) (s->state, s->fdf, s->x, s->f, s->J, s->dx);
 134 }
 135
 136 int
 137 gsl_multiroot_fdfsolver_iterate (gsl_multiroot_fdfsolver * s)
 138 {
 139   return (s->type->iterate) (s->state, s->fdf, s->x, s->f, s->J, s->dx);
 140 }
 141
 142 void
 143 gsl_multiroot_fdfsolver_free (gsl_multiroot_fdfsolver * s)
 144 {
 145   (s->type->free) (s->state);
 146   free (s->state);
 147   gsl_vector_free (s->dx);
 148   gsl_vector_free (s->x);
 149   gsl_vector_free (s->f);
 150   gsl_matrix_free (s->J);
 151   free (s);
 152 }
 153
 154 const char *
 155 gsl_multiroot_fdfsolver_name (const gsl_multiroot_fdfsolver * s)
 156 {
 157   return s->type->name;
 158 }
 159
 160 gsl_vector *
 161 gsl_multiroot_fdfsolver_root (const gsl_multiroot_fdfsolver * s)
 162 {
 163   return s->x;
 164 }
 165
 166 gsl_vector *
 167 gsl_multiroot_fdfsolver_dx (const gsl_multiroot_fdfsolver * s)
 168 {
 169   return s->dx;
 170 }
 171
 172 gsl_vector *
 173 gsl_multiroot_fdfsolver_f (const gsl_multiroot_fdfsolver * s)
 174 {
 175   return s->f;
 176 }