help-gsl
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

changing and rebuilding the library


From: Aleja Carmento
Subject: changing and rebuilding the library
Date: Mon, 12 Jun 2023 14:03:54 +0200

 <https://stackoverflow.com/posts/76271073/timeline>

If you want to change the gsl library to print the trust radius at every
iteration by adding a fprintf function to the trust.c source file like so :
fprintf(stdout, "trust radius is : %g \n", state->delta);

I added the header file for it like so : #include <stdio.h> and modified
the trust_iterate function like so :

static int trust_iterate(void *vstate, const gsl_vector *swts,
              gsl_multifit_nlinear_fdf *fdf, gsl_vector *x,
              gsl_vector *f, gsl_matrix *J, gsl_vector *g,
              gsl_vector *dx)
{
  int status;
  trust_state_t *state = (trust_state_t *) vstate;
  const gsl_multifit_nlinear_parameters *params = &(state->params);
  const gsl_multifit_nlinear_trs *trs = params->trs;
  gsl_multifit_nlinear_trust_state trust_state;
  gsl_vector *x_trial = state->x_trial;       /* trial x + dx */
  gsl_vector *f_trial = state->f_trial;       /* trial f(x + dx) */
  gsl_vector *diag = state->diag;             /* diag(D) */
  double rho;                                 /* ratio
actual_reduction/predicted_reduction */
  int foundstep = 0;                          /* found step dx */
  int bad_steps = 0;                          /* consecutive rejected steps */

  /* store all state parameters needed by low level methods */
  trust_state.x = x;
  trust_state.f = f;
  trust_state.g = g;
  trust_state.J = J;
  trust_state.diag = state->diag;
  trust_state.sqrt_wts = swts;
  trust_state.mu = &(state->mu);
  trust_state.params = params;
  trust_state.solver_state = state->solver_state;
  trust_state.fdf = fdf;
  trust_state.avratio = &(state->avratio);

  /* initialize trust region subproblem with this Jacobian */
  status = (trs->preloop)(&trust_state, state->trs_state);
  if (status)
    return status;

  /* loop until we find an acceptable step dx */
  while (!foundstep)
    {
      /* calculate new step */
      status = (trs->step)(&trust_state, state->delta, dx, state->trs_state);

      /* occasionally the iterative methods (ie: CG Steihaug) can fail
to find a step,
       * so in this case skip rho calculation and count it as a rejected step */

      if (status == GSL_SUCCESS)
        {
          /* compute x_trial = x + dx */
          trust_trial_step(x, dx, x_trial);

          /* compute f_trial = f(x + dx) */
          status = gsl_multifit_nlinear_eval_f(fdf, x_trial, swts, f_trial);
          if (status)
            return status;

          /* check if step should be accepted or rejected */
          status = trust_eval_step(f, f_trial, g, J, dx, &rho, state);
          if (status == GSL_SUCCESS)
            foundstep = 1;
        }
      else
        {
          /* an iterative TRS method failed to find a step vector */
          rho = -1.0;
        }

      /*
       * update trust region radius: if rho is large,
       * then the quadratic model is a good approximation
       * to the objective function, enlarge trust region.
       * If rho is small (or negative), the model function
       * is a poor approximation so decrease trust region. This
       * can happen even if the step is accepted.
       */
      if (rho > 0.75)
        state->delta *= params->factor_up;
      else if (rho < 0.25)
        state->delta /= params->factor_down;

      if (foundstep)
        {
          /* step was accepted */

          /* compute J <- J(x + dx) */
          status = gsl_multifit_nlinear_eval_df(x_trial, f_trial, swts,
                                                params->h_df, params->fdtype,
                                                fdf, J, state->workn);
          if (status)
            return status;

          /* update x <- x + dx */
          gsl_vector_memcpy(x, x_trial);

          /* update f <- f(x + dx) */
          gsl_vector_memcpy(f, f_trial);

          /* compute new g = J^T f */
          gsl_blas_dgemv(CblasTrans, 1.0, J, f, 0.0, g);

          /* update scaling matrix D */
          (params->scale->update)(J, diag);

          /* step accepted, decrease LM parameter */
          status = nielsen_accept(rho, &(state->mu), &(state->nu));
          if (status)
            return status;

          bad_steps = 0;
        }
      else
        {
          /* step rejected, increase LM parameter */
          status = nielsen_reject(&(state->mu), &(state->nu));
          if (status)
            return status;

          if (++bad_steps > 15)
            {
              /* if more than 15 consecutive rejected steps, report no
progress */
              return GSL_ENOPROG;
            }
        }

      fprintf(stdout, "trust radius is : %g \n", state->delta);

    }

  return GSL_SUCCESS;
} /* trust_iterate() */



I tried to recompile the library using cmake but when i start debugging ,
the function isn't printed out . I have downloaded the gsl library using
vcpkg.


reply via email to

[Prev in Thread] Current Thread [Next in Thread]