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

   1 /* multiroots/dnewton.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
  22 #include <stddef.h>
  23 #include <stdlib.h>
  24 #include <stdio.h>
  25 #include <math.h>
  26 #include <float.h>
  27
  28 #include <gsl/gsl_math.h>
  29 #include <gsl/gsl_errno.h>
  30 #include <gsl/gsl_multiroots.h>
  31 #include <gsl/gsl_linalg.h>
  32
  33 /* Newton method using a finite difference approximation to the jacobian.
  34    The derivatives are estimated using a step size of
  35
  36    h_i = sqrt(DBL_EPSILON) * x_i
  37
  38    */
  39
  40 typedef struct
  41   {
  42     gsl_matrix * J;
  43     gsl_matrix * lu;
  44     gsl_permutation * permutation;
  45   }
  46 dnewton_state_t;
  47
  48 static int dnewton_alloc (void * vstate, size_t n);
  49 static int dnewton_set (void * vstate, gsl_multiroot_function * function, gsl_vector * x, gsl_vector * f, gsl_vector * dx);
  50 static int dnewton_iterate (void * vstate, gsl_multiroot_function * function, gsl_vector * x, gsl_vector * f, gsl_vector * dx);
  51 static void dnewton_free (void * vstate);
  52
  53 static int
  54 dnewton_alloc (void * vstate, size_t n)
  55 {
  56   dnewton_state_t * state = (dnewton_state_t *) vstate;
  57   gsl_permutation * p;
  58   gsl_matrix * m, * J;
  59
  60   m = gsl_matrix_calloc (n,n);
  61
  62   if (m == 0)
  63     {
  64       GSL_ERROR ("failed to allocate space for lu", GSL_ENOMEM);
  65     }
  66
  67   state->lu = m ;
  68
  69   p = gsl_permutation_calloc (n);
  70
  71   if (p == 0)
  72     {
  73       gsl_matrix_free(m);
  74
  75       GSL_ERROR ("failed to allocate space for permutation", GSL_ENOMEM);
  76     }
  77
  78   state->permutation = p ;
  79
  80   J = gsl_matrix_calloc (n,n);
  81
  82   if (J == 0)
  83     {
  84       gsl_permutation_free(p);
  85       gsl_matrix_free(m);
  86
  87       GSL_ERROR ("failed to allocate space for d", GSL_ENOMEM);
  88     }
  89
  90   state->J = J;
  91
  92   return GSL_SUCCESS;
  93 }
  94
  95 static int
  96 dnewton_set (void * vstate, gsl_multiroot_function * function, gsl_vector * x, gsl_vector * f, gsl_vector * dx)
  97 {
  98   dnewton_state_t * state = (dnewton_state_t *) vstate;
  99   size_t i, n = function->n ;
 100   int status;
 101
 102   status = GSL_MULTIROOT_FN_EVAL (function, x, f);
 103   if (status)
 104     return status;
 105
 106   status = gsl_multiroot_fdjacobian (function, x, f, GSL_SQRT_DBL_EPSILON,
 107       state->J);
 108   if (status)
 109     return status;
 110
 111   for (i = 0; i < n; i++)
 112     {
 113       gsl_vector_set (dx, i, 0.0);
 114     }
 115
 116   return GSL_SUCCESS;
 117 }
 118
 119 static int
 120 dnewton_iterate (void * vstate, gsl_multiroot_function * function, gsl_vector * x, gsl_vector * f, gsl_vector * dx)
 121 {
 122   dnewton_state_t * state = (dnewton_state_t *) vstate;
 123
 124   int signum ;
 125
 126   size_t i;
 127
 128   size_t n = function->n ;
 129
 130   gsl_matrix_memcpy (state->lu, state->J);
 131
 132   {
 133     int status = gsl_linalg_LU_decomp (state->lu, state->permutation, &signum);
 134     if (status)
 135       return status;
 136   }
 137
 138   {
 139     int status = gsl_linalg_LU_solve (state->lu, state->permutation, f, dx);
 140     if (status)
 141       return status;
 142   }
 143
 144   for (i = 0; i < n; i++)
 145     {
 146       double e = gsl_vector_get (dx, i);
 147       double y = gsl_vector_get (x, i);
 148       gsl_vector_set (dx, i, -e);
 149       gsl_vector_set (x, i, y - e);
 150     }
 151
 152   {
 153     int status = GSL_MULTIROOT_FN_EVAL (function, x, f);
 154
 155     if (status != GSL_SUCCESS)
 156       {
 157         return GSL_EBADFUNC;
 158       }
 159   }
 160
 161   gsl_multiroot_fdjacobian (function, x, f, GSL_SQRT_DBL_EPSILON, state->J);
 162
 163   return GSL_SUCCESS;
 164 }
 165
 166
 167 static void
 168 dnewton_free (void * vstate)
 169 {
 170   dnewton_state_t * state = (dnewton_state_t *) vstate;
 171
 172   gsl_matrix_free(state->J);
 173   gsl_matrix_free(state->lu);
 174   gsl_permutation_free(state->permutation);
 175 }
 176
 177
 178 static const gsl_multiroot_fsolver_type dnewton_type =
 179 {"dnewton",                             /* name */
 180  sizeof (dnewton_state_t),
 181  &dnewton_alloc,
 182  &dnewton_set,
 183  &dnewton_iterate,
 184  &dnewton_free};
 185
 186 const gsl_multiroot_fsolver_type  * gsl_multiroot_fsolver_dnewton = &dnewton_type;