static void compute_diag (const gsl_matrix * J, gsl_vector * diag);

static void update_diag (const gsl_matrix * J, gsl_vector * diag);

static double compute_delta (gsl_vector * diag, gsl_vector * x);

static double scaled_enorm (const gsl_vector * d, const gsl_vector * f);

static double enorm (const gsl_vector * f);

static double

enorm (const gsl_vector * f)

{

  return gsl_blas_dnrm2 (f);

}

static double

scaled_enorm (const gsl_vector * d, const gsl_vector * f)

{

  double e2 = 0;

  size_t i, n = f->size;

  for (i = 0; i < n; i++)

    {

      double fi = gsl_vector_get (f, i);

      double di = gsl_vector_get (d, i);

      double u = di * fi;

      e2 += u * u;

    }

  return sqrt (e2);

}

static double

compute_delta (gsl_vector * diag, gsl_vector * x)

{

  double Dx = scaled_enorm (diag, x);

  double factor = 100;  /* generally recommended value from MINPACK */

  return (Dx > 0) ? factor * Dx : factor;

}

static double

compute_actual_reduction (double fnorm, double fnorm1)

{

  double actred;

  if (0.1 * fnorm1 < fnorm)

    {

      double u = fnorm1 / fnorm;

      actred = 1 - u * u;

    }

  else

    {

      actred = -1;

    }

  return actred;

}

static void

compute_diag (const gsl_matrix * J, gsl_vector * diag)

{

  size_t j, p = J->size2;

  for (j = 0; j < p; j++)

    {

      gsl_vector_const_view v = gsl_matrix_const_column(J, j);

      double norm = gsl_blas_dnrm2(&v.vector);

      if (norm == 0)

        norm = 1.0;

      gsl_vector_set (diag, j, norm);

    }

}

static void

update_diag (const gsl_matrix * J, gsl_vector * diag)

{

  size_t j, p = J->size2;

  for (j = 0; j < p; j++)

    {

      gsl_vector_const_view v = gsl_matrix_const_column(J, j);

      double norm = gsl_blas_dnrm2(&v.vector);

      double *diagj = gsl_vector_ptr(diag, j);

      if (norm == 0)

        norm = 1.0;

      *diagj = GSL_MAX(*diagj, norm);

    }

}

static void

compute_rptdx (const gsl_matrix * r, const gsl_permutation * p,

               const gsl_vector * dx, gsl_vector * rptdx)

{

  size_t i, j, N = dx->size;

  for (i = 0; i < N; i++)

    {

      double sum = 0;

      for (j = i; j < N; j++)

        {

          size_t pj = gsl_permutation_get (p, j);

          sum += gsl_matrix_get (r, i, j) * gsl_vector_get (dx, pj);

        }

      gsl_vector_set (rptdx, i, sum);

    }

}

static void

compute_trial_step (gsl_vector * x, gsl_vector * dx, gsl_vector * x_trial)

{

  size_t i, N = x->size;

  for (i = 0; i < N; i++)

    {

      double pi = gsl_vector_get (dx, i);

      double xi = gsl_vector_get (x, i);

      gsl_vector_set (x_trial, i, xi + pi);

    }

}