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

   1 /* multiroots/gnewton.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 #include "enorm.c"
  34
  35 /* Simple globally convergent Newton method (rejects uphill steps) */
  36
  37 typedef struct
  38   {
  39     double phi;
  40     gsl_vector * x_trial;
  41     gsl_vector * d;
  42     gsl_matrix * lu;
  43     gsl_permutation * permutation;
  44   }
  45 gnewton_state_t;
  46
  47 static int gnewton_alloc (void * vstate, size_t n);
  48 static int gnewton_set (void * vstate, gsl_multiroot_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx);
  49 static int gnewton_iterate (void * vstate, gsl_multiroot_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx);
  50 static void gnewton_free (void * vstate);
  51
  52 static int
  53 gnewton_alloc (void * vstate, size_t n)
  54 {
  55   gnewton_state_t * state = (gnewton_state_t *) vstate;
  56   gsl_vector * d, * x_trial ;
  57   gsl_permutation * p;
  58   gsl_matrix * m;
  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   d = gsl_vector_calloc (n);
  81
  82   if (d == 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->d = d;
  91
  92   x_trial = gsl_vector_calloc (n);
  93
  94   if (x_trial == 0)
  95     {
  96       gsl_vector_free(d);
  97       gsl_permutation_free(p);
  98       gsl_matrix_free(m);
  99
 100       GSL_ERROR ("failed to allocate space for x_trial", GSL_ENOMEM);
 101     }
 102
 103   state->x_trial = x_trial;
 104
 105   return GSL_SUCCESS;
 106 }
 107
 108
 109 static int
 110 gnewton_set (void * vstate, gsl_multiroot_function_fdf * FDF, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx)
 111 {
 112   gnewton_state_t * state = (gnewton_state_t *) vstate;
 113   size_t i, n = FDF->n ;
 114
 115   GSL_MULTIROOT_FN_EVAL_F_DF (FDF, x, f, J);
 116
 117   for (i = 0; i < n; i++)
 118     {
 119       gsl_vector_set (dx, i, 0.0);
 120     }
 121
 122   state->phi = enorm(f);
 123
 124   return GSL_SUCCESS;
 125 }
 126
 127 static int
 128 gnewton_iterate (void * vstate, gsl_multiroot_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx)
 129 {
 130   gnewton_state_t * state = (gnewton_state_t *) vstate;
 131
 132   int signum ;
 133   double t, phi0, phi1;
 134
 135   size_t i;
 136
 137   size_t n = fdf->n ;
 138
 139   gsl_matrix_memcpy (state->lu, J);
 140
 141   gsl_linalg_LU_decomp (state->lu, state->permutation, &signum);
 142
 143   gsl_linalg_LU_solve (state->lu, state->permutation, f, state->d);
 144
 145   t = 1;
 146
 147   phi0 = state->phi;
 148
 149 new_step:
 150
 151   for (i = 0; i < n; i++)
 152     {
 153       double di = gsl_vector_get (state->d, i);
 154       double xi = gsl_vector_get (x, i);
 155       gsl_vector_set (state->x_trial, i, xi - t*di);
 156     }
 157
 158   {
 159     int status = GSL_MULTIROOT_FN_EVAL_F (fdf, state->x_trial, f);
 160
 161     if (status != GSL_SUCCESS)
 162       {
 163         return GSL_EBADFUNC;
 164       }
 165   }
 166
 167   phi1 = enorm (f);
 168
 169   if (phi1 > phi0 && t > GSL_DBL_EPSILON)
 170     {
 171       /* full step goes uphill, take a reduced step instead */
 172
 173       double theta = phi1 / phi0;
 174       double u = (sqrt(1.0 + 6.0 * theta) - 1.0) / (3.0 * theta);
 175
 176       t *= u ;
 177
 178       goto new_step;
 179     }
 180
 181   /* copy x_trial into x */
 182
 183   gsl_vector_memcpy (x, state->x_trial);
 184
 185   for (i = 0; i < n; i++)
 186     {
 187       double di = gsl_vector_get (state->d, i);
 188       gsl_vector_set (dx, i, -t*di);
 189     }
 190
 191   {
 192     int status = GSL_MULTIROOT_FN_EVAL_DF (fdf, x, J);
 193
 194     if (status != GSL_SUCCESS)
 195       {
 196         return GSL_EBADFUNC;
 197       }
 198   }
 199
 200   state->phi = phi1;
 201
 202   return GSL_SUCCESS;
 203 }
 204
 205
 206 static void
 207 gnewton_free (void * vstate)
 208 {
 209   gnewton_state_t * state = (gnewton_state_t *) vstate;
 210
 211   gsl_vector_free(state->d);
 212   gsl_vector_free(state->x_trial);
 213   gsl_matrix_free(state->lu);
 214   gsl_permutation_free(state->permutation);
 215 }
 216
 217
 218 static const gsl_multiroot_fdfsolver_type gnewton_type =
 219 {"gnewton",                             /* name */
 220  sizeof (gnewton_state_t),
 221  &gnewton_alloc,
 222  &gnewton_set,
 223  &gnewton_iterate,
 224  &gnewton_free};
 225
 226 const gsl_multiroot_fdfsolver_type  * gsl_multiroot_fdfsolver_gnewton = &gnewton_type;