text: #define wnlin_N         40
  #define wnlin_P         3
  
  /* initial guess should be chosen so that Jacobian has full rank,
   * or some solvers will fail */
  static double wnlin_x0[wnlin_P] = { 1.0, 0.9, 0.0 };
  static double wnlin_epsrel = 1.0e-7;
  
  static int wnlin_internal_weight = 1;
  
  /* data */
  static double wnlin_Y[wnlin_N] = {
    6.08035e+00, 5.47552e+00, 5.94654e+00, 5.04920e+00, 4.78568e+00,
    3.51748e+00, 2.84671e+00, 3.24634e+00, 3.23395e+00, 3.30385e+00,
    2.83439e+00, 2.31891e+00, 2.33858e+00, 2.40559e+00, 2.41856e+00,
    1.99966e+00, 1.88127e+00, 1.91477e+00, 1.70415e+00, 1.60316e+00,
    1.77937e+00, 1.55302e+00, 1.50903e+00, 1.36364e+00, 1.36873e+00,
    1.41954e+00, 1.37778e+00, 1.23573e+00, 1.28524e+00, 1.46327e+00,
    1.22315e+00, 1.19330e+00, 1.18717e+00, 8.83172e-01, 1.23424e+00,
    1.14683e+00, 1.11091e+00, 1.20396e+00, 1.28722e+00, 1.05801e+00
  };
  
  /* weights */
  static double wnlin_W[wnlin_N] = {
    2.77778e+00, 3.27690e+00, 3.85426e+00, 4.51906e+00, 5.28083e+00,
    6.14919e+00, 7.13370e+00, 8.24349e+00, 9.48703e+00, 1.08717e+01,
    1.24036e+01, 1.40869e+01, 1.59238e+01, 1.79142e+01, 2.00553e+01,
    2.23415e+01, 2.47646e+01, 2.73137e+01, 2.99753e+01, 3.27337e+01,
    3.55714e+01, 3.84696e+01, 4.14085e+01, 4.43678e+01, 4.73278e+01,
    5.02690e+01, 5.31731e+01, 5.60234e+01, 5.88046e+01, 6.15036e+01,
    6.41092e+01, 6.66121e+01, 6.90054e+01, 7.12839e+01, 7.34442e+01,
    7.54848e+01, 7.74053e+01, 7.92069e+01, 8.08918e+01, 8.24632e+01
  };
  
  static void
  wnlin_checksol(const double x[], const double sumsq,
                 const double epsrel, const char *sname,
                 const char *pname)
  {
    size_t i;
    const double sumsq_exact = 29.7481259665713758;
    const double wnlin_x[wnlin_P] = { 5.17378551196259195,
                                      0.111041758006851149,
                                      1.05282724070446099 };
  
    gsl_test_rel(sumsq, sumsq_exact, epsrel, "%s/%s sumsq",
                 sname, pname);
  
    for (i = 0; i < wnlin_P; ++i)
      {
        gsl_test_rel(x[i], wnlin_x[i], epsrel, "%s/%s i=%zu",
                     sname, pname, i);
      }
  }
  
  static int
  wnlin_f (const gsl_vector *x, void *params, gsl_vector *f)
  {
    int *iptr = (int *) params;
    int doweight = iptr ? *iptr : 0;
    double A = gsl_vector_get (x, 0);
    double lambda = gsl_vector_get (x, 1);
    double b = gsl_vector_get (x, 2);
    size_t i;
  
    /* model Yi = A * exp(-lambda * i) + b */
    for (i = 0; i < wnlin_N; i++)
      {
        double ti = i;
        double yi = wnlin_Y[i];
        double swi = sqrt(wnlin_W[i]);
        double Mi = A * exp (-lambda * ti) + b;
  
        if (doweight)
          gsl_vector_set (f, i, swi * (Mi - yi));
        else
          gsl_vector_set (f, i, Mi - yi);
      }
  
    return GSL_SUCCESS;
  }
  
  static int
  wnlin_df (const gsl_vector *x, void *params, gsl_matrix *df)
  {
    int *iptr = (int *) params;
    int doweight = iptr ? *iptr : 0;
    double A = gsl_vector_get (x, 0);
    double lambda = gsl_vector_get (x, 1);
    size_t i;
  
    for (i = 0; i < wnlin_N; i++)
      {
        gsl_vector_view v = gsl_matrix_row(df, i);
        double ti = i;
        double swi = sqrt(wnlin_W[i]);
        double e = exp(-lambda * ti);
  
        gsl_vector_set(&v.vector, 0, e);
        gsl_vector_set(&v.vector, 1, -ti * A * e);
        gsl_vector_set(&v.vector, 2, 1.0);
  
        if (doweight)
          gsl_vector_scale(&v.vector, swi);
      }
  
    return GSL_SUCCESS;
  }
  
  static int
  wnlin_fvv (const gsl_vector * x, const gsl_vector * v,
             void *params, gsl_vector * fvv)
  {
    int *iptr = (int *) params;
    int doweight = iptr ? *iptr : 0;
    double A = gsl_vector_get (x, 0);
    double lambda = gsl_vector_get (x, 1);
    double v1 = gsl_vector_get(v, 0);
    double v2 = gsl_vector_get(v, 1);
    size_t i;
  
    for (i = 0; i < wnlin_N; i++)
      {
        double ti = i;
        double swi = sqrt(wnlin_W[i]);
        double fvvi;
  
        fvvi = exp(-ti*lambda)*ti*v2 * (-2*v1 + ti*v2*A);
  
        if (doweight)
          gsl_vector_set(fvv, i, swi * fvvi);
        else
          gsl_vector_set(fvv, i, fvvi);
      }
  
    return GSL_SUCCESS;
  }
  
  static gsl_multifit_nlinear_fdf wnlin_func1 =
  {
    wnlin_f,
    wnlin_df,
    wnlin_fvv,
    wnlin_N,
    wnlin_P,
    (void *) &wnlin_internal_weight,
    0,
    0,
    0
  };
  
  static gsl_multifit_nlinear_fdf wnlin_func2 =
  {
    wnlin_f,
    wnlin_df,
    wnlin_fvv,
    wnlin_N,
    wnlin_P,
    NULL,
    0,
    0,
    0
  };
  
  static test_fdf_problem wnlin_problem1 =
  {
    "wnlin_internal_weights",
    wnlin_x0,
    NULL,
    NULL,
    &wnlin_epsrel,
    &wnlin_checksol,
    &wnlin_func1
  };
  
  static test_fdf_problem wnlin_problem2 =
  {
    "wnlin_external_weights",
    wnlin_x0,
    wnlin_W,
    NULL,
    &wnlin_epsrel,
    &wnlin_checksol,
    &wnlin_func2
  };
  