|
Packit |
67cb25 |
/* multiroots/hybridj.c
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Brian Gough
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is free software; you can redistribute it and/or modify
|
|
Packit |
67cb25 |
* it under the terms of the GNU General Public License as published by
|
|
Packit |
67cb25 |
* the Free Software Foundation; either version 3 of the License, or (at
|
|
Packit |
67cb25 |
* your option) any later version.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is distributed in the hope that it will be useful, but
|
|
Packit |
67cb25 |
* WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
Packit |
67cb25 |
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
Packit |
67cb25 |
* General Public License for more details.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* You should have received a copy of the GNU General Public License
|
|
Packit |
67cb25 |
* along with this program; if not, write to the Free Software
|
|
Packit |
67cb25 |
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include <config.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include <stddef.h>
|
|
Packit |
67cb25 |
#include <stdlib.h>
|
|
Packit |
67cb25 |
#include <stdio.h>
|
|
Packit |
67cb25 |
#include <math.h>
|
|
Packit |
67cb25 |
#include <float.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include <gsl/gsl_math.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_errno.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_multiroots.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_linalg.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include "dogleg.c"
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
typedef struct
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t iter;
|
|
Packit |
67cb25 |
size_t ncfail;
|
|
Packit |
67cb25 |
size_t ncsuc;
|
|
Packit |
67cb25 |
size_t nslow1;
|
|
Packit |
67cb25 |
size_t nslow2;
|
|
Packit |
67cb25 |
double fnorm;
|
|
Packit |
67cb25 |
double delta;
|
|
Packit |
67cb25 |
gsl_matrix *q;
|
|
Packit |
67cb25 |
gsl_matrix *r;
|
|
Packit |
67cb25 |
gsl_vector *tau;
|
|
Packit |
67cb25 |
gsl_vector *diag;
|
|
Packit |
67cb25 |
gsl_vector *qtf;
|
|
Packit |
67cb25 |
gsl_vector *newton;
|
|
Packit |
67cb25 |
gsl_vector *gradient;
|
|
Packit |
67cb25 |
gsl_vector *x_trial;
|
|
Packit |
67cb25 |
gsl_vector *f_trial;
|
|
Packit |
67cb25 |
gsl_vector *df;
|
|
Packit |
67cb25 |
gsl_vector *qtdf;
|
|
Packit |
67cb25 |
gsl_vector *rdx;
|
|
Packit |
67cb25 |
gsl_vector *w;
|
|
Packit |
67cb25 |
gsl_vector *v;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
hybridj_state_t;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static int hybridj_alloc (void *vstate, size_t n);
|
|
Packit |
67cb25 |
static int hybridj_set (void *vstate, gsl_multiroot_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx);
|
|
Packit |
67cb25 |
static int hybridsj_set (void *vstate, gsl_multiroot_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx);
|
|
Packit |
67cb25 |
static int hybridj_set_impl (void *vstate, gsl_multiroot_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx, int scale);
|
|
Packit |
67cb25 |
static int hybridj_iterate (void *vstate, gsl_multiroot_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx);
|
|
Packit |
67cb25 |
static void hybridj_free (void *vstate);
|
|
Packit |
67cb25 |
static int hybridj_iterate_impl (void *vstate, gsl_multiroot_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx, int scale);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static int
|
|
Packit |
67cb25 |
hybridj_alloc (void *vstate, size_t n)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
hybridj_state_t *state = (hybridj_state_t *) vstate;
|
|
Packit |
67cb25 |
gsl_matrix *q, *r;
|
|
Packit |
67cb25 |
gsl_vector *tau, *diag, *qtf, *newton, *gradient, *x_trial, *f_trial,
|
|
Packit |
67cb25 |
*df, *qtdf, *rdx, *w, *v;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
q = gsl_matrix_calloc (n, n);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (q == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("failed to allocate space for q", GSL_ENOMEM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
state->q = q;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
r = gsl_matrix_calloc (n, n);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (r == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix_free (q);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
GSL_ERROR ("failed to allocate space for r", GSL_ENOMEM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
state->r = r;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
tau = gsl_vector_calloc (n);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (tau == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix_free (q);
|
|
Packit |
67cb25 |
gsl_matrix_free (r);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
GSL_ERROR ("failed to allocate space for tau", GSL_ENOMEM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
state->tau = tau;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
diag = gsl_vector_calloc (n);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (diag == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix_free (q);
|
|
Packit |
67cb25 |
gsl_matrix_free (r);
|
|
Packit |
67cb25 |
gsl_vector_free (tau);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
GSL_ERROR ("failed to allocate space for diag", GSL_ENOMEM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
state->diag = diag;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
qtf = gsl_vector_calloc (n);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (qtf == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix_free (q);
|
|
Packit |
67cb25 |
gsl_matrix_free (r);
|
|
Packit |
67cb25 |
gsl_vector_free (tau);
|
|
Packit |
67cb25 |
gsl_vector_free (diag);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
GSL_ERROR ("failed to allocate space for qtf", GSL_ENOMEM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
state->qtf = qtf;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
newton = gsl_vector_calloc (n);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (newton == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix_free (q);
|
|
Packit |
67cb25 |
gsl_matrix_free (r);
|
|
Packit |
67cb25 |
gsl_vector_free (tau);
|
|
Packit |
67cb25 |
gsl_vector_free (diag);
|
|
Packit |
67cb25 |
gsl_vector_free (qtf);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
GSL_ERROR ("failed to allocate space for newton", GSL_ENOMEM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
state->newton = newton;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gradient = gsl_vector_calloc (n);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (gradient == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix_free (q);
|
|
Packit |
67cb25 |
gsl_matrix_free (r);
|
|
Packit |
67cb25 |
gsl_vector_free (tau);
|
|
Packit |
67cb25 |
gsl_vector_free (diag);
|
|
Packit |
67cb25 |
gsl_vector_free (qtf);
|
|
Packit |
67cb25 |
gsl_vector_free (newton);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
GSL_ERROR ("failed to allocate space for gradient", GSL_ENOMEM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
state->gradient = gradient;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
x_trial = gsl_vector_calloc (n);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (x_trial == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix_free (q);
|
|
Packit |
67cb25 |
gsl_matrix_free (r);
|
|
Packit |
67cb25 |
gsl_vector_free (tau);
|
|
Packit |
67cb25 |
gsl_vector_free (diag);
|
|
Packit |
67cb25 |
gsl_vector_free (qtf);
|
|
Packit |
67cb25 |
gsl_vector_free (newton);
|
|
Packit |
67cb25 |
gsl_vector_free (gradient);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
GSL_ERROR ("failed to allocate space for x_trial", GSL_ENOMEM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
state->x_trial = x_trial;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
f_trial = gsl_vector_calloc (n);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (f_trial == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix_free (q);
|
|
Packit |
67cb25 |
gsl_matrix_free (r);
|
|
Packit |
67cb25 |
gsl_vector_free (tau);
|
|
Packit |
67cb25 |
gsl_vector_free (diag);
|
|
Packit |
67cb25 |
gsl_vector_free (qtf);
|
|
Packit |
67cb25 |
gsl_vector_free (newton);
|
|
Packit |
67cb25 |
gsl_vector_free (gradient);
|
|
Packit |
67cb25 |
gsl_vector_free (x_trial);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
GSL_ERROR ("failed to allocate space for f_trial", GSL_ENOMEM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
state->f_trial = f_trial;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
df = gsl_vector_calloc (n);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (df == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix_free (q);
|
|
Packit |
67cb25 |
gsl_matrix_free (r);
|
|
Packit |
67cb25 |
gsl_vector_free (tau);
|
|
Packit |
67cb25 |
gsl_vector_free (diag);
|
|
Packit |
67cb25 |
gsl_vector_free (qtf);
|
|
Packit |
67cb25 |
gsl_vector_free (newton);
|
|
Packit |
67cb25 |
gsl_vector_free (gradient);
|
|
Packit |
67cb25 |
gsl_vector_free (x_trial);
|
|
Packit |
67cb25 |
gsl_vector_free (f_trial);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
GSL_ERROR ("failed to allocate space for df", GSL_ENOMEM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
state->df = df;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
qtdf = gsl_vector_calloc (n);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (qtdf == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix_free (q);
|
|
Packit |
67cb25 |
gsl_matrix_free (r);
|
|
Packit |
67cb25 |
gsl_vector_free (tau);
|
|
Packit |
67cb25 |
gsl_vector_free (diag);
|
|
Packit |
67cb25 |
gsl_vector_free (qtf);
|
|
Packit |
67cb25 |
gsl_vector_free (newton);
|
|
Packit |
67cb25 |
gsl_vector_free (gradient);
|
|
Packit |
67cb25 |
gsl_vector_free (x_trial);
|
|
Packit |
67cb25 |
gsl_vector_free (f_trial);
|
|
Packit |
67cb25 |
gsl_vector_free (df);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
GSL_ERROR ("failed to allocate space for qtdf", GSL_ENOMEM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
state->qtdf = qtdf;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
rdx = gsl_vector_calloc (n);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (rdx == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix_free (q);
|
|
Packit |
67cb25 |
gsl_matrix_free (r);
|
|
Packit |
67cb25 |
gsl_vector_free (tau);
|
|
Packit |
67cb25 |
gsl_vector_free (diag);
|
|
Packit |
67cb25 |
gsl_vector_free (qtf);
|
|
Packit |
67cb25 |
gsl_vector_free (newton);
|
|
Packit |
67cb25 |
gsl_vector_free (gradient);
|
|
Packit |
67cb25 |
gsl_vector_free (x_trial);
|
|
Packit |
67cb25 |
gsl_vector_free (f_trial);
|
|
Packit |
67cb25 |
gsl_vector_free (df);
|
|
Packit |
67cb25 |
gsl_vector_free (qtdf);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
GSL_ERROR ("failed to allocate space for rdx", GSL_ENOMEM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
state->rdx = rdx;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
w = gsl_vector_calloc (n);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (w == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix_free (q);
|
|
Packit |
67cb25 |
gsl_matrix_free (r);
|
|
Packit |
67cb25 |
gsl_vector_free (tau);
|
|
Packit |
67cb25 |
gsl_vector_free (diag);
|
|
Packit |
67cb25 |
gsl_vector_free (qtf);
|
|
Packit |
67cb25 |
gsl_vector_free (newton);
|
|
Packit |
67cb25 |
gsl_vector_free (gradient);
|
|
Packit |
67cb25 |
gsl_vector_free (x_trial);
|
|
Packit |
67cb25 |
gsl_vector_free (f_trial);
|
|
Packit |
67cb25 |
gsl_vector_free (df);
|
|
Packit |
67cb25 |
gsl_vector_free (qtdf);
|
|
Packit |
67cb25 |
gsl_vector_free (rdx);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
GSL_ERROR ("failed to allocate space for w", GSL_ENOMEM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
state->w = w;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
v = gsl_vector_calloc (n);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (v == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix_free (q);
|
|
Packit |
67cb25 |
gsl_matrix_free (r);
|
|
Packit |
67cb25 |
gsl_vector_free (tau);
|
|
Packit |
67cb25 |
gsl_vector_free (diag);
|
|
Packit |
67cb25 |
gsl_vector_free (qtf);
|
|
Packit |
67cb25 |
gsl_vector_free (newton);
|
|
Packit |
67cb25 |
gsl_vector_free (gradient);
|
|
Packit |
67cb25 |
gsl_vector_free (x_trial);
|
|
Packit |
67cb25 |
gsl_vector_free (f_trial);
|
|
Packit |
67cb25 |
gsl_vector_free (df);
|
|
Packit |
67cb25 |
gsl_vector_free (qtdf);
|
|
Packit |
67cb25 |
gsl_vector_free (rdx);
|
|
Packit |
67cb25 |
gsl_vector_free (w);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
GSL_ERROR ("failed to allocate space for v", GSL_ENOMEM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
state->v = v;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static int
|
|
Packit |
67cb25 |
hybridj_set (void *vstate, gsl_multiroot_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int status = hybridj_set_impl (vstate, fdf, x, f, J, dx, 0);
|
|
Packit |
67cb25 |
return status ;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static int
|
|
Packit |
67cb25 |
hybridsj_set (void *vstate, gsl_multiroot_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int status = hybridj_set_impl (vstate, fdf, x, f, J, dx, 1);
|
|
Packit |
67cb25 |
return status ;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static int
|
|
Packit |
67cb25 |
hybridj_set_impl (void *vstate, gsl_multiroot_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx, int scale)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
hybridj_state_t *state = (hybridj_state_t *) vstate;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix *q = state->q;
|
|
Packit |
67cb25 |
gsl_matrix *r = state->r;
|
|
Packit |
67cb25 |
gsl_vector *tau = state->tau;
|
|
Packit |
67cb25 |
gsl_vector *diag = state->diag;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
GSL_MULTIROOT_FN_EVAL_F_DF (fdf, x, f, J);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
state->iter = 1;
|
|
Packit |
67cb25 |
state->fnorm = enorm (f);
|
|
Packit |
67cb25 |
state->ncfail = 0;
|
|
Packit |
67cb25 |
state->ncsuc = 0;
|
|
Packit |
67cb25 |
state->nslow1 = 0;
|
|
Packit |
67cb25 |
state->nslow2 = 0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_set_all (dx, 0.0);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Store column norms in diag */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (scale)
|
|
Packit |
67cb25 |
compute_diag (J, diag);
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
gsl_vector_set_all (diag, 1.0);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Set delta to factor |D x| or to factor if |D x| is zero */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
state->delta = compute_delta (diag, x);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Factorize J into QR decomposition */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_linalg_QR_decomp (J, tau);
|
|
Packit |
67cb25 |
gsl_linalg_QR_unpack (J, tau, q, r);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static int
|
|
Packit |
67cb25 |
hybridj_iterate (void *vstate, gsl_multiroot_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int status = hybridj_iterate_impl (vstate, fdf, x, f, J, dx, 0);
|
|
Packit |
67cb25 |
return status;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static int
|
|
Packit |
67cb25 |
hybridsj_iterate (void *vstate, gsl_multiroot_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int status = hybridj_iterate_impl (vstate, fdf, x, f, J, dx, 1);
|
|
Packit |
67cb25 |
return status;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static int
|
|
Packit |
67cb25 |
hybridj_iterate_impl (void *vstate, gsl_multiroot_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx, int scale)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
hybridj_state_t *state = (hybridj_state_t *) vstate;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
const double fnorm = state->fnorm;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix *q = state->q;
|
|
Packit |
67cb25 |
gsl_matrix *r = state->r;
|
|
Packit |
67cb25 |
gsl_vector *tau = state->tau;
|
|
Packit |
67cb25 |
gsl_vector *diag = state->diag;
|
|
Packit |
67cb25 |
gsl_vector *qtf = state->qtf;
|
|
Packit |
67cb25 |
gsl_vector *x_trial = state->x_trial;
|
|
Packit |
67cb25 |
gsl_vector *f_trial = state->f_trial;
|
|
Packit |
67cb25 |
gsl_vector *df = state->df;
|
|
Packit |
67cb25 |
gsl_vector *qtdf = state->qtdf;
|
|
Packit |
67cb25 |
gsl_vector *rdx = state->rdx;
|
|
Packit |
67cb25 |
gsl_vector *w = state->w;
|
|
Packit |
67cb25 |
gsl_vector *v = state->v;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double prered, actred;
|
|
Packit |
67cb25 |
double pnorm, fnorm1, fnorm1p;
|
|
Packit |
67cb25 |
double ratio;
|
|
Packit |
67cb25 |
double p1 = 0.1, p5 = 0.5, p001 = 0.001, p0001 = 0.0001;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Compute qtf = Q^T f */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
compute_qtf (q, f, qtf);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Compute dogleg step */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
dogleg (r, qtf, diag, state->delta, state->newton, state->gradient, dx);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Take a trial step */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
compute_trial_step (x, dx, state->x_trial);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
pnorm = scaled_enorm (diag, dx);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (state->iter == 1)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (pnorm < state->delta)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
state->delta = pnorm;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Evaluate function at x + p */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int status = GSL_MULTIROOT_FN_EVAL_F (fdf, x_trial, f_trial);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (status != GSL_SUCCESS)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return GSL_EBADFUNC;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Set df = f_trial - f */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
compute_df (f_trial, f, df);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Compute the scaled actual reduction */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
fnorm1 = enorm (f_trial);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
actred = compute_actual_reduction (fnorm, fnorm1);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Compute rdx = R dx */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
compute_rdx (r, dx, rdx);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Compute the scaled predicted reduction phi1p = |Q^T f + R dx| */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
fnorm1p = enorm_sum (qtf, rdx);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
prered = compute_predicted_reduction (fnorm, fnorm1p);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Compute the ratio of the actual to predicted reduction */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (prered > 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
ratio = actred / prered;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
ratio = 0;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Update the step bound */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (ratio < p1)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
state->ncsuc = 0;
|
|
Packit |
67cb25 |
state->ncfail++;
|
|
Packit |
67cb25 |
state->delta *= p5;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
state->ncfail = 0;
|
|
Packit |
67cb25 |
state->ncsuc++;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (ratio >= p5 || state->ncsuc > 1)
|
|
Packit |
67cb25 |
state->delta = GSL_MAX (state->delta, pnorm / p5);
|
|
Packit |
67cb25 |
if (fabs (ratio - 1) <= p1)
|
|
Packit |
67cb25 |
state->delta = pnorm / p5;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Test for successful iteration */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (ratio >= p0001)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_vector_memcpy (x, x_trial);
|
|
Packit |
67cb25 |
gsl_vector_memcpy (f, f_trial);
|
|
Packit |
67cb25 |
state->fnorm = fnorm1;
|
|
Packit |
67cb25 |
state->iter++;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Determine the progress of the iteration */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
state->nslow1++;
|
|
Packit |
67cb25 |
if (actred >= p001)
|
|
Packit |
67cb25 |
state->nslow1 = 0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (actred >= p1)
|
|
Packit |
67cb25 |
state->nslow2 = 0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (state->ncfail == 2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int status = GSL_MULTIROOT_FN_EVAL_DF (fdf, x, J);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (status != GSL_SUCCESS)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return GSL_EBADFUNC;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
state->nslow2++;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (state->iter == 1)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (scale)
|
|
Packit |
67cb25 |
compute_diag (J, diag);
|
|
Packit |
67cb25 |
state->delta = compute_delta (diag, x);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (scale)
|
|
Packit |
67cb25 |
update_diag (J, diag);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Factorize J into QR decomposition */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_linalg_QR_decomp (J, tau);
|
|
Packit |
67cb25 |
gsl_linalg_QR_unpack (J, tau, q, r);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Compute qtdf = Q^T df, w = (Q^T df - R dx)/|dx|, v = D^2 dx/|dx| */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
compute_qtf (q, df, qtdf);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
compute_wv (qtdf, rdx, dx, diag, pnorm, w, v);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Rank-1 update of the jacobian Q'R' = Q(R + w v^T) */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_linalg_QR_update (q, r, w, v);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* No progress as measured by jacobian evaluations */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (state->nslow2 == 5)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return GSL_ENOPROGJ;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* No progress as measured by function evaluations */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (state->nslow1 == 10)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return GSL_ENOPROG;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static void
|
|
Packit |
67cb25 |
hybridj_free (void *vstate)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
hybridj_state_t *state = (hybridj_state_t *) vstate;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_free (state->v);
|
|
Packit |
67cb25 |
gsl_vector_free (state->w);
|
|
Packit |
67cb25 |
gsl_vector_free (state->rdx);
|
|
Packit |
67cb25 |
gsl_vector_free (state->qtdf);
|
|
Packit |
67cb25 |
gsl_vector_free (state->df);
|
|
Packit |
67cb25 |
gsl_vector_free (state->f_trial);
|
|
Packit |
67cb25 |
gsl_vector_free (state->x_trial);
|
|
Packit |
67cb25 |
gsl_vector_free (state->gradient);
|
|
Packit |
67cb25 |
gsl_vector_free (state->newton);
|
|
Packit |
67cb25 |
gsl_vector_free (state->qtf);
|
|
Packit |
67cb25 |
gsl_vector_free (state->diag);
|
|
Packit |
67cb25 |
gsl_vector_free (state->tau);
|
|
Packit |
67cb25 |
gsl_matrix_free (state->r);
|
|
Packit |
67cb25 |
gsl_matrix_free (state->q);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static const gsl_multiroot_fdfsolver_type hybridj_type =
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
"hybridj", /* name */
|
|
Packit |
67cb25 |
sizeof (hybridj_state_t),
|
|
Packit |
67cb25 |
&hybridj_alloc,
|
|
Packit |
67cb25 |
&hybridj_set,
|
|
Packit |
67cb25 |
&hybridj_iterate,
|
|
Packit |
67cb25 |
&hybridj_free
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static const gsl_multiroot_fdfsolver_type hybridsj_type =
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
"hybridsj", /* name */
|
|
Packit |
67cb25 |
sizeof (hybridj_state_t),
|
|
Packit |
67cb25 |
&hybridj_alloc,
|
|
Packit |
67cb25 |
&hybridsj_set,
|
|
Packit |
67cb25 |
&hybridsj_iterate,
|
|
Packit |
67cb25 |
&hybridj_free
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
const gsl_multiroot_fdfsolver_type *gsl_multiroot_fdfsolver_hybridj = &hybridj_type;
|
|
Packit |
67cb25 |
const gsl_multiroot_fdfsolver_type *gsl_multiroot_fdfsolver_hybridsj = &hybridsj_type;
|