htswanalysis/MACS/lib/gsl/gsl-1.11/multimin/conjugate_fr.c

   1 /* multimin/conjugate_fr.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Fabrice Rossi
   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 /* conjugate_fr.c -- Conjugate gradient Fletcher-Reeve algorithm */
  21
  22 /* Modified by Brian Gough to use single iteration structure */
  23
  24 #include <config.h>
  25 #include <gsl/gsl_multimin.h>
  26 #include <gsl/gsl_blas.h>
  27
  28 #include "directional_minimize.c"
  29
  30 typedef struct
  31 {
  32   int iter;
  33   double step;
  34   double max_step;
  35   double tol;
  36   gsl_vector *x1;
  37   gsl_vector *dx1;
  38   gsl_vector *x2;
  39   double pnorm;
  40   gsl_vector *p;
  41   double g0norm;
  42   gsl_vector *g0;
  43 }
  44 conjugate_fr_state_t;
  45
  46 static int
  47 conjugate_fr_alloc (void *vstate, size_t n)
  48 {
  49   conjugate_fr_state_t *state = (conjugate_fr_state_t *) vstate;
  50
  51   state->x1 = gsl_vector_calloc (n);
  52
  53   if (state->x1 == 0)
  54     {
  55       GSL_ERROR ("failed to allocate space for x1", GSL_ENOMEM);
  56     }
  57
  58   state->dx1 = gsl_vector_calloc (n);
  59
  60   if (state->dx1 == 0)
  61     {
  62       gsl_vector_free (state->x1);
  63       GSL_ERROR ("failed to allocate space for dx1", GSL_ENOMEM);
  64     }
  65
  66   state->x2 = gsl_vector_calloc (n);
  67
  68   if (state->x2 == 0)
  69     {
  70       gsl_vector_free (state->dx1);
  71       gsl_vector_free (state->x1);
  72       GSL_ERROR ("failed to allocate space for x2", GSL_ENOMEM);
  73     }
  74
  75   state->p = gsl_vector_calloc (n);
  76
  77   if (state->p == 0)
  78     {
  79       gsl_vector_free (state->x2);
  80       gsl_vector_free (state->dx1);
  81       gsl_vector_free (state->x1);
  82       GSL_ERROR ("failed to allocate space for p", GSL_ENOMEM);
  83     }
  84
  85   state->g0 = gsl_vector_calloc (n);
  86
  87   if (state->g0 == 0)
  88     {
  89       gsl_vector_free (state->p);
  90       gsl_vector_free (state->x2);
  91       gsl_vector_free (state->dx1);
  92       gsl_vector_free (state->x1);
  93       GSL_ERROR ("failed to allocate space for g0", GSL_ENOMEM);
  94     }
  95
  96   return GSL_SUCCESS;
  97 }
  98
  99 static int
 100 conjugate_fr_set (void *vstate, gsl_multimin_function_fdf * fdf,
 101                   const gsl_vector * x, double *f, gsl_vector * gradient,
 102                   double step_size, double tol)
 103 {
 104   conjugate_fr_state_t *state = (conjugate_fr_state_t *) vstate;
 105
 106   state->iter = 0;
 107   state->step = step_size;
 108   state->max_step = step_size;
 109   state->tol  = tol;
 110
 111   GSL_MULTIMIN_FN_EVAL_F_DF (fdf, x, f, gradient);
 112
 113   /* Use the gradient as the initial direction */
 114
 115   gsl_vector_memcpy (state->p, gradient);
 116   gsl_vector_memcpy (state->g0, gradient);
 117
 118   {
 119     double gnorm = gsl_blas_dnrm2 (gradient);
 120
 121     state->pnorm = gnorm;
 122     state->g0norm = gnorm;
 123   }
 124
 125   return GSL_SUCCESS;
 126 }
 127
 128 static void
 129 conjugate_fr_free (void *vstate)
 130 {
 131   conjugate_fr_state_t *state = (conjugate_fr_state_t *) vstate;
 132
 133   gsl_vector_free (state->g0);
 134   gsl_vector_free (state->p);
 135   gsl_vector_free (state->x2);
 136   gsl_vector_free (state->dx1);
 137   gsl_vector_free (state->x1);
 138 }
 139
 140 static int
 141 conjugate_fr_restart (void *vstate)
 142 {
 143   conjugate_fr_state_t *state = (conjugate_fr_state_t *) vstate;
 144
 145   state->iter = 0;
 146   return GSL_SUCCESS;
 147 }
 148
 149 static int
 150 conjugate_fr_iterate (void *vstate, gsl_multimin_function_fdf * fdf,
 151                       gsl_vector * x, double *f,
 152                       gsl_vector * gradient, gsl_vector * dx)
 153 {
 154   conjugate_fr_state_t *state = (conjugate_fr_state_t *) vstate;
 155
 156   gsl_vector *x1 = state->x1;
 157   gsl_vector *dx1 = state->dx1;
 158   gsl_vector *x2 = state->x2;
 159   gsl_vector *p = state->p;
 160   gsl_vector *g0 = state->g0;
 161
 162   double pnorm = state->pnorm;
 163   double g0norm = state->g0norm;
 164
 165   double fa = *f, fb, fc;
 166   double dir;
 167   double stepa = 0.0, stepb, stepc = state->step, tol = state->tol;
 168
 169   double g1norm;
 170   double pg;
 171
 172   if (pnorm == 0.0 || g0norm == 0.0)
 173     {
 174       gsl_vector_set_zero (dx);
 175       return GSL_ENOPROG;
 176     }
 177
 178   /* Determine which direction is downhill, +p or -p */
 179
 180   gsl_blas_ddot (p, gradient, &pg);
 181
 182   dir = (pg >= 0.0) ? +1.0 : -1.0;
 183
 184   /* Compute new trial point at x_c= x - step * p, where p is the
 185      current direction */
 186
 187   take_step (x, p, stepc, dir / pnorm, x1, dx);
 188
 189   /* Evaluate function and gradient at new point xc */
 190
 191   fc = GSL_MULTIMIN_FN_EVAL_F (fdf, x1);
 192
 193   if (fc < fa)
 194     {
 195       /* Success, reduced the function value */
 196       state->step = stepc * 2.0;
 197       *f = fc;
 198       gsl_vector_memcpy (x, x1);
 199       GSL_MULTIMIN_FN_EVAL_DF (fdf, x1, gradient);
 200       return GSL_SUCCESS;
 201     }
 202
 203 #ifdef DEBUG
 204   printf ("got stepc = %g fc = %g\n", stepc, fc);
 205 #endif
 206
 207   /* Do a line minimisation in the region (xa,fa) (xc,fc) to find an
 208      intermediate (xb,fb) satisifying fa > fb < fc.  Choose an initial
 209      xb based on parabolic interpolation */
 210
 211   intermediate_point (fdf, x, p, dir / pnorm, pg,
 212                       stepa, stepc, fa, fc, x1, dx1, gradient, &stepb, &fb);
 213
 214   if (stepb == 0.0)
 215     {
 216       return GSL_ENOPROG;
 217     }
 218
 219   minimize (fdf, x, p, dir / pnorm,
 220             stepa, stepb, stepc, fa, fb, fc, tol,
 221             x1, dx1, x2, dx, gradient, &(state->step), f, &g1norm);
 222
 223   gsl_vector_memcpy (x, x2);
 224
 225   /* Choose a new conjugate direction for the next step */
 226
 227   state->iter = (state->iter + 1) % x->size;
 228
 229   if (state->iter == 0)
 230     {
 231       gsl_vector_memcpy (p, gradient);
 232       state->pnorm = g1norm;
 233     }
 234   else
 235     {
 236       /* p' = g1 - beta * p */
 237
 238       double beta = -pow (g1norm / g0norm, 2.0);
 239       gsl_blas_dscal (-beta, p);
 240       gsl_blas_daxpy (1.0, gradient, p);
 241       state->pnorm = gsl_blas_dnrm2 (p);
 242     }
 243
 244   state->g0norm = g1norm;
 245   gsl_vector_memcpy (g0, gradient);
 246
 247 #ifdef DEBUG
 248   printf ("updated conjugate directions\n");
 249   printf ("p: ");
 250   gsl_vector_fprintf (stdout, p, "%g");
 251   printf ("g: ");
 252   gsl_vector_fprintf (stdout, gradient, "%g");
 253 #endif
 254
 255   return GSL_SUCCESS;
 256 }
 257
 258
 259
 260 static const gsl_multimin_fdfminimizer_type conjugate_fr_type = {
 261   "conjugate_fr",               /* name */
 262   sizeof (conjugate_fr_state_t),
 263   &conjugate_fr_alloc,
 264   &conjugate_fr_set,
 265   &conjugate_fr_iterate,
 266   &conjugate_fr_restart,
 267   &conjugate_fr_free
 268 };
 269
 270 const gsl_multimin_fdfminimizer_type
 271   * gsl_multimin_fdfminimizer_conjugate_fr = &conjugate_fr_type;