From f2184d0549b231319995e7f95d64cf89be4f7745 Mon Sep 17 00:00:00 2001
From: Packit <packit>
Date: Sep 21 2020 09:03:27 +0000
Subject: Apply patch gsl-tol.patch


patch_name: gsl-tol.patch
present_in_specfile: true

---

diff --git a/multifit/test_nonlinear.c b/multifit/test_nonlinear.c
index d01f8f3..8c60edc 100644
--- a/multifit/test_nonlinear.c
+++ b/multifit/test_nonlinear.c
@@ -175,8 +175,8 @@ static test_fdf_problem *test_fdf_nist[] = {
 static void
 test_nonlinear(void)
 {
-  const double xtol = pow(GSL_DBL_EPSILON, 0.9);
-  const double gtol = pow(GSL_DBL_EPSILON, 0.9);
+  const double xtol = pow(GSL_DBL_EPSILON, 0.9) * 10.0;
+  const double gtol = pow(GSL_DBL_EPSILON, 0.9) * 10.0;
   const double ftol = 0.0;
   size_t i, j;
 
diff --git a/multifit/test_nonlinear.c.tol b/multifit/test_nonlinear.c.tol
new file mode 100644
index 0000000..d01f8f3
--- /dev/null
+++ b/multifit/test_nonlinear.c.tol
@@ -0,0 +1,580 @@
+/* multifit/test_nonlinear.c
+ * 
+ * Copyright (C) 2007, 2013, 2014 Brian Gough, Patrick Alken
+ * 
+ * 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.
+ */
+
+typedef struct
+{
+  const char *name;
+  double *x0;       /* initial parameters (size p) */
+  double *sigma;
+  double *epsrel;   /* relative tolerance for solution checking */
+  size_t ntries;
+  void (*checksol) (const double x[], const double sumsq,
+                    const double epsrel, const char *sname,
+                    const char *pname);
+  gsl_multifit_function_fdf *fdf;
+} test_fdf_problem;
+
+#include "test_bard.c"
+#include "test_beale.c"
+#include "test_biggs.c"
+#include "test_box.c"
+#include "test_boxbod.c"
+#include "test_brown1.c"
+#include "test_brown2.c"
+#include "test_brown3.c"
+#include "test_eckerle.c"
+#include "test_enso.c"
+#include "test_exp1.c"
+#include "test_gaussian.c"
+#include "test_hahn1.c"
+#include "test_helical.c"
+#include "test_jennrich.c"
+#include "test_kirby2.c"
+#include "test_kowalik.c"
+#include "test_lin1.c"
+#include "test_lin2.c"
+#include "test_lin3.c"
+#include "test_meyer.c"
+#include "test_meyerscal.c"
+#include "test_osborne.c"
+#include "test_penalty1.c"
+#include "test_penalty2.c"
+#include "test_powell1.c"
+#include "test_powell2.c"
+#include "test_powell3.c"
+#include "test_rat42.c"
+#include "test_rat43.c"
+#include "test_rosenbrock.c"
+#include "test_rosenbrocke.c"
+#include "test_roth.c"
+#include "test_thurber.c"
+#include "test_vardim.c"
+#include "test_watson.c"
+#include "test_wood.c"
+
+#include "test_wnlin.c"
+
+static void test_fdf(const gsl_multifit_fdfsolver_type * T,
+                     const double xtol, const double gtol,
+                     const double ftol, const double epsrel,
+                     const double x0_scale, test_fdf_problem *problem,
+                     const double *wts);
+static void test_fdfridge(const gsl_multifit_fdfsolver_type * T,
+                          const double xtol, const double gtol,
+                          const double ftol, const double epsrel,
+                          const double x0_scale, test_fdf_problem *problem,
+                          const double *wts);
+static void test_fdf_checksol(const char *sname, const char *pname,
+                              const double epsrel,
+                              gsl_multifit_fdfsolver *s,
+                              test_fdf_problem *problem);
+static void test_scale_x0(gsl_vector *x0, const double scale);
+
+/*
+ * These test problems are taken from
+ *
+ * H. B. Nielsen, UCTP test problems for unconstrained optimization,
+ * IMM Department of Mathematical Modeling, Tech. Report IMM-REP-2000-17,
+ * 2000.
+ */
+static test_fdf_problem *test_fdf_nielsen[] = {
+  &lin1_problem,       /* 1 */
+  &lin2_problem,       /* 2 */
+  &lin3_problem,       /* 3 */
+  &rosenbrock_problem, /* 4 */
+  &helical_problem,    /* 5 */
+  &powell1_problem,    /* 6 */
+  &roth_problem,       /* 7 */
+  &bard_problem,       /* 8 */
+  &kowalik_problem,    /* 9 */
+  &meyer_problem,      /* 10 */
+  &watson_problem,     /* 11 */
+  &box_problem,        /* 12 */
+  &jennrich_problem,   /* 13 */
+  &brown1_problem,     /* 14 */
+  &brown2_problem,     /* 16 */
+  &osborne_problem,    /* 17 */
+  &exp1_problem,       /* 18 */
+  &meyerscal_problem,  /* 20 */
+
+  &powell2_problem,
+
+  NULL
+};
+
+/*
+ * These tests are from
+ *
+ * J. J. More, B. S. Garbow and K. E. Hillstrom, Testing
+ * Unconstrained Optimization Software, ACM Trans. Math. Soft.
+ * Vol 7, No 1, 1981.
+ *
+ * Many of these overlap with the Nielsen tests
+ */
+static test_fdf_problem *test_fdf_more[] = {
+  &rosenbrock_problem,   /* 1 */
+  &roth_problem,         /* 2 */
+  &powell3_problem,      /* 3 */
+  &brown3_problem,       /* 4 */
+  &beale_problem,        /* 5 */
+  &jennrich_problem,     /* 6 */
+  &helical_problem,      /* 7 */
+  &bard_problem,         /* 8 */
+  &gaussian_problem,     /* 9 */
+  &meyer_problem,        /* 10 */
+  &box_problem,          /* 12 */
+  &powell1_problem,      /* 13 */
+  &wood_problem,         /* 14 */
+  &kowalik_problem,      /* 15 */
+  &brown1_problem,       /* 16 */
+  &osborne_problem,      /* 17 */
+  &biggs_problem,        /* 18 */
+  &watson_problem,       /* 20 */
+  &rosenbrocke_problem,  /* 21 */
+  &penalty1_problem,     /* 23 */
+  &penalty2_problem,     /* 24 */
+  &vardim_problem,       /* 25 */
+  &brown2_problem,       /* 27 */
+  &lin1_problem,         /* 32 */
+  &lin2_problem,         /* 33 */
+  &lin3_problem,         /* 34 */
+
+  NULL
+};
+
+/* NIST test cases */
+static test_fdf_problem *test_fdf_nist[] = {
+  &kirby2_problem,
+  &hahn1_problem,
+  &enso_problem,
+  &thurber_problem,
+  &boxbod_problem,
+  &rat42_problem,
+  &eckerle_problem,
+  &rat43_problem,
+
+  NULL
+};
+
+static void
+test_nonlinear(void)
+{
+  const double xtol = pow(GSL_DBL_EPSILON, 0.9);
+  const double gtol = pow(GSL_DBL_EPSILON, 0.9);
+  const double ftol = 0.0;
+  size_t i, j;
+
+  /* Nielsen tests */
+  for (i = 0; test_fdf_nielsen[i] != NULL; ++i)
+    {
+      test_fdf_problem *problem = test_fdf_nielsen[i];
+      double epsrel = *(problem->epsrel);
+      double scale = 1.0;
+
+      for (j = 0; j < problem->ntries; ++j)
+        {
+          double eps_scale = epsrel * scale;
+
+          test_fdf(gsl_multifit_fdfsolver_lmsder, xtol, gtol, ftol,
+                   eps_scale, scale, problem, NULL);
+          test_fdfridge(gsl_multifit_fdfsolver_lmsder, xtol, gtol, ftol,
+                        eps_scale, scale, problem, NULL);
+
+          /* test finite difference Jacobian */
+          {
+            gsl_multifit_function_fdf fdf;
+            fdf.df = problem->fdf->df;
+            problem->fdf->df = NULL;
+            test_fdf(gsl_multifit_fdfsolver_lmsder, xtol, gtol, ftol,
+                     1.0e5 * eps_scale, 1.0, problem, NULL);
+            test_fdfridge(gsl_multifit_fdfsolver_lmsder, xtol, gtol, ftol,
+                          1.0e5 * eps_scale, 1.0, problem, NULL);
+            problem->fdf->df = fdf.df;
+          }
+
+          scale *= 10.0;
+        }
+
+      test_fdf(gsl_multifit_fdfsolver_lmniel, xtol, gtol, ftol,
+               10.0 * epsrel, 1.0, problem, NULL);
+    }
+
+  /* More tests */
+  for (i = 0; test_fdf_more[i] != NULL; ++i)
+    {
+      test_fdf_problem *problem = test_fdf_more[i];
+      double epsrel = *(problem->epsrel);
+      double scale = 1.0;
+
+      for (j = 0; j < problem->ntries; ++j)
+        {
+          double eps_scale = epsrel * scale;
+
+          test_fdf(gsl_multifit_fdfsolver_lmsder, xtol, gtol, ftol,
+                   eps_scale, scale, problem, NULL);
+          test_fdfridge(gsl_multifit_fdfsolver_lmsder, xtol, gtol, ftol,
+                        eps_scale, scale, problem, NULL);
+
+          /* test finite difference Jacobian */
+          {
+            gsl_multifit_function_fdf fdf;
+            fdf.df = problem->fdf->df;
+            problem->fdf->df = NULL;
+            test_fdf(gsl_multifit_fdfsolver_lmsder, xtol, gtol, ftol,
+                     1.0e5 * eps_scale, 1.0, problem, NULL);
+            test_fdfridge(gsl_multifit_fdfsolver_lmsder, xtol, gtol, ftol,
+                          1.0e5 * eps_scale, 1.0, problem, NULL);
+            problem->fdf->df = fdf.df;
+          }
+
+          scale *= 10.0;
+        }
+
+      test_fdf(gsl_multifit_fdfsolver_lmniel, xtol, gtol, ftol,
+               10.0 * epsrel, 1.0, problem, NULL);
+    }
+
+  /* NIST tests */
+  for (i = 0; test_fdf_nist[i] != NULL; ++i)
+    {
+      test_fdf_problem *problem = test_fdf_nist[i];
+      double epsrel = *(problem->epsrel);
+      double scale = 1.0;
+
+      for (j = 0; j < problem->ntries; ++j)
+        {
+          double eps_scale = epsrel * scale;
+
+          test_fdf(gsl_multifit_fdfsolver_lmsder, xtol, gtol, ftol,
+                   eps_scale, scale, problem, NULL);
+          test_fdf(gsl_multifit_fdfsolver_lmder, xtol, gtol, ftol,
+                   eps_scale, scale, problem, NULL);
+
+          /* test finite difference Jacobian */
+          {
+            gsl_multifit_function_fdf fdf;
+            fdf.df = problem->fdf->df;
+            problem->fdf->df = NULL;
+            test_fdf(gsl_multifit_fdfsolver_lmsder, xtol, gtol, ftol,
+                     eps_scale, 1.0, problem, NULL);
+            test_fdf(gsl_multifit_fdfsolver_lmder, xtol, gtol, ftol,
+                     eps_scale, scale, problem, NULL);
+            problem->fdf->df = fdf.df;
+          }
+
+          scale *= 10.0;
+        }
+
+      test_fdf(gsl_multifit_fdfsolver_lmniel, xtol, gtol, ftol,
+               epsrel, 1.0, problem, NULL);
+    }
+
+  /* test weighted nonlinear least squares */
+
+  /* internal weighting in _f and _df functions */
+  test_fdf(gsl_multifit_fdfsolver_lmsder, xtol, gtol, ftol,
+           wnlin_epsrel, 1.0, &wnlin_problem1, NULL);
+  test_fdf(gsl_multifit_fdfsolver_lmniel, xtol, gtol, ftol,
+           wnlin_epsrel, 1.0, &wnlin_problem1, NULL);
+  test_fdfridge(gsl_multifit_fdfsolver_lmsder, xtol, gtol, ftol,
+                wnlin_epsrel, 1.0, &wnlin_problem1, NULL);
+  test_fdfridge(gsl_multifit_fdfsolver_lmniel, xtol, gtol, ftol,
+                wnlin_epsrel, 1.0, &wnlin_problem1, NULL);
+
+  /* weighting through fdfsolver_wset */
+  test_fdf(gsl_multifit_fdfsolver_lmsder, xtol, gtol, ftol,
+           wnlin_epsrel, 1.0, &wnlin_problem2, wnlin_W);
+  test_fdf(gsl_multifit_fdfsolver_lmniel, xtol, gtol, ftol,
+           wnlin_epsrel, 1.0, &wnlin_problem2, wnlin_W);
+  test_fdfridge(gsl_multifit_fdfsolver_lmsder, xtol, gtol, ftol,
+                wnlin_epsrel, 1.0, &wnlin_problem2, wnlin_W);
+  test_fdfridge(gsl_multifit_fdfsolver_lmniel, xtol, gtol, ftol,
+                wnlin_epsrel, 1.0, &wnlin_problem2, wnlin_W);
+}
+
+/*
+test_fdf()
+  Test a weighted nonlinear least squares problem
+
+Inputs: T        - solver to use
+        xtol     - tolerance in x
+        gtol     - tolerance in gradient
+        ftol     - tolerance in residual vector
+        epsrel   - relative error tolerance in solution
+        x0_scale - to test robustness against starting points,
+                   the standard starting point in 'problem' is
+                   multiplied by this scale factor:
+                   x0 <- x0 * x0_scale
+                   If x0 = 0, then all components of x0 are set to
+                   x0_scale
+        problem  - contains the nonlinear problem and solution point
+        wts      - weight vector (NULL for unweighted)
+*/
+
+static void
+test_fdf(const gsl_multifit_fdfsolver_type * T, const double xtol,
+         const double gtol, const double ftol,
+         const double epsrel, const double x0_scale,
+         test_fdf_problem *problem,
+         const double *wts)
+{
+  gsl_multifit_function_fdf *fdf = problem->fdf;
+  const size_t n = fdf->n;
+  const size_t p = fdf->p;
+  const size_t max_iter = 1500;
+  gsl_vector *x0 = gsl_vector_alloc(p);
+  gsl_vector_view x0v = gsl_vector_view_array(problem->x0, p);
+  gsl_multifit_fdfsolver *s = gsl_multifit_fdfsolver_alloc (T, n, p);
+  const char *pname = problem->name;
+  char sname[2048];
+  int status, info;
+
+  sprintf(sname, "%s/scale=%g%s",
+    gsl_multifit_fdfsolver_name(s), x0_scale,
+    problem->fdf->df ? "" : "/fdiff");
+
+  /* scale starting point x0 */
+  gsl_vector_memcpy(x0, &x0v.vector);
+  test_scale_x0(x0, x0_scale);
+
+  if (wts)
+    {
+      gsl_vector_const_view wv = gsl_vector_const_view_array(wts, n);
+      gsl_multifit_fdfsolver_wset(s, fdf, x0, &wv.vector);
+    }
+  else
+    gsl_multifit_fdfsolver_set(s, fdf, x0);
+
+  status = gsl_multifit_fdfsolver_driver(s, max_iter, xtol, gtol,
+                                         ftol, &info);
+  gsl_test(status, "%s/%s did not converge, status=%s",
+           sname, pname, gsl_strerror(status));
+
+  /* check solution */
+  test_fdf_checksol(sname, pname, epsrel, s, problem);
+
+  if (wts == NULL)
+    {
+      /* test again with weighting matrix W = I */
+      gsl_vector *wv = gsl_vector_alloc(n);
+
+      sprintf(sname, "%s/scale=%g%s/weights",
+        gsl_multifit_fdfsolver_name(s), x0_scale,
+        problem->fdf->df ? "" : "/fdiff");
+
+      gsl_vector_memcpy(x0, &x0v.vector);
+      test_scale_x0(x0, x0_scale);
+
+      gsl_vector_set_all(wv, 1.0);
+      gsl_multifit_fdfsolver_wset(s, fdf, x0, wv);
+  
+      status = gsl_multifit_fdfsolver_driver(s, max_iter, xtol, gtol,
+                                             ftol, &info);
+      gsl_test(status, "%s/%s did not converge, status=%s",
+               sname, pname, gsl_strerror(status));
+
+      test_fdf_checksol(sname, pname, epsrel, s, problem);
+
+      gsl_vector_free(wv);
+    }
+
+  gsl_multifit_fdfsolver_free(s);
+  gsl_vector_free(x0);
+}
+
+/*
+test_fdfridge()
+  Test a nonlinear least squares problem
+
+Inputs: T        - solver to use
+        xtol     - tolerance in x
+        gtol     - tolerance in gradient
+        ftol     - tolerance in residual vector
+        epsrel   - relative error tolerance in solution
+        x0_scale - to test robustness against starting points,
+                   the standard starting point in 'problem' is
+                   multiplied by this scale factor:
+                   x0 <- x0 * x0_scale
+                   If x0 = 0, then all components of x0 are set to
+                   x0_scale
+        problem  - contains the nonlinear problem and solution point
+        wts      - weight vector
+*/
+
+static void
+test_fdfridge(const gsl_multifit_fdfsolver_type * T, const double xtol,
+              const double gtol, const double ftol,
+              const double epsrel, const double x0_scale,
+              test_fdf_problem *problem, const double *wts)
+{
+  gsl_multifit_function_fdf *fdf = problem->fdf;
+  const size_t n = fdf->n;
+  const size_t p = fdf->p;
+  const size_t max_iter = 1500;
+  gsl_vector *x0 = gsl_vector_alloc(p);
+  gsl_vector_view x0v = gsl_vector_view_array(problem->x0, p);
+  gsl_multifit_fdfridge *w = gsl_multifit_fdfridge_alloc (T, n, p);
+  const char *pname = problem->name;
+  char sname[2048];
+  int status, info;
+  double lambda = 0.0;
+
+  sprintf(sname, "ridge/%s", gsl_multifit_fdfridge_name(w));
+
+  /* scale starting point x0 */
+  gsl_vector_memcpy(x0, &x0v.vector);
+  test_scale_x0(x0, x0_scale);
+
+  /* test undamped case with lambda = 0 */
+  if (wts)
+    {
+      gsl_vector_const_view wv = gsl_vector_const_view_array(wts, n);
+      gsl_multifit_fdfridge_wset(w, fdf, x0, lambda, &wv.vector);
+    }
+  else
+    gsl_multifit_fdfridge_set(w, fdf, x0, lambda);
+
+  status = gsl_multifit_fdfridge_driver(w, max_iter, xtol, gtol,
+                                        ftol, &info);
+  gsl_test(status, "%s/%s did not converge, status=%s",
+           sname, pname, gsl_strerror(status));
+
+  /* check solution */
+  test_fdf_checksol(sname, pname, epsrel, w->s, problem);
+
+  /* test for self consisent solution with L = \lambda I */
+  {
+    const double eps = 1.0e-10;
+    gsl_matrix *L = gsl_matrix_calloc(p, p);
+    gsl_vector_view diag = gsl_matrix_diagonal(L);
+    gsl_multifit_fdfridge *w2 = gsl_multifit_fdfridge_alloc (T, n, p);
+    gsl_vector *y0 = gsl_vector_alloc(p);
+    size_t i;
+
+    /* pick some value for lambda and set L = \lambda I */
+    lambda = 5.0;
+    gsl_vector_set_all(&diag.vector, lambda);
+
+    /* scale initial vector */
+    gsl_vector_memcpy(x0, &x0v.vector);
+    test_scale_x0(x0, x0_scale);
+    gsl_vector_memcpy(y0, x0);
+
+    if (wts)
+      {
+        gsl_vector_const_view wv = gsl_vector_const_view_array(wts, n);
+        gsl_multifit_fdfridge_wset(w, fdf, x0, lambda, &wv.vector);
+        gsl_multifit_fdfridge_wset3(w2, fdf, y0, L, &wv.vector);
+      }
+    else
+      {
+        gsl_multifit_fdfridge_set(w, fdf, x0, lambda);
+        gsl_multifit_fdfridge_set3(w2, fdf, y0, L);
+      }
+
+    /* solve with scalar lambda routine */
+    status = gsl_multifit_fdfridge_driver(w, max_iter, xtol, gtol,
+                                          ftol, &info);
+    gsl_test(status, "%s/lambda/%s did not converge, status=%s",
+             sname, pname, gsl_strerror(status));
+
+    /* solve with general matrix routine */
+    status = gsl_multifit_fdfridge_driver(w2, max_iter, xtol, gtol,
+                                          ftol, &info);
+    gsl_test(status, "%s/L/%s did not converge, status=%s",
+             sname, pname, gsl_strerror(status));
+
+    /* test x = y */
+    for (i = 0; i < p; ++i)
+      {
+        double xi = gsl_vector_get(w->s->x, i);
+        double yi = gsl_vector_get(w2->s->x, i);
+
+        if (fabs(xi) < eps)
+          {
+            gsl_test_abs(yi, xi, eps, "%s/%s ridge lambda=%g i=%zu",
+                         sname, pname, lambda, i);
+          }
+        else
+          {
+            gsl_test_rel(yi, xi, eps, "%s/%s ridge lambda=%g i=%zu",
+                         sname, pname, lambda, i);
+          }
+      }
+
+    gsl_matrix_free(L);
+    gsl_vector_free(y0);
+    gsl_multifit_fdfridge_free(w2);
+  }
+
+  gsl_multifit_fdfridge_free(w);
+  gsl_vector_free(x0);
+}
+
+static void
+test_fdf_checksol(const char *sname, const char *pname,
+                  const double epsrel, gsl_multifit_fdfsolver *s,
+                  test_fdf_problem *problem)
+{
+  gsl_multifit_function_fdf *fdf = problem->fdf;
+  const double *sigma = problem->sigma;
+  gsl_vector *f = gsl_multifit_fdfsolver_residual(s);
+  gsl_vector *x = gsl_multifit_fdfsolver_position(s);
+  double sumsq;
+
+  /* check solution vector x and sumsq = ||f||^2 */
+  gsl_blas_ddot(f, f, &sumsq);
+  (problem->checksol)(x->data, sumsq, epsrel, sname, pname);
+
+#if 1
+  /* check variances */
+  if (sigma)
+    {
+      const size_t n = fdf->n;
+      const size_t p = fdf->p;
+      size_t i;
+      gsl_matrix * J = gsl_matrix_alloc(n, p);
+      gsl_matrix * covar = gsl_matrix_alloc (p, p);
+
+      gsl_multifit_fdfsolver_jac (s, J);
+      gsl_multifit_covar(J, 0.0, covar);
+
+      for (i = 0; i < p; i++) 
+        {
+          double ei = sqrt(sumsq/(n-p))*sqrt(gsl_matrix_get(covar,i,i));
+          gsl_test_rel (ei, sigma[i], epsrel, 
+                        "%s/%s, sigma(%d)", sname, pname, i) ;
+        }
+
+      gsl_matrix_free (J);
+      gsl_matrix_free (covar);
+    }
+#endif
+}
+
+static void
+test_scale_x0(gsl_vector *x0, const double scale)
+{
+  double nx = gsl_blas_dnrm2(x0);
+
+  if (nx == 0.0)
+    gsl_vector_set_all(x0, scale);
+  else
+    gsl_vector_scale(x0, scale);
+} /* test_scale_x0() */
diff --git a/multifit/test_powell1.c b/multifit/test_powell1.c
index 572de3b..76dcc1a 100644
--- a/multifit/test_powell1.c
+++ b/multifit/test_powell1.c
@@ -4,7 +4,7 @@
 #define powell1_NTRIES   4
 
 static double powell1_x0[powell1_P] = { 3.0, -1.0, 0.0, 1.0 };
-static double powell1_epsrel = 1.0e-5;
+static double powell1_epsrel = 2.0e-5;
 
 static void
 powell1_checksol(const double x[], const double sumsq,
diff --git a/multifit/test_powell1.c.tol b/multifit/test_powell1.c.tol
new file mode 100644
index 0000000..572de3b
--- /dev/null
+++ b/multifit/test_powell1.c.tol
@@ -0,0 +1,97 @@
+#define powell1_N        4
+#define powell1_P        4
+
+#define powell1_NTRIES   4
+
+static double powell1_x0[powell1_P] = { 3.0, -1.0, 0.0, 1.0 };
+static double powell1_epsrel = 1.0e-5;
+
+static void
+powell1_checksol(const double x[], const double sumsq,
+                 const double epsrel, const char *sname,
+                 const char *pname)
+{
+  size_t i;
+  const double sumsq_exact = 0.0;
+
+  gsl_test_rel(sumsq, sumsq_exact, epsrel, "%s/%s sumsq",
+               sname, pname);
+
+  for (i = 0; i < powell1_P; ++i)
+    {
+      gsl_test_rel(x[i], 0.0, epsrel, "%s/%s i=%zu",
+                   sname, pname, i);
+    }
+}
+
+static int
+powell1_f (const gsl_vector * x, void *params, gsl_vector * f)
+{
+  double x1 = gsl_vector_get (x, 0);
+  double x2 = gsl_vector_get (x, 1);
+  double x3 = gsl_vector_get (x, 2);
+  double x4 = gsl_vector_get (x, 3);
+
+  gsl_vector_set(f, 0, x1 + 10.0*x2);
+  gsl_vector_set(f, 1, sqrt(5.0) * (x3 - x4));
+  gsl_vector_set(f, 2, pow(x2 - 2.0*x3, 2.0));
+  gsl_vector_set(f, 3, sqrt(10.0) * pow((x1 - x4), 2.0));
+
+  return GSL_SUCCESS;
+}
+
+static int
+powell1_df (const gsl_vector * x, void *params, gsl_matrix * J)
+{
+  double x1 = gsl_vector_get (x, 0);
+  double x2 = gsl_vector_get (x, 1);
+  double x3 = gsl_vector_get (x, 2);
+  double x4 = gsl_vector_get (x, 3);
+  double term1 = x2 - 2.0*x3;
+  double term2 = x1 - x4;
+
+  gsl_matrix_set(J, 0, 0, 1.0);
+  gsl_matrix_set(J, 0, 1, 10.0);
+  gsl_matrix_set(J, 0, 2, 0.0);
+  gsl_matrix_set(J, 0, 3, 0.0);
+
+  gsl_matrix_set(J, 1, 0, 0.0);
+  gsl_matrix_set(J, 1, 1, 0.0);
+  gsl_matrix_set(J, 1, 2, sqrt(5.0));
+  gsl_matrix_set(J, 1, 3, -sqrt(5.0));
+
+  gsl_matrix_set(J, 2, 0, 0.0);
+  gsl_matrix_set(J, 2, 1, 2.0*term1);
+  gsl_matrix_set(J, 2, 2, -4.0*term1);
+  gsl_matrix_set(J, 2, 3, 0.0);
+
+  gsl_matrix_set(J, 3, 0, 2.0*sqrt(10.0)*term2);
+  gsl_matrix_set(J, 3, 1, 0.0);
+  gsl_matrix_set(J, 3, 2, 0.0);
+  gsl_matrix_set(J, 3, 3, -2.0*sqrt(10.0)*term2);
+
+  return GSL_SUCCESS;
+}
+
+static gsl_multifit_function_fdf powell1_func =
+{
+  &powell1_f,
+  &powell1_df,
+  NULL,
+  powell1_N,
+  powell1_P,
+  NULL,
+  0,
+  0
+};
+
+static test_fdf_problem powell1_problem =
+{
+  "powell_singular",
+  powell1_x0,
+  NULL,
+  &powell1_epsrel,
+  powell1_NTRIES,
+  &powell1_checksol,
+  &powell1_func
+};
diff --git a/multifit/test_powell3.c b/multifit/test_powell3.c
index 7c83384..a3f6dd2 100644
--- a/multifit/test_powell3.c
+++ b/multifit/test_powell3.c
@@ -4,7 +4,7 @@
 #define powell3_NTRIES    1
 
 static double powell3_x0[powell3_P] = { 0.0, 1.0 };
-static double powell3_epsrel = 1.0e-12;
+static double powell3_epsrel = 1.0e-8;
 
 static void
 powell3_checksol(const double x[], const double sumsq,
diff --git a/multifit/test_powell3.c.tol b/multifit/test_powell3.c.tol
new file mode 100644
index 0000000..7c83384
--- /dev/null
+++ b/multifit/test_powell3.c.tol
@@ -0,0 +1,77 @@
+#define powell3_N         2
+#define powell3_P         2
+
+#define powell3_NTRIES    1
+
+static double powell3_x0[powell3_P] = { 0.0, 1.0 };
+static double powell3_epsrel = 1.0e-12;
+
+static void
+powell3_checksol(const double x[], const double sumsq,
+                 const double epsrel, const char *sname,
+                 const char *pname)
+{
+  size_t i;
+  const double sumsq_exact = 0.0;
+  const double powell3_x[powell3_P] = { 1.09815932969975976e-05,
+                                        9.10614673986700218 };
+
+  gsl_test_rel(sumsq, sumsq_exact, epsrel, "%s/%s sumsq",
+               sname, pname);
+
+  for (i = 0; i < powell3_P; ++i)
+    {
+      gsl_test_rel(x[i], powell3_x[i], epsrel, "%s/%s i=%zu",
+                   sname, pname, i);
+    }
+}
+
+static int
+powell3_f (const gsl_vector * x, void *params, gsl_vector * f)
+{
+  double x1 = gsl_vector_get(x, 0);
+  double x2 = gsl_vector_get(x, 1);
+
+  gsl_vector_set(f, 0, 1.0e4*x1*x2 - 1.0);
+  gsl_vector_set(f, 1, exp(-x1) + exp(-x2) - 1.0001);
+
+  return GSL_SUCCESS;
+}
+
+static int
+powell3_df (const gsl_vector * x, void *params, gsl_matrix * J)
+{
+  double x1 = gsl_vector_get(x, 0);
+  double x2 = gsl_vector_get(x, 1);
+
+  gsl_matrix_set(J, 0, 0, 1.0e4*x2);
+  gsl_matrix_set(J, 0, 1, 1.0e4*x1);
+
+  gsl_matrix_set(J, 1, 0, -exp(-x1));
+  gsl_matrix_set(J, 1, 1, -exp(-x2));
+
+  return GSL_SUCCESS;
+}
+
+static gsl_multifit_function_fdf powell3_func =
+{
+  &powell3_f,
+  &powell3_df,
+  NULL,
+  powell3_N,
+  powell3_P,
+  NULL,
+  0,
+  0
+};
+
+static test_fdf_problem powell3_problem =
+{
+  "powell_badly_scaled",
+  powell3_x0,
+  NULL,
+  &powell3_epsrel,
+  powell3_NTRIES,
+  &powell3_checksol,
+  &powell3_func
+};