summaryrefslogtreecommitdiff
path: root/gsl-1.9/multifit/lmder.c
diff options
context:
space:
mode:
Diffstat (limited to 'gsl-1.9/multifit/lmder.c')
-rw-r--r--gsl-1.9/multifit/lmder.c387
1 files changed, 387 insertions, 0 deletions
diff --git a/gsl-1.9/multifit/lmder.c b/gsl-1.9/multifit/lmder.c
new file mode 100644
index 0000000..0dbeca8
--- /dev/null
+++ b/gsl-1.9/multifit/lmder.c
@@ -0,0 +1,387 @@
+/* multfit/lmder.c
+ *
+ * Copyright (C) 1996, 1997, 1998, 1999, 2000 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 2 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 <config.h>
+
+#include <stddef.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <float.h>
+
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_errno.h>
+#include <gsl/gsl_multifit_nlin.h>
+#include <gsl/gsl_blas.h>
+#include <gsl/gsl_linalg.h>
+#include <gsl/gsl_permutation.h>
+
+
+typedef struct
+ {
+ size_t iter;
+ double xnorm;
+ double fnorm;
+ double delta;
+ double par;
+ gsl_matrix *r;
+ gsl_vector *tau;
+ gsl_vector *diag;
+ gsl_vector *qtf;
+ gsl_vector *newton;
+ gsl_vector *gradient;
+ gsl_vector *x_trial;
+ gsl_vector *f_trial;
+ gsl_vector *df;
+ gsl_vector *sdiag;
+ gsl_vector *rptdx;
+ gsl_vector *w;
+ gsl_vector *work1;
+ gsl_permutation * perm;
+ }
+lmder_state_t;
+
+static int lmder_alloc (void *vstate, size_t n, size_t p);
+static int lmder_set (void *vstate, gsl_multifit_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx);
+static int lmsder_set (void *vstate, gsl_multifit_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx);
+static int set (void *vstate, gsl_multifit_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx, int scale);
+static int lmder_iterate (void *vstate, gsl_multifit_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx);
+static void lmder_free (void *vstate);
+static int iterate (void *vstate, gsl_multifit_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx, int scale);
+
+#include "lmutil.c"
+#include "lmpar.c"
+#include "lmset.c"
+#include "lmiterate.c"
+
+
+static int
+lmder_alloc (void *vstate, size_t n, size_t p)
+{
+ lmder_state_t *state = (lmder_state_t *) vstate;
+ gsl_matrix *r;
+ gsl_vector *tau, *diag, *qtf, *newton, *gradient, *x_trial, *f_trial,
+ *df, *sdiag, *rptdx, *w, *work1;
+ gsl_permutation *perm;
+
+ r = gsl_matrix_calloc (n, p);
+
+ if (r == 0)
+ {
+ GSL_ERROR ("failed to allocate space for r", GSL_ENOMEM);
+ }
+
+ state->r = r;
+
+ tau = gsl_vector_calloc (GSL_MIN(n, p));
+
+ if (tau == 0)
+ {
+ gsl_matrix_free (r);
+
+ GSL_ERROR ("failed to allocate space for tau", GSL_ENOMEM);
+ }
+
+ state->tau = tau;
+
+ diag = gsl_vector_calloc (p);
+
+ if (diag == 0)
+ {
+ gsl_matrix_free (r);
+ gsl_vector_free (tau);
+
+ GSL_ERROR ("failed to allocate space for diag", GSL_ENOMEM);
+ }
+
+ state->diag = diag;
+
+ qtf = gsl_vector_calloc (n);
+
+ if (qtf == 0)
+ {
+ gsl_matrix_free (r);
+ gsl_vector_free (tau);
+ gsl_vector_free (diag);
+
+ GSL_ERROR ("failed to allocate space for qtf", GSL_ENOMEM);
+ }
+
+ state->qtf = qtf;
+
+ newton = gsl_vector_calloc (p);
+
+ if (newton == 0)
+ {
+ gsl_matrix_free (r);
+ gsl_vector_free (tau);
+ gsl_vector_free (diag);
+ gsl_vector_free (qtf);
+
+ GSL_ERROR ("failed to allocate space for newton", GSL_ENOMEM);
+ }
+
+ state->newton = newton;
+
+ gradient = gsl_vector_calloc (p);
+
+ if (gradient == 0)
+ {
+ gsl_matrix_free (r);
+ gsl_vector_free (tau);
+ gsl_vector_free (diag);
+ gsl_vector_free (qtf);
+ gsl_vector_free (newton);
+
+ GSL_ERROR ("failed to allocate space for gradient", GSL_ENOMEM);
+ }
+
+ state->gradient = gradient;
+
+ x_trial = gsl_vector_calloc (p);
+
+ if (x_trial == 0)
+ {
+ gsl_matrix_free (r);
+ gsl_vector_free (tau);
+ gsl_vector_free (diag);
+ gsl_vector_free (qtf);
+ gsl_vector_free (newton);
+ gsl_vector_free (gradient);
+
+ GSL_ERROR ("failed to allocate space for x_trial", GSL_ENOMEM);
+ }
+
+ state->x_trial = x_trial;
+
+ f_trial = gsl_vector_calloc (n);
+
+ if (f_trial == 0)
+ {
+ gsl_matrix_free (r);
+ gsl_vector_free (tau);
+ gsl_vector_free (diag);
+ gsl_vector_free (qtf);
+ gsl_vector_free (newton);
+ gsl_vector_free (gradient);
+ gsl_vector_free (x_trial);
+
+ GSL_ERROR ("failed to allocate space for f_trial", GSL_ENOMEM);
+ }
+
+ state->f_trial = f_trial;
+
+ df = gsl_vector_calloc (n);
+
+ if (df == 0)
+ {
+ gsl_matrix_free (r);
+ gsl_vector_free (tau);
+ gsl_vector_free (diag);
+ gsl_vector_free (qtf);
+ gsl_vector_free (newton);
+ gsl_vector_free (gradient);
+ gsl_vector_free (x_trial);
+ gsl_vector_free (f_trial);
+
+ GSL_ERROR ("failed to allocate space for df", GSL_ENOMEM);
+ }
+
+ state->df = df;
+
+ sdiag = gsl_vector_calloc (p);
+
+ if (sdiag == 0)
+ {
+ gsl_matrix_free (r);
+ gsl_vector_free (tau);
+ gsl_vector_free (diag);
+ gsl_vector_free (qtf);
+ gsl_vector_free (newton);
+ gsl_vector_free (gradient);
+ gsl_vector_free (x_trial);
+ gsl_vector_free (f_trial);
+ gsl_vector_free (df);
+
+ GSL_ERROR ("failed to allocate space for sdiag", GSL_ENOMEM);
+ }
+
+ state->sdiag = sdiag;
+
+
+ rptdx = gsl_vector_calloc (n);
+
+ if (rptdx == 0)
+ {
+ gsl_matrix_free (r);
+ gsl_vector_free (tau);
+ gsl_vector_free (diag);
+ gsl_vector_free (qtf);
+ gsl_vector_free (newton);
+ gsl_vector_free (gradient);
+ gsl_vector_free (x_trial);
+ gsl_vector_free (f_trial);
+ gsl_vector_free (df);
+ gsl_vector_free (sdiag);
+
+ GSL_ERROR ("failed to allocate space for rptdx", GSL_ENOMEM);
+ }
+
+ state->rptdx = rptdx;
+
+ w = gsl_vector_calloc (n);
+
+ if (w == 0)
+ {
+ gsl_matrix_free (r);
+ gsl_vector_free (tau);
+ gsl_vector_free (diag);
+ gsl_vector_free (qtf);
+ gsl_vector_free (newton);
+ gsl_vector_free (gradient);
+ gsl_vector_free (x_trial);
+ gsl_vector_free (f_trial);
+ gsl_vector_free (df);
+ gsl_vector_free (sdiag);
+ gsl_vector_free (rptdx);
+
+ GSL_ERROR ("failed to allocate space for w", GSL_ENOMEM);
+ }
+
+ state->w = w;
+
+ work1 = gsl_vector_calloc (p);
+
+ if (work1 == 0)
+ {
+ gsl_matrix_free (r);
+ gsl_vector_free (tau);
+ gsl_vector_free (diag);
+ gsl_vector_free (qtf);
+ gsl_vector_free (newton);
+ gsl_vector_free (gradient);
+ gsl_vector_free (x_trial);
+ gsl_vector_free (f_trial);
+ gsl_vector_free (df);
+ gsl_vector_free (sdiag);
+ gsl_vector_free (rptdx);
+ gsl_vector_free (w);
+
+ GSL_ERROR ("failed to allocate space for work1", GSL_ENOMEM);
+ }
+
+ state->work1 = work1;
+
+ perm = gsl_permutation_calloc (p);
+
+ if (perm == 0)
+ {
+ gsl_matrix_free (r);
+ gsl_vector_free (tau);
+ gsl_vector_free (diag);
+ gsl_vector_free (qtf);
+ gsl_vector_free (newton);
+ gsl_vector_free (gradient);
+ gsl_vector_free (x_trial);
+ gsl_vector_free (f_trial);
+ gsl_vector_free (df);
+ gsl_vector_free (sdiag);
+ gsl_vector_free (rptdx);
+ gsl_vector_free (w);
+ gsl_vector_free (work1);
+
+ GSL_ERROR ("failed to allocate space for perm", GSL_ENOMEM);
+ }
+
+ state->perm = perm;
+
+ return GSL_SUCCESS;
+}
+
+static int
+lmder_set (void *vstate, gsl_multifit_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx)
+{
+ int status = set (vstate, fdf, x, f, J, dx, 0);
+ return status ;
+}
+
+static int
+lmsder_set (void *vstate, gsl_multifit_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx)
+{
+ int status = set (vstate, fdf, x, f, J, dx, 1);
+ return status ;
+}
+
+static int
+lmder_iterate (void *vstate, gsl_multifit_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx)
+{
+ int status = iterate (vstate, fdf, x, f, J, dx, 0);
+ return status;
+}
+
+static int
+lmsder_iterate (void *vstate, gsl_multifit_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx)
+{
+ int status = iterate (vstate, fdf, x, f, J, dx, 1);
+ return status;
+}
+
+static void
+lmder_free (void *vstate)
+{
+ lmder_state_t *state = (lmder_state_t *) vstate;
+
+ gsl_permutation_free (state->perm);
+ gsl_vector_free (state->work1);
+ gsl_vector_free (state->w);
+ gsl_vector_free (state->rptdx);
+ gsl_vector_free (state->sdiag);
+ gsl_vector_free (state->df);
+ gsl_vector_free (state->f_trial);
+ gsl_vector_free (state->x_trial);
+ gsl_vector_free (state->gradient);
+ gsl_vector_free (state->newton);
+ gsl_vector_free (state->qtf);
+ gsl_vector_free (state->diag);
+ gsl_vector_free (state->tau);
+ gsl_matrix_free (state->r);
+}
+
+static const gsl_multifit_fdfsolver_type lmder_type =
+{
+ "lmder", /* name */
+ sizeof (lmder_state_t),
+ &lmder_alloc,
+ &lmder_set,
+ &lmder_iterate,
+ &lmder_free
+};
+
+static const gsl_multifit_fdfsolver_type lmsder_type =
+{
+ "lmsder", /* name */
+ sizeof (lmder_state_t),
+ &lmder_alloc,
+ &lmsder_set,
+ &lmsder_iterate,
+ &lmder_free
+};
+
+const gsl_multifit_fdfsolver_type *gsl_multifit_fdfsolver_lmder = &lmder_type;
+const gsl_multifit_fdfsolver_type *gsl_multifit_fdfsolver_lmsder = &lmsder_type;
generated by cgit v1.2.3 (git 2.25.1) at 2024-05-01 10:47:23 +0000