123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457 |
- /* multiroots/broyden.c
- *
- * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Brian Gough
- *
- * This program is free software; you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation; either version 3 of the License, or (at
- * your option) any later version.
- *
- * This program is distributed in the hope that it will be useful, but
- * WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
- * General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program; if not, write to the Free Software
- * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
- */
- #include "gsl__config.h"
- #include <stddef.h>
- #include <stdlib.h>
- #include <stdio.h>
- #include <math.h>
- #include <float.h>
- #include "gsl_math.h"
- #include "gsl_errno.h"
- #include "gsl_multiroots.h"
- #include "gsl_linalg.h"
- #include "gsl_multiroots__enorm.c"
- /* Broyden's method. It is not an efficient or modern algorithm but
- gives an example of a rank-1 update.
- C.G. Broyden, "A Class of Methods for Solving Nonlinear
- Simultaneous Equations", Mathematics of Computation, vol 19 (1965),
- p 577-593
- */
- typedef struct
- {
- gsl_matrix *H;
- gsl_matrix *lu;
- gsl_permutation *permutation;
- gsl_vector *v;
- gsl_vector *w;
- gsl_vector *y;
- gsl_vector *p;
- gsl_vector *fnew;
- gsl_vector *x_trial;
- double phi;
- }
- broyden_state_t;
- static int broyden_alloc (void *vstate, size_t n);
- static int broyden_set (void *vstate, gsl_multiroot_function * function, gsl_vector * x, gsl_vector * f, gsl_vector * dx);
- static int broyden_iterate (void *vstate, gsl_multiroot_function * function, gsl_vector * x, gsl_vector * f, gsl_vector * dx);
- static void broyden_free (void *vstate);
- static int
- broyden_alloc (void *vstate, size_t n)
- {
- broyden_state_t *state = (broyden_state_t *) vstate;
- gsl_vector *v, *w, *y, *fnew, *x_trial, *p;
- gsl_permutation *perm;
- gsl_matrix *m, *H;
- m = gsl_matrix_calloc (n, n);
- if (m == 0)
- {
- GSL_ERROR ("failed to allocate space for lu", GSL_ENOMEM);
- }
- state->lu = m;
- perm = gsl_permutation_calloc (n);
- if (perm == 0)
- {
- gsl_matrix_free (m);
- GSL_ERROR ("failed to allocate space for permutation", GSL_ENOMEM);
- }
- state->permutation = perm;
- H = gsl_matrix_calloc (n, n);
- if (H == 0)
- {
- gsl_permutation_free (perm);
- gsl_matrix_free (m);
- GSL_ERROR ("failed to allocate space for d", GSL_ENOMEM);
- }
- state->H = H;
- v = gsl_vector_calloc (n);
- if (v == 0)
- {
- gsl_matrix_free (H);
- gsl_permutation_free (perm);
- gsl_matrix_free (m);
- GSL_ERROR ("failed to allocate space for v", GSL_ENOMEM);
- }
- state->v = v;
- w = gsl_vector_calloc (n);
- if (w == 0)
- {
- gsl_vector_free (v);
- gsl_matrix_free (H);
- gsl_permutation_free (perm);
- gsl_matrix_free (m);
- GSL_ERROR ("failed to allocate space for w", GSL_ENOMEM);
- }
- state->w = w;
- y = gsl_vector_calloc (n);
- if (y == 0)
- {
- gsl_vector_free (w);
- gsl_vector_free (v);
- gsl_matrix_free (H);
- gsl_permutation_free (perm);
- gsl_matrix_free (m);
- GSL_ERROR ("failed to allocate space for y", GSL_ENOMEM);
- }
- state->y = y;
- fnew = gsl_vector_calloc (n);
- if (fnew == 0)
- {
- gsl_vector_free (y);
- gsl_vector_free (w);
- gsl_vector_free (v);
- gsl_matrix_free (H);
- gsl_permutation_free (perm);
- gsl_matrix_free (m);
- GSL_ERROR ("failed to allocate space for fnew", GSL_ENOMEM);
- }
- state->fnew = fnew;
- x_trial = gsl_vector_calloc (n);
- if (x_trial == 0)
- {
- gsl_vector_free (fnew);
- gsl_vector_free (y);
- gsl_vector_free (w);
- gsl_vector_free (v);
- gsl_matrix_free (H);
- gsl_permutation_free (perm);
- gsl_matrix_free (m);
- GSL_ERROR ("failed to allocate space for x_trial", GSL_ENOMEM);
- }
- state->x_trial = x_trial;
- p = gsl_vector_calloc (n);
- if (p == 0)
- {
- gsl_vector_free (x_trial);
- gsl_vector_free (fnew);
- gsl_vector_free (y);
- gsl_vector_free (w);
- gsl_vector_free (v);
- gsl_matrix_free (H);
- gsl_permutation_free (perm);
- gsl_matrix_free (m);
- GSL_ERROR ("failed to allocate space for p", GSL_ENOMEM);
- }
- state->p = p;
- return GSL_SUCCESS;
- }
- static int
- broyden_set (void *vstate, gsl_multiroot_function * function, gsl_vector * x, gsl_vector * f, gsl_vector * dx)
- {
- broyden_state_t *state = (broyden_state_t *) vstate;
- size_t i, j, n = function->n;
- int signum = 0;
- GSL_MULTIROOT_FN_EVAL (function, x, f);
- gsl_multiroot_fdjacobian (function, x, f, GSL_SQRT_DBL_EPSILON, state->lu);
- gsl_linalg_LU_decomp (state->lu, state->permutation, &signum);
- gsl_linalg_LU_invert (state->lu, state->permutation, state->H);
- for (i = 0; i < n; i++)
- for (j = 0; j < n; j++)
- gsl_matrix_set(state->H,i,j,-gsl_matrix_get(state->H,i,j));
- for (i = 0; i < n; i++)
- {
- gsl_vector_set (dx, i, 0.0);
- }
- state->phi = enorm (f);
- return GSL_SUCCESS;
- }
- static int
- broyden_iterate (void *vstate, gsl_multiroot_function * function, gsl_vector * x, gsl_vector * f, gsl_vector * dx)
- {
- broyden_state_t *state = (broyden_state_t *) vstate;
- double phi0, phi1, t, lambda;
- gsl_matrix *H = state->H;
- gsl_vector *p = state->p;
- gsl_vector *y = state->y;
- gsl_vector *v = state->v;
- gsl_vector *w = state->w;
- gsl_vector *fnew = state->fnew;
- gsl_vector *x_trial = state->x_trial;
- gsl_matrix *lu = state->lu;
- gsl_permutation *perm = state->permutation;
- size_t i, j, iter;
- size_t n = function->n;
- /* p = H f */
- for (i = 0; i < n; i++)
- {
- double sum = 0;
- for (j = 0; j < n; j++)
- {
- sum += gsl_matrix_get (H, i, j) * gsl_vector_get (f, j);
- }
- gsl_vector_set (p, i, sum);
- }
- t = 1;
- iter = 0;
- phi0 = state->phi;
- new_step:
- for (i = 0; i < n; i++)
- {
- double pi = gsl_vector_get (p, i);
- double xi = gsl_vector_get (x, i);
- gsl_vector_set (x_trial, i, xi + t * pi);
- }
- {
- int status = GSL_MULTIROOT_FN_EVAL (function, x_trial, fnew);
- if (status != GSL_SUCCESS)
- {
- return GSL_EBADFUNC;
- }
- }
- phi1 = enorm (fnew);
- iter++ ;
- if (phi1 > phi0 && iter < 10 && t > 0.1)
- {
- /* full step goes uphill, take a reduced step instead */
-
- double theta = phi1 / phi0;
- t *= (sqrt (1.0 + 6.0 * theta) - 1.0) / (3.0 * theta);
- goto new_step;
- }
- if (phi1 > phi0)
- {
- /* need to recompute Jacobian */
- int signum = 0;
-
- gsl_multiroot_fdjacobian (function, x, f, GSL_SQRT_DBL_EPSILON, lu);
-
- for (i = 0; i < n; i++)
- for (j = 0; j < n; j++)
- gsl_matrix_set(lu,i,j,-gsl_matrix_get(lu,i,j));
-
- gsl_linalg_LU_decomp (lu, perm, &signum);
- gsl_linalg_LU_invert (lu, perm, H);
-
- gsl_linalg_LU_solve (lu, perm, f, p);
- t = 1;
- for (i = 0; i < n; i++)
- {
- double pi = gsl_vector_get (p, i);
- double xi = gsl_vector_get (x, i);
- gsl_vector_set (x_trial, i, xi + t * pi);
- }
-
- {
- int status = GSL_MULTIROOT_FN_EVAL (function, x_trial, fnew);
-
- if (status != GSL_SUCCESS)
- {
- return GSL_EBADFUNC;
- }
- }
-
- phi1 = enorm (fnew);
- }
-
- /* y = f' - f */
- for (i = 0; i < n; i++)
- {
- double yi = gsl_vector_get (fnew, i) - gsl_vector_get (f, i);
- gsl_vector_set (y, i, yi);
- }
- /* v = H y */
- for (i = 0; i < n; i++)
- {
- double sum = 0;
- for (j = 0; j < n; j++)
- {
- sum += gsl_matrix_get (H, i, j) * gsl_vector_get (y, j);
- }
- gsl_vector_set (v, i, sum);
- }
- /* lambda = p . v */
- lambda = 0;
- for (i = 0; i < n; i++)
- {
- lambda += gsl_vector_get (p, i) * gsl_vector_get (v, i);
- }
- if (lambda == 0)
- {
- GSL_ERROR ("approximation to Jacobian has collapsed", GSL_EZERODIV) ;
- }
- /* v' = v + t * p */
- for (i = 0; i < n; i++)
- {
- double vi = gsl_vector_get (v, i) + t * gsl_vector_get (p, i);
- gsl_vector_set (v, i, vi);
- }
- /* w^T = p^T H */
- for (i = 0; i < n; i++)
- {
- double sum = 0;
- for (j = 0; j < n; j++)
- {
- sum += gsl_matrix_get (H, j, i) * gsl_vector_get (p, j);
- }
- gsl_vector_set (w, i, sum);
- }
- /* Hij -> Hij - (vi wj / lambda) */
- for (i = 0; i < n; i++)
- {
- double vi = gsl_vector_get (v, i);
- for (j = 0; j < n; j++)
- {
- double wj = gsl_vector_get (w, j);
- double Hij = gsl_matrix_get (H, i, j) - vi * wj / lambda;
- gsl_matrix_set (H, i, j, Hij);
- }
- }
- /* copy fnew into f */
- gsl_vector_memcpy (f, fnew);
- /* copy x_trial into x */
- gsl_vector_memcpy (x, x_trial);
- for (i = 0; i < n; i++)
- {
- double pi = gsl_vector_get (p, i);
- gsl_vector_set (dx, i, t * pi);
- }
- state->phi = phi1;
- return GSL_SUCCESS;
- }
- static void
- broyden_free (void *vstate)
- {
- broyden_state_t *state = (broyden_state_t *) vstate;
- gsl_matrix_free (state->H);
- gsl_matrix_free (state->lu);
- gsl_permutation_free (state->permutation);
-
- gsl_vector_free (state->v);
- gsl_vector_free (state->w);
- gsl_vector_free (state->y);
- gsl_vector_free (state->p);
- gsl_vector_free (state->fnew);
- gsl_vector_free (state->x_trial);
-
- }
- static const gsl_multiroot_fsolver_type broyden_type =
- {"broyden", /* name */
- sizeof (broyden_state_t),
- &broyden_alloc,
- &broyden_set,
- &broyden_iterate,
- &broyden_free};
- const gsl_multiroot_fsolver_type *gsl_multiroot_fsolver_broyden = &broyden_type;
|