/* multimin/directional_minimize.c *  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Fabrice Rossi *  * 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 2 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., 675 Mass Ave, Cambridge, MA 02139, USA. */static voidtake_step (const gsl_vector * x, const gsl_vector * p,           double step, double lambda, gsl_vector * x1, gsl_vector * dx){  gsl_vector_set_zero (dx);  gsl_blas_daxpy (-step * lambda, p, dx);  gsl_vector_memcpy (x1, x);  gsl_blas_daxpy (1.0, dx, x1);}static void intermediate_point (gsl_multimin_function_fdf * fdf,                    const gsl_vector * x, const gsl_vector * p,                    double lambda,                     double pg,                    double stepa, double stepc,                    double fa, double fc,                    gsl_vector * x1, gsl_vector * dx, gsl_vector * gradient,                    double * step, double * f){  double stepb, fb;trial:  {    double u = fabs (pg * lambda * stepc);    stepb = 0.5 * stepc * u / ((fc - fa) + u);  }  take_step (x, p, stepb, lambda, x1, dx);  fb = GSL_MULTIMIN_FN_EVAL_F (fdf, x1);#ifdef DEBUG  printf ("trying stepb = %g  fb = %.18e\n", stepb, fb);#endif  if (fb >= fa  && stepb > 0.0)    {      /* downhill step failed, reduce step-size and try again */      fc = fb;      stepc = stepb;      goto trial;    }#ifdef DEBUG  printf ("ok!\n");#endif  *step = stepb;  *f = fb;  GSL_MULTIMIN_FN_EVAL_DF(fdf, x1, gradient);}static voidminimize (gsl_multimin_function_fdf * fdf,          const gsl_vector * x, const gsl_vector * p,          double lambda,          double stepa, double stepb, double stepc,          double fa, double fb, double fc, double tol,          gsl_vector * x1, gsl_vector * dx1,           gsl_vector * x2, gsl_vector * dx2, gsl_vector * gradient,                    double * step, double * f, double * gnorm){  /* Starting at (x0, f0) move along the direction p to find a minimum     f(x0 - lambda * p), returning the new point x1 = x0-lambda*p,     f1=f(x1) and g1 = grad(f) at x1.  */  double u = stepb;  double v = stepa;  double w = stepc;  double fu = fb;  double fv = fa;  double fw = fc;  double old2 = fabs(w - v);  double old1 = fabs(v - u);  double stepm, fm, pg, gnorm1;  double iter = 0;  gsl_vector_memcpy (x2, x1);  gsl_vector_memcpy (dx2, dx1);  *f = fb;  *step = stepb;  *gnorm = gsl_blas_dnrm2 (gradient);mid_trial:  iter++;  if (iter > 10)    {      return;  /* MAX ITERATIONS */    }  {    double dw = w - u;    double dv = v - u;    double du = 0.0;    double e1 = ((fv - fu) * dw * dw + (fu - fw) * dv * dv);    double e2 = 2.0 * ((fv - fu) * dw + (fu - fw) * dv);    if (e2 != 0.0)      {        du = e1 / e2;      }    if (du > 0 && du < (stepc - stepb) && fabs(du) < 0.5 * old2)      {        stepm = u + du;      }    else if (du < 0 && du > (stepa - stepb) && fabs(du) < 0.5 * old2)      {        stepm = u + du;      }    else if ((stepc - stepb) > (stepb - stepa))      {        stepm = 0.38 * (stepc - stepb) + stepb;      }    else      {        stepm = stepb - 0.38 * (stepb - stepa);      }  }  take_step (x, p, stepm, lambda, x1, dx1);  fm = GSL_MULTIMIN_FN_EVAL_F (fdf, x1);#ifdef DEBUG  printf ("trying stepm = %g  fm = %.18e\n", stepm, fm);#endif  if (fm > fb)    {      if (fm < fv)        {          w = v;          v = stepm;          fw = fv;          fv = fm;        }      else if (fm < fw)        {          w = stepm;          fw = fm;        }      if (stepm < stepb)        {          stepa = stepm;          fa = fm;        }      else        {          stepc = stepm;          fc = fm;        }      goto mid_trial;    }  else if (fm <= fb)    {      old2 = old1;      old1 = fabs(u - stepm);      w = v;      v = u;      u = stepm;      fw = fv;      fv = fu;      fu = fm;      gsl_vector_memcpy (x2, x1);      gsl_vector_memcpy (dx2, dx1);      GSL_MULTIMIN_FN_EVAL_DF (fdf, x1, gradient);      gsl_blas_ddot (p, gradient, &pg);      gnorm1 = gsl_blas_dnrm2 (gradient);#ifdef DEBUG      printf ("p: "); gsl_vector_fprintf(stdout, p, "%g");      printf ("g: "); gsl_vector_fprintf(stdout, gradient, "%g");      printf ("gnorm: %.18e\n", gnorm1);      printf ("pg: %.18e\n", pg);      printf ("orth: %g\n", fabs (pg * lambda/ gnorm1));#endif      *f = fm;      *step = stepm;      *gnorm = gnorm1;      if (fabs (pg * lambda / gnorm1) < tol)        {#ifdef DEBUG          printf("ok!\n");#endif          return; /* SUCCESS */        }      if (stepm < stepb)        {          stepc = stepb;          fc = fb;          stepb = stepm;          fb = fm;        }      else        {          stepa = stepb;          fa = fb;          stepb = stepm;          fb = fm;        }      goto mid_trial;    }}