const size_t nelson_N = 128;
const size_t nelson_P = 3;

/* double nelson_x0[3] = { 2, 0.0001, -0.01}; */

double nelson_x0[3] = { 2.5 , 0.000000005, -0.01}; 

double nelson_x[3] = {
  2.5906836021E+00,
  5.6177717026E-09,
  -5.7701013174E-02
};

double nelson_sumsq = 3.7976833176E+00;

double nelson_sigma[3] = {
  1.9149996413E-02,
  6.1124096540E-09,
  3.9572366543E-03
};

double nelson_F[128][3] = {
  {  15.00E0,      1E0,      180E0},
  {  17.00E0,      1E0,      180E0},
  {  15.50E0,      1E0,      180E0},
  {  16.50E0,      1E0,      180E0},
  {  15.50E0,      1E0,      225E0},
  {  15.00E0,      1E0,      225E0},
  {  16.00E0,      1E0,      225E0},
  {  14.50E0,      1E0,      225E0},
  {  15.00E0,      1E0,      250E0},
  {  14.50E0,      1E0,      250E0},
  {  12.50E0,      1E0,      250E0},
  {  11.00E0,      1E0,      250E0},
  {  14.00E0,      1E0,      275E0},
  {  13.00E0,      1E0,      275E0},
  {  14.00E0,      1E0,      275E0},
  {  11.50E0,      1E0,      275E0},
  {  14.00E0,      2E0,      180E0},
  {  16.00E0,      2E0,      180E0},
  {  13.00E0,      2E0,      180E0},
  {  13.50E0,      2E0,      180E0},
  {  13.00E0,      2E0,      225E0},
  {  13.50E0,      2E0,      225E0},
  {  12.50E0,      2E0,      225E0},
  {  12.50E0,      2E0,      225E0},
  {  12.50E0,      2E0,      250E0},
  {  12.00E0,      2E0,      250E0},
  {  11.50E0,      2E0,      250E0},
  {  12.00E0,      2E0,      250E0},
  {  13.00E0,      2E0,      275E0},
  {  11.50E0,      2E0,      275E0},
  {  13.00E0,      2E0,      275E0},
  {  12.50E0,      2E0,      275E0},
  {  13.50E0,      4E0,      180E0},
  {  17.50E0,      4E0,      180E0},
  {  17.50E0,      4E0,      180E0},
  {  13.50E0,      4E0,      180E0},
  {  12.50E0,      4E0,      225E0},
  {  12.50E0,      4E0,      225E0},
  {  15.00E0,      4E0,      225E0},
  {  13.00E0,      4E0,      225E0},
  {  12.00E0,      4E0,      250E0},
  {  13.00E0,      4E0,      250E0},
  {  12.00E0,      4E0,      250E0},
  {  13.50E0,      4E0,      250E0},
  {  10.00E0,      4E0,      275E0},
  {  11.50E0,      4E0,      275E0},
  {  11.00E0,      4E0,      275E0},
  {   9.50E0,      4E0,      275E0},
  {  15.00E0,      8E0,      180E0},
  {  15.00E0,      8E0,      180E0},
  {  15.50E0,      8E0,      180E0},
  {  16.00E0,      8E0,      180E0},
  {  13.00E0,      8E0,      225E0},
  {  10.50E0,      8E0,      225E0},
  {  13.50E0,      8E0,      225E0},
  {  14.00E0,      8E0,      225E0},
  {  12.50E0,      8E0,      250E0},
  {  12.00E0,      8E0,      250E0},
  {  11.50E0,      8E0,      250E0},
  {  11.50E0,      8E0,      250E0},
  {   6.50E0,      8E0,      275E0},
  {   5.50E0,      8E0,      275E0},
  {   6.00E0,      8E0,      275E0},
  {   6.00E0,      8E0,      275E0},
  {  18.50E0,     16E0,      180E0},
  {  17.00E0,     16E0,      180E0},
  {  15.30E0,     16E0,      180E0},
  {  16.00E0,     16E0,      180E0},
  {  13.00E0,     16E0,      225E0},
  {  14.00E0,     16E0,      225E0},
  {  12.50E0,     16E0,      225E0},
  {  11.00E0,     16E0,      225E0},
  {  12.00E0,     16E0,      250E0},
  {  12.00E0,     16E0,      250E0},
  {  11.50E0,     16E0,      250E0},
  {  12.00E0,     16E0,      250E0},
  {   6.00E0,     16E0,      275E0},
  {   6.00E0,     16E0,      275E0},
  {   5.00E0,     16E0,      275E0},
  {   5.50E0,     16E0,      275E0},
  {  12.50E0,     32E0,      180E0},
  {  13.00E0,     32E0,      180E0},
  {  16.00E0,     32E0,      180E0},
  {  12.00E0,     32E0,      180E0},
  {  11.00E0,     32E0,      225E0},
  {   9.50E0,     32E0,      225E0},
  {  11.00E0,     32E0,      225E0},
  {  11.00E0,     32E0,      225E0},
  {  11.00E0,     32E0,      250E0},
  {  10.00E0,     32E0,      250E0},
  {  10.50E0,     32E0,      250E0},
  {  10.50E0,     32E0,      250E0},
  {   2.70E0,     32E0,      275E0},
  {   2.70E0,     32E0,      275E0},
  {   2.50E0,     32E0,      275E0},
  {   2.40E0,     32E0,      275E0},
  {  13.00E0,     48E0,      180E0},
  {  13.50E0,     48E0,      180E0},
  {  16.50E0,     48E0,      180E0},
  {  13.60E0,     48E0,      180E0},
  {  11.50E0,     48E0,      225E0},
  {  10.50E0,     48E0,      225E0},
  {  13.50E0,     48E0,      225E0},
  {  12.00E0,     48E0,      225E0},
  {   7.00E0,     48E0,      250E0},
  {   6.90E0,     48E0,      250E0},
  {   8.80E0,     48E0,      250E0},
  {   7.90E0,     48E0,      250E0},
  {   1.20E0,     48E0,      275E0},
  {   1.50E0,     48E0,      275E0},
  {   1.00E0,     48E0,      275E0},
  {   1.50E0,     48E0,      275E0},
  {  13.00E0,     64E0,      180E0},
  {  12.50E0,     64E0,      180E0},
  {  16.50E0,     64E0,      180E0},
  {  16.00E0,     64E0,      180E0},
  {  11.00E0,     64E0,      225E0},
  {  11.50E0,     64E0,      225E0},
  {  10.50E0,     64E0,      225E0},
  {  10.00E0,     64E0,      225E0},
  {   7.27E0,     64E0,      250E0},
  {   7.50E0,     64E0,      250E0},
  {   6.70E0,     64E0,      250E0},
  {   7.60E0,     64E0,      250E0},
  {   1.50E0,     64E0,      275E0},
  {   1.00E0,     64E0,      275E0},
  {   1.20E0,     64E0,      275E0},
  {   1.20E0,     64E0,      275E0}
};


int
nelson_f (const gsl_vector * x, void *params, gsl_vector * f)
{
  double b[3];
  size_t i;

  for (i = 0; i < 3; i++)
    {
      b[i] = gsl_vector_get(x, i);
    }


  for (i = 0; i < 128; i++)
    {
      double x1 = nelson_F[i][1];
      double x2 = nelson_F[i][2];

      double y = b[0] - b[1] * x1 * (b[2]*x2 < -100) ? 0.0 : exp(-b[2] * x2);

      gsl_vector_set (f, i, log(nelson_F[i][0]) - y);
    }

  return GSL_SUCCESS;
}

int
nelson_df (const gsl_vector * x, void *params, gsl_matrix * df)
{
  double b[3];
  size_t i;

  for (i = 0; i < 3; i++)
    {
      b[i] = gsl_vector_get(x, i);
    }

  for (i = 0; i < 128; i++)
    {
      double x1 = nelson_F[i][1];
      double x2 = nelson_F[i][2];

      gsl_matrix_set (df, i, 0, -1.0);
      gsl_matrix_set (df, i, 1, x1 * exp(-b[2] * x2));
      gsl_matrix_set (df, i, 2, -b[1] * x1 * x2 * exp(-b[2] * x2));
    }

  return GSL_SUCCESS;
}

int
nelson_fdf (const gsl_vector * x, void *params,
           gsl_vector * f, gsl_matrix * df)
{
  nelson_f (x, params, f);
  nelson_df (x, params, df);

  return GSL_SUCCESS;
}