htswanalysis/MACS/lib/gsl/gsl-1.11/multifit/lmder.c

   1 /* multfit/lmder.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_multifit_nlin.h>
  31 #include <gsl/gsl_blas.h>
  32 #include <gsl/gsl_linalg.h>
  33 #include <gsl/gsl_permutation.h>
  34
  35
  36 typedef struct
  37   {
  38     size_t iter;
  39     double xnorm;
  40     double fnorm;
  41     double delta;
  42     double par;
  43     gsl_matrix *r;
  44     gsl_vector *tau;
  45     gsl_vector *diag;
  46     gsl_vector *qtf;
  47     gsl_vector *newton;
  48     gsl_vector *gradient;
  49     gsl_vector *x_trial;
  50     gsl_vector *f_trial;
  51     gsl_vector *df;
  52     gsl_vector *sdiag;
  53     gsl_vector *rptdx;
  54     gsl_vector *w;
  55     gsl_vector *work1;
  56     gsl_permutation * perm;
  57   }
  58 lmder_state_t;
  59
  60 static int lmder_alloc (void *vstate, size_t n, size_t p);
  61 static int lmder_set (void *vstate, gsl_multifit_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx);
  62 static int lmsder_set (void *vstate, gsl_multifit_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx);
  63 static int set (void *vstate, gsl_multifit_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx, int scale);
  64 static int lmder_iterate (void *vstate, gsl_multifit_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx);
  65 static void lmder_free (void *vstate);
  66 static int iterate (void *vstate, gsl_multifit_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx, int scale);
  67
  68 #include "lmutil.c"
  69 #include "lmpar.c"
  70 #include "lmset.c"
  71 #include "lmiterate.c"
  72
  73
  74 static int
  75 lmder_alloc (void *vstate, size_t n, size_t p)
  76 {
  77   lmder_state_t *state = (lmder_state_t *) vstate;
  78   gsl_matrix *r;
  79   gsl_vector *tau, *diag, *qtf, *newton, *gradient, *x_trial, *f_trial,
  80    *df, *sdiag, *rptdx, *w, *work1;
  81   gsl_permutation *perm;
  82
  83   r = gsl_matrix_calloc (n, p);
  84
  85   if (r == 0)
  86     {
  87       GSL_ERROR ("failed to allocate space for r", GSL_ENOMEM);
  88     }
  89
  90   state->r = r;
  91
  92   tau = gsl_vector_calloc (GSL_MIN(n, p));
  93
  94   if (tau == 0)
  95     {
  96       gsl_matrix_free (r);
  97
  98       GSL_ERROR ("failed to allocate space for tau", GSL_ENOMEM);
  99     }
 100
 101   state->tau = tau;
 102
 103   diag = gsl_vector_calloc (p);
 104
 105   if (diag == 0)
 106     {
 107       gsl_matrix_free (r);
 108       gsl_vector_free (tau);
 109
 110       GSL_ERROR ("failed to allocate space for diag", GSL_ENOMEM);
 111     }
 112
 113   state->diag = diag;
 114
 115   qtf = gsl_vector_calloc (n);
 116
 117   if (qtf == 0)
 118     {
 119       gsl_matrix_free (r);
 120       gsl_vector_free (tau);
 121       gsl_vector_free (diag);
 122
 123       GSL_ERROR ("failed to allocate space for qtf", GSL_ENOMEM);
 124     }
 125
 126   state->qtf = qtf;
 127
 128   newton = gsl_vector_calloc (p);
 129
 130   if (newton == 0)
 131     {
 132       gsl_matrix_free (r);
 133       gsl_vector_free (tau);
 134       gsl_vector_free (diag);
 135       gsl_vector_free (qtf);
 136
 137       GSL_ERROR ("failed to allocate space for newton", GSL_ENOMEM);
 138     }
 139
 140   state->newton = newton;
 141
 142   gradient = gsl_vector_calloc (p);
 143
 144   if (gradient == 0)
 145     {
 146       gsl_matrix_free (r);
 147       gsl_vector_free (tau);
 148       gsl_vector_free (diag);
 149       gsl_vector_free (qtf);
 150       gsl_vector_free (newton);
 151
 152       GSL_ERROR ("failed to allocate space for gradient", GSL_ENOMEM);
 153     }
 154
 155   state->gradient = gradient;
 156
 157   x_trial = gsl_vector_calloc (p);
 158
 159   if (x_trial == 0)
 160     {
 161       gsl_matrix_free (r);
 162       gsl_vector_free (tau);
 163       gsl_vector_free (diag);
 164       gsl_vector_free (qtf);
 165       gsl_vector_free (newton);
 166       gsl_vector_free (gradient);
 167
 168       GSL_ERROR ("failed to allocate space for x_trial", GSL_ENOMEM);
 169     }
 170
 171   state->x_trial = x_trial;
 172
 173   f_trial = gsl_vector_calloc (n);
 174
 175   if (f_trial == 0)
 176     {
 177       gsl_matrix_free (r);
 178       gsl_vector_free (tau);
 179       gsl_vector_free (diag);
 180       gsl_vector_free (qtf);
 181       gsl_vector_free (newton);
 182       gsl_vector_free (gradient);
 183       gsl_vector_free (x_trial);
 184
 185       GSL_ERROR ("failed to allocate space for f_trial", GSL_ENOMEM);
 186     }
 187
 188   state->f_trial = f_trial;
 189
 190   df = gsl_vector_calloc (n);
 191
 192   if (df == 0)
 193     {
 194       gsl_matrix_free (r);
 195       gsl_vector_free (tau);
 196       gsl_vector_free (diag);
 197       gsl_vector_free (qtf);
 198       gsl_vector_free (newton);
 199       gsl_vector_free (gradient);
 200       gsl_vector_free (x_trial);
 201       gsl_vector_free (f_trial);
 202
 203       GSL_ERROR ("failed to allocate space for df", GSL_ENOMEM);
 204     }
 205
 206   state->df = df;
 207
 208   sdiag = gsl_vector_calloc (p);
 209
 210   if (sdiag == 0)
 211     {
 212       gsl_matrix_free (r);
 213       gsl_vector_free (tau);
 214       gsl_vector_free (diag);
 215       gsl_vector_free (qtf);
 216       gsl_vector_free (newton);
 217       gsl_vector_free (gradient);
 218       gsl_vector_free (x_trial);
 219       gsl_vector_free (f_trial);
 220       gsl_vector_free (df);
 221
 222       GSL_ERROR ("failed to allocate space for sdiag", GSL_ENOMEM);
 223     }
 224
 225   state->sdiag = sdiag;
 226
 227
 228   rptdx = gsl_vector_calloc (n);
 229
 230   if (rptdx == 0)
 231     {
 232       gsl_matrix_free (r);
 233       gsl_vector_free (tau);
 234       gsl_vector_free (diag);
 235       gsl_vector_free (qtf);
 236       gsl_vector_free (newton);
 237       gsl_vector_free (gradient);
 238       gsl_vector_free (x_trial);
 239       gsl_vector_free (f_trial);
 240       gsl_vector_free (df);
 241       gsl_vector_free (sdiag);
 242
 243       GSL_ERROR ("failed to allocate space for rptdx", GSL_ENOMEM);
 244     }
 245
 246   state->rptdx = rptdx;
 247
 248   w = gsl_vector_calloc (n);
 249
 250   if (w == 0)
 251     {
 252       gsl_matrix_free (r);
 253       gsl_vector_free (tau);
 254       gsl_vector_free (diag);
 255       gsl_vector_free (qtf);
 256       gsl_vector_free (newton);
 257       gsl_vector_free (gradient);
 258       gsl_vector_free (x_trial);
 259       gsl_vector_free (f_trial);
 260       gsl_vector_free (df);
 261       gsl_vector_free (sdiag);
 262       gsl_vector_free (rptdx);
 263
 264       GSL_ERROR ("failed to allocate space for w", GSL_ENOMEM);
 265     }
 266
 267   state->w = w;
 268
 269   work1 = gsl_vector_calloc (p);
 270
 271   if (work1 == 0)
 272     {
 273       gsl_matrix_free (r);
 274       gsl_vector_free (tau);
 275       gsl_vector_free (diag);
 276       gsl_vector_free (qtf);
 277       gsl_vector_free (newton);
 278       gsl_vector_free (gradient);
 279       gsl_vector_free (x_trial);
 280       gsl_vector_free (f_trial);
 281       gsl_vector_free (df);
 282       gsl_vector_free (sdiag);
 283       gsl_vector_free (rptdx);
 284       gsl_vector_free (w);
 285
 286       GSL_ERROR ("failed to allocate space for work1", GSL_ENOMEM);
 287     }
 288
 289   state->work1 = work1;
 290
 291   perm = gsl_permutation_calloc (p);
 292
 293   if (perm == 0)
 294     {
 295       gsl_matrix_free (r);
 296       gsl_vector_free (tau);
 297       gsl_vector_free (diag);
 298       gsl_vector_free (qtf);
 299       gsl_vector_free (newton);
 300       gsl_vector_free (gradient);
 301       gsl_vector_free (x_trial);
 302       gsl_vector_free (f_trial);
 303       gsl_vector_free (df);
 304       gsl_vector_free (sdiag);
 305       gsl_vector_free (rptdx);
 306       gsl_vector_free (w);
 307       gsl_vector_free (work1);
 308
 309       GSL_ERROR ("failed to allocate space for perm", GSL_ENOMEM);
 310     }
 311
 312   state->perm = perm;
 313
 314   return GSL_SUCCESS;
 315 }
 316
 317 static int
 318 lmder_set (void *vstate, gsl_multifit_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx)
 319 {
 320   int status = set (vstate, fdf, x, f, J, dx, 0);
 321   return status ;
 322 }
 323
 324 static int
 325 lmsder_set (void *vstate, gsl_multifit_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx)
 326 {
 327   int status = set (vstate, fdf, x, f, J, dx, 1);
 328   return status ;
 329 }
 330
 331 static int
 332 lmder_iterate (void *vstate, gsl_multifit_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx)
 333 {
 334   int status = iterate (vstate, fdf, x, f, J, dx, 0);
 335   return status;
 336 }
 337
 338 static int
 339 lmsder_iterate (void *vstate, gsl_multifit_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx)
 340 {
 341   int status = iterate (vstate, fdf, x, f, J, dx, 1);
 342   return status;
 343 }
 344
 345 static void
 346 lmder_free (void *vstate)
 347 {
 348   lmder_state_t *state = (lmder_state_t *) vstate;
 349
 350   gsl_permutation_free (state->perm);
 351   gsl_vector_free (state->work1);
 352   gsl_vector_free (state->w);
 353   gsl_vector_free (state->rptdx);
 354   gsl_vector_free (state->sdiag);
 355   gsl_vector_free (state->df);
 356   gsl_vector_free (state->f_trial);
 357   gsl_vector_free (state->x_trial);
 358   gsl_vector_free (state->gradient);
 359   gsl_vector_free (state->newton);
 360   gsl_vector_free (state->qtf);
 361   gsl_vector_free (state->diag);
 362   gsl_vector_free (state->tau);
 363   gsl_matrix_free (state->r);
 364 }
 365
 366 static const gsl_multifit_fdfsolver_type lmder_type =
 367 {
 368   "lmder",                      /* name */
 369   sizeof (lmder_state_t),
 370   &lmder_alloc,
 371   &lmder_set,
 372   &lmder_iterate,
 373   &lmder_free
 374 };
 375
 376 static const gsl_multifit_fdfsolver_type lmsder_type =
 377 {
 378   "lmsder",                     /* name */
 379   sizeof (lmder_state_t),
 380   &lmder_alloc,
 381   &lmsder_set,
 382   &lmsder_iterate,
 383   &lmder_free
 384 };
 385
 386 const gsl_multifit_fdfsolver_type *gsl_multifit_fdfsolver_lmder = &lmder_type;
 387 const gsl_multifit_fdfsolver_type *gsl_multifit_fdfsolver_lmsder = &lmsder_type;