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

   1 /* multiroots/newton.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 typedef struct
  34   {
  35     gsl_matrix * lu;
  36     gsl_permutation * permutation;
  37   }
  38 newton_state_t;
  39
  40 static int newton_alloc (void * vstate, size_t n);
  41 static int newton_set (void * vstate, gsl_multiroot_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx);
  42 static int newton_iterate (void * vstate, gsl_multiroot_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx);
  43 static void newton_free (void * vstate);
  44
  45 static int
  46 newton_alloc (void * vstate, size_t n)
  47 {
  48   newton_state_t * state = (newton_state_t *) vstate;
  49   gsl_permutation * p;
  50   gsl_matrix * m;
  51
  52   m = gsl_matrix_calloc (n,n);
  53
  54   if (m == 0)
  55     {
  56       GSL_ERROR ("failed to allocate space for lu", GSL_ENOMEM);
  57     }
  58
  59   state->lu = m ;
  60
  61   p = gsl_permutation_calloc (n);
  62
  63   if (p == 0)
  64     {
  65       gsl_matrix_free(m);
  66
  67       GSL_ERROR ("failed to allocate space for permutation", GSL_ENOMEM);
  68     }
  69
  70   state->permutation = p ;
  71
  72   return GSL_SUCCESS;
  73 }
  74
  75 static int
  76 newton_set (void * vstate, gsl_multiroot_function_fdf * FDF, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx)
  77 {
  78   newton_state_t * state = (newton_state_t *) vstate;
  79
  80   size_t i, n = FDF->n ;
  81
  82   state = 0 ; /* avoid warnings about unused parameters */
  83
  84   GSL_MULTIROOT_FN_EVAL_F_DF (FDF, x, f, J);
  85
  86   for (i = 0; i < n; i++)
  87     {
  88       gsl_vector_set (dx, i, 0.0);
  89     }
  90
  91   return GSL_SUCCESS;
  92 }
  93
  94 static int
  95 newton_iterate (void * vstate, gsl_multiroot_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx)
  96 {
  97   newton_state_t * state = (newton_state_t *) vstate;
  98
  99   int signum;
 100
 101   size_t i;
 102
 103   size_t n = fdf->n ;
 104
 105   gsl_matrix_memcpy (state->lu, J);
 106
 107   gsl_linalg_LU_decomp (state->lu, state->permutation, &signum);
 108
 109   gsl_linalg_LU_solve (state->lu, state->permutation, f, dx);
 110
 111   for (i = 0; i < n; i++)
 112     {
 113       double e = gsl_vector_get (dx, i);
 114       double y = gsl_vector_get (x, i);
 115       gsl_vector_set (dx, i, -e);
 116       gsl_vector_set (x, i, y - e);
 117     }
 118
 119   {
 120     int status = GSL_MULTIROOT_FN_EVAL_F_DF (fdf, x, f, J);
 121
 122     if (status != GSL_SUCCESS)
 123       {
 124         return GSL_EBADFUNC;
 125       }
 126   }
 127
 128   return GSL_SUCCESS;
 129 }
 130
 131
 132 static void
 133 newton_free (void * vstate)
 134 {
 135   newton_state_t * state = (newton_state_t *) vstate;
 136
 137   gsl_matrix_free(state->lu);
 138
 139   gsl_permutation_free(state->permutation);
 140 }
 141
 142
 143 static const gsl_multiroot_fdfsolver_type newton_type =
 144 {"newton",                              /* name */
 145  sizeof (newton_state_t),
 146  &newton_alloc,
 147  &newton_set,
 148  &newton_iterate,
 149  &newton_free};
 150
 151 const gsl_multiroot_fdfsolver_type  * gsl_multiroot_fdfsolver_newton = &newton_type;