gsl_mroot_hybrids.h

/*
  -------------------------------------------------------------------
  
  Copyright (C) 2006, 2007, 2008, Andrew W. Steiner
  
  This file is part of O2scl.
  
  O2scl 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 3 of the License, or
  (at your option) any later version.
  
  O2scl 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 O2scl. If not, see <http://www.gnu.org/licenses/>.

  -------------------------------------------------------------------
*/
#ifndef O2SCL_GSL_MROOT_HYBRIDS_H
#define O2SCL_GSL_MROOT_HYBRIDS_H

#include <gsl/gsl_multiroots.h>
#include <gsl/gsl_linalg.h>

#include <o2scl/ovector_tlate.h>
#include <o2scl/gsl_mroot_hybrids_b.h>
#include <o2scl/mroot.h>
#include <o2scl/jacobian.h>
#include <o2scl/qr.h>

#ifndef DOXYGENP
namespace o2scl {
#endif

  /** \brief State class for \ref gsl_mroot_hybrids

      \todo Improve the documentation in this class
   */
  template<class vec_t=ovector_view, class alloc_vec_t=ovector, 
    class alloc_t=ovector_alloc, class mat_t=omatrix_view, 
    class alloc_mat_t=omatrix, class mat_alloc_t=omatrix_alloc>
    class o2scl_hybrid_state_t {

  public:

  o2scl_hybrid_state_t() {
    dim2=0;
  }

  /// Allocate memory for a solver with \c n variables
  int allocate(size_t n) {

    ma.allocate(J,n,n);

    q=gsl_matrix_calloc(n, n);
    if (q==0) {
      ma.free(J,n);
      O2SCL_ERR2_RET("Failed to allocate space for q",
                     " in o2scl_hybrid_state_t::allocate().",gsl_enomem);
    }
    
    r=gsl_matrix_calloc(n, n);
    if (r==0) {
      ma.free(J,n);
      gsl_matrix_free(q);
      O2SCL_ERR2_RET("Failed to allocate space for r",
                     " in o2scl_hybrid_state_t::allocate().",gsl_enomem);
    }
        
    tau=gsl_vector_calloc(n);
    if (tau==0) {
      ma.free(J,n);
      gsl_matrix_free(q);
      gsl_matrix_free(r);
      O2SCL_ERR2_RET("Failed to allocate space for tau",
                     " in o2scl_hybrid_state_t::allocate().",gsl_enomem);
    }
        
    diag=gsl_vector_calloc(n);
    if (diag==0) {
      ma.free(J,n);
      gsl_matrix_free(q);
      gsl_matrix_free(r);
      gsl_vector_free(tau);
      O2SCL_ERR2_RET("Failed to allocate space for diag",
                     " in o2scl_hybrid_state_t::allocate().",gsl_enomem);
    }
        
    qtf=gsl_vector_calloc(n);
    if (qtf==0) {
      ma.free(J,n);
      gsl_matrix_free(q);
      gsl_matrix_free(r);
      gsl_vector_free(tau);
      gsl_vector_free(diag);
      O2SCL_ERR2_RET("Failed to allocate space for qtf",
                     " in o2scl_hybrid_state_t::allocate().",gsl_enomem);
    }
        
    newton=gsl_vector_calloc(n);
    if (newton==0) {
      ma.free(J,n);
      gsl_matrix_free(q);
      gsl_matrix_free(r);
      gsl_vector_free(tau);
      gsl_vector_free(diag);
      gsl_vector_free(qtf);
      O2SCL_ERR2_RET("Failed to allocate space for newton",
                     " in o2scl_hybrid_state_t::allocate().",gsl_enomem);
    }
        
    gradient=gsl_vector_calloc(n);
    if (gradient==0) {
      ma.free(J,n);
      gsl_matrix_free(q);
      gsl_matrix_free(r);
      gsl_vector_free(tau);
      gsl_vector_free(diag);
      gsl_vector_free(qtf);
      gsl_vector_free(newton);
      O2SCL_ERR2_RET("Failed to allocate space for gradient",
                     " in o2scl_hybrid_state_t::allocate().",gsl_enomem);
    }
        
    df=gsl_vector_calloc(n);
    if (df==0) {
      ma.free(J,n);
      gsl_matrix_free(q);
      gsl_matrix_free(r);
      gsl_vector_free(tau);
      gsl_vector_free(diag);
      gsl_vector_free(qtf);
      gsl_vector_free(newton);
      gsl_vector_free(gradient);
      O2SCL_ERR2_RET("Failed to allocate space for df",
                     " in o2scl_hybrid_state_t::allocate().",gsl_enomem);
    }
        
    qtdf=gsl_vector_calloc(n);
    if (qtdf==0) {
      ma.free(J,n);
      gsl_matrix_free(q);
      gsl_matrix_free(r);
      gsl_vector_free(tau);
      gsl_vector_free(diag);
      gsl_vector_free(qtf);
      gsl_vector_free(newton);
      gsl_vector_free(gradient);
      gsl_vector_free(df);
      O2SCL_ERR2_RET("Failed to allocate space for qtdf",
                     " in o2scl_hybrid_state_t::allocate().",gsl_enomem);
    }

    rdx=gsl_vector_calloc(n);
    if (rdx==0) {
      ma.free(J,n);
      gsl_matrix_free(q);
      gsl_matrix_free(r);
      gsl_vector_free(tau);
      gsl_vector_free(diag);
      gsl_vector_free(qtf);
      gsl_vector_free(newton);
      gsl_vector_free(gradient);
      gsl_vector_free(df);
      gsl_vector_free(qtdf);
      O2SCL_ERR2_RET("Failed to allocate space for rdx",
                     " in o2scl_hybrid_state_t::allocate().",gsl_enomem);
    }

    w=gsl_vector_calloc(n);
    if (w==0) {
      ma.free(J,n);
      gsl_matrix_free(q);
      gsl_matrix_free(r);
      gsl_vector_free(tau);
      gsl_vector_free(diag);
      gsl_vector_free(qtf);
      gsl_vector_free(newton);
      gsl_vector_free(gradient);
      gsl_vector_free(df);
      gsl_vector_free(qtdf);
      gsl_vector_free(rdx);
      O2SCL_ERR2_RET("Failed to allocate space for w",
                     " in o2scl_hybrid_state_t::allocate().",gsl_enomem);
    }

    v=gsl_vector_calloc(n);
    if (v==0) {
      ma.free(J,n);
      gsl_matrix_free(q);
      gsl_matrix_free(r);
      gsl_vector_free(tau);
      gsl_vector_free(diag);
      gsl_vector_free(qtf);
      gsl_vector_free(newton);
      gsl_vector_free(gradient);
      gsl_vector_free(df);
      gsl_vector_free(qtdf);
      gsl_vector_free(rdx);
      gsl_vector_free(w);
      O2SCL_ERR2_RET("Failed to allocate space for v",
                     " in o2scl_hybrid_state_t::allocate().",gsl_enomem);
    }

    dim2=n;
    return GSL_SUCCESS;
  }

  /// Free allocated memory
  int free() {
    if (dim2>0) {
      gsl_vector_free(v);
      gsl_vector_free(w);
      gsl_vector_free(rdx);
      gsl_vector_free(qtdf);
      gsl_vector_free(df);
      gsl_vector_free(gradient);
      gsl_vector_free(newton);
      gsl_vector_free(qtf);
      gsl_vector_free(diag);
      gsl_vector_free(tau);
      gsl_matrix_free(r);
      gsl_matrix_free(q);
      ma.free(J,dim2);
    }
    return 0;
  }

  /// Vector allocator
  alloc_t va;
  /// Matrix allocator
  mat_alloc_t ma;
  /// Number of iterations
  size_t iter;
  /// Desc
  size_t ncfail;
  /// Desc
  size_t ncsuc;
  /// Desc
  size_t nslow1;
  /// Desc
  size_t nslow2;
  /// Desc
  double fnorm;
  /// Desc
  double delta;
  /// Jacobian
  alloc_mat_t J;
  /// Q matrix from QR decomposition
  gsl_matrix *q;
  /// R matrix from QR decomposition
  gsl_matrix *r;
  /// tau vector from QR decomposition
  gsl_vector *tau;
  /// Desc
  gsl_vector *diag;
  /// Desc
  gsl_vector *qtf;
  /// Desc
  gsl_vector *newton;
  /// Desc
  gsl_vector *gradient;
  /// Desc
  gsl_vector *df;
  /// Desc
  gsl_vector *qtdf;
  /// Desc
  gsl_vector *rdx;
  /// Desc
  gsl_vector *w;
  /// Desc
  gsl_vector *v;
  /// Number of variables
  size_t dim2;
  //@}
  };

  /** 
      \brief Multidimensional root-finding algorithm using
      Powell's Hybrid method (GSL)
      
      This is a recasted version of the GSL routines which use a
      modified version of Powell's Hybrid method as implemented in the
      HYBRJ algorithm in MINPACK. Both the scaled and unscaled options
      are available by setting \ref int_scaling (the scaled version is
      the default). If derivatives are not provided, they will be
      computed automatically. This class provides the GSL-like
      interface using allocate(), set() (or set_de() in case where
      derivatives are available), iterate(), and free() and
      higher-level interfaces, msolve() and msolve_de(), which perform
      the solution and the memory allocation automatically. Some
      additional checking is performed in case the user calls the
      functions out of order (i.e. set() without allocate()).
      
      The functions msolve() and msolve_de() use the condition \f$
      \sum_i |f_i| < \f$ \ref mroot::tolf to determine if the solver has
      succeeded.
      
      The original GSL algorithm has been modified to shrink the
      stepsize if a proposed step causes the function to return a
      non-zero value. This allows the routine to automatically try to
      avoid regions where the function is not defined. To return to
      the default GSL behavior, set \ref shrink_step to false.

      The default method for numerically computing the Jacobian
      is from \ref simple_jacobian. This default is identical
      to the GSL approach, except that the default
      value of \ref simple_jacobian::epsmin is non-zero. See
      \ref simple_jacobian for more details.

      There is an example for the usage of the multidimensional solver
      classes given in <tt>examples/ex_mroot.cpp</tt>, see \ref
      ex_mroot_sect .

  */
  template<class param_t, class func_t=mm_funct<param_t>, 
    class vec_t=ovector_view, class alloc_vec_t=ovector, 
    class alloc_t=ovector_alloc, class mat_t=omatrix_view, 
    class alloc_mat_t=omatrix, class mat_alloc_t=omatrix_alloc,
    class jfunc_t=jac_funct<param_t,vec_t,mat_t> > 
    class gsl_mroot_hybrids : 
  public mroot<param_t,func_t,vec_t,jfunc_t>, public hybrids_base {
    
#ifndef DOXYGEN_INTERNAL
    
  protected:

  /// Desc
  void compute_Rg(size_t N, const gsl_matrix * r, 
                  const gsl_vector * gradient, vec_t &Rg) {
    size_t i, j;
      
    for (i = 0; i < N; i++) {
      double sum = 0;
        
      for (j = i; j < N; j++) {
        double gj = gsl_vector_get (gradient, j);
        double rij = gsl_matrix_get (r, i, j);
        sum += rij * gj;
      }
        
      Rg[i]=sum;
    }
  }

  /// Desc
  void compute_wv(size_t n, const gsl_vector * qtdf, 
                  const gsl_vector *rdx, const vec_t &dxx, 
                  const gsl_vector *diag, double pnorm, 
                  gsl_vector * w, gsl_vector * v) {

    size_t i;
      
    for (i = 0; i < n; i++) {
        
      double qtdfi = gsl_vector_get (qtdf, i);
      double rdxi = gsl_vector_get (rdx, i);
      double dxi = dxx[i];
      double diagi = gsl_vector_get (diag, i);
        
      gsl_vector_set (w, i, (qtdfi - rdxi) / pnorm);
      gsl_vector_set (v, i, diagi * diagi * dxi / pnorm);
    }
  }

  /// Desc
  void compute_rdx(size_t N, const gsl_matrix * r, 
                   const vec_t &dxx, gsl_vector * rdx) {
    size_t i, j;
      
    for (i = 0; i < N; i++) {
      double sum = 0;
      for (j = i; j < N; j++) {
        sum += gsl_matrix_get (r, i, j) * dxx[j];
      }
      gsl_vector_set (rdx, i, sum);
    }
  }

  /// Desc
  double scaled_enorm_tvec(size_t n, const gsl_vector * d, const vec_t &ff) {
    double e2 = 0;
    size_t i;
    for (i = 0; i < n ; i++) {
      double fi=ff[i];
      double di= gsl_vector_get(d, i);
      double u = di * fi;
      e2 += u * u ;
    }
    return sqrt(e2);
  }

  /// Desc
  double compute_delta(size_t n, gsl_vector * diag, vec_t &xx) {
    double Dx = scaled_enorm_tvec(n,diag,xx);
    double factor = 100;
      
    return (Dx > 0) ? factor * Dx : factor;
  }

  /// Desc
  double enorm_tvec(const vec_t &ff) {
    double e2 = 0;
    size_t i;
    for (i = 0; i < dim ; i++) {
      double fi= ff[i];
      e2 += fi * fi;
    }
    return sqrt(e2);
  }

  /// Desc
  int compute_trial_step_tvec(size_t N, vec_t &xl, vec_t &dxl, 
                              vec_t &xx_trial) {
    size_t i;
      
    for (i = 0; i < N; i++) {
      xx_trial[i]=xl[i]+dxl[i];
    }
    return 0;
  }

  /// Desc
  int compute_df_tvec(size_t n, const vec_t &ff_trial, 
                      const vec_t &fl, gsl_vector * dfl) {
    size_t i;
      
    for (i = 0; i < n; i++) {
      double dfi=ff_trial[i]-fl[i];
      gsl_vector_set(dfl,i,dfi);
    }

    return 0;
  }
    
  /// Desc
  void compute_diag_tvec(size_t n, const mat_t &J, gsl_vector *diag) {
    size_t i, j;
      
    for (j = 0; j < n; j++) {
      double sum = 0;
      for (i = 0; i < n; i++) {
        double Jij = J[i][j];
        sum += Jij * Jij;
      }
      if (sum == 0) {
        sum = 1.0;
      }
        
      gsl_vector_set(diag,j,sqrt(sum));
    }
  }
    
  /// Desc
  void compute_qtf_tvec(size_t N, const gsl_matrix * q, const vec_t &ff,
                        gsl_vector * qtf) {
    size_t i, j;
      
    for (j = 0; j < N; j++) {
      double sum = 0;
      for (i = 0; i < N; i++) {
        sum += gsl_matrix_get (q, i, j)*ff[i];
      }
        
      gsl_vector_set (qtf, j, sum);
    }
  }

  /// Desc
  void update_diag_tvec(size_t n, const mat_t &J, gsl_vector * diag) {
    size_t i, j;
      
    for (j = 0; j < n; j++) {
      double cnorm, diagj, sum = 0;
      for (i = 0; i < n; i++) {
        double Jij = J[i][j];
        sum += Jij * Jij;
      }
      if (sum == 0) {
        sum = 1.0;
      }
        
      cnorm = sqrt (sum);
      diagj = gsl_vector_get (diag, j);
        
      if (cnorm > diagj) {
        gsl_vector_set (diag, j, cnorm);
      }
    }
  }
    
  /// Desc
  void scaled_addition_tvec(size_t N, double alpha, gsl_vector *newton, 
                            double beta, gsl_vector *gradient, vec_t &pp) {
    size_t i;
      
    for (i = 0; i < N; i++) {
      double ni = gsl_vector_get (newton, i);
      double gi = gsl_vector_get (gradient, i);
      pp[i]=alpha*ni + beta*gi;
    }
  }

  /// Take a dogleg step
  int dogleg(size_t n, const gsl_matrix *r, const gsl_vector *qtf,
             const gsl_vector *diag, double delta,
             gsl_vector *newton, gsl_vector *gradient, 
             vec_t &p) {

    double qnorm, gnorm, sgnorm, bnorm, temp;
      
    newton_direction(r, qtf, newton);
      
    qnorm = scaled_enorm (diag, newton);
      
    if (qnorm <= delta) {
      for(size_t i=0;i<n;i++) p[i]=gsl_vector_get(newton,i);
      return GSL_SUCCESS;
    }
      
    gradient_direction (r, qtf, diag, gradient);
      
    gnorm = enorm (gradient);
      
    if (gnorm == 0) {
      double alpha = delta / qnorm;
      double beta = 0;
      scaled_addition_tvec(n,alpha, newton, beta, gradient, p);
      return GSL_SUCCESS;
    }
      
    minimum_step (gnorm, diag, gradient);
      
    /* Use p as temporary space to compute Rg */
    compute_Rg(n,r, gradient, p);  
      
    temp = enorm_tvec(p);
    sgnorm = (gnorm / temp) / temp;
      
    if (sgnorm > delta)  {
      double alpha = 0;
      double beta = delta;
      scaled_addition_tvec(n,alpha, newton, beta, gradient, p);
      return GSL_SUCCESS;
    }
      
    bnorm = enorm (qtf);
      
    double bg = bnorm / gnorm;
    double bq = bnorm / qnorm;
    double dq = delta / qnorm;
    double dq2 = dq * dq;
    double sd = sgnorm / delta;
    double sd2 = sd * sd;
      
    double t1 = bg * bq * sd;
    double u = t1 - dq;
    double t2 = t1 - dq * sd2 + sqrt (u * u + (1-dq2) * (1 - sd2));
      
    double alpha = dq * (1 - sd2) / t2;
    double beta = (1 - alpha) * sgnorm;
      
    scaled_addition_tvec(n,alpha, newton, beta, gradient, p);
      
    return GSL_SUCCESS;
  }
    
  /// Pointer to the user-specified Jacobian object
  jfunc_t *jac;
    
  /// Pointer to the automatic Jacobian object
  jacobian<param_t,func_t,vec_t,mat_t> *ajac;

  /// Memory allocator for objects of type \c alloc_vec_t
  alloc_t ao;

  /// The value of the derivative
  alloc_vec_t dx;

  /// Trial root
  alloc_vec_t x_trial;
    
  /// Trial function value
  alloc_vec_t f_trial;

  /// The solver state
  o2scl_hybrid_state_t<vec_t,alloc_vec_t,alloc_t,
  mat_t,alloc_mat_t,mat_alloc_t> state;

  /// The function parameters
  param_t *params;

  /// The number of equations and unknowns
  size_t dim;

  /// True if the jacobian has been given
  bool jac_given;

  /// Pointer to the user-specified function
  func_t *fnewp;

  /// True if "set" has been called
  bool set_called;

  /// Finish the solution after set() or set_de() has been called
  int solve_set(size_t nn, vec_t &xx, param_t &pa, func_t &ufunc) {
    int iter=0, status;

    do {
      iter++;
        
      status=iterate();
      if (status) break;

      // ------------------------------------------------------
      // The effective equivalent of the statement:
      // 
      // status=gsl_multiroot_test_residual(f,this->tolf);

      double resid=0.0;
      for(size_t i=0;i<nn;i++) {
        resid+=fabs(f[i]);
      }
      if (resid<this->tolf) status=gsl_success;
      else status=gsl_continue;
        
      // ------------------------------------------------------
        
      if (this->verbose>0) {
        this->print_iter(nn,x,f,iter,resid,this->tolf,
                         "gsl_mroot_hybrids");
      }
        
    } while (status==GSL_CONTINUE && iter < this->ntrial);
      
    for(size_t i=0;i<nn;i++) {
      xx[i]=x[i];
    }

    this->last_ntrial=iter;

    if (iter>=this->ntrial) {
      if (this->last_conv==0) this->last_conv=gsl_emaxiter;
      O2SCL_CONV2_RET("Function gsl_mroot_hybrids::msolve() ",
                      "exceeded max. number of iterations.",
                      gsl_emaxiter,this->err_nonconv);
    }
      
    if (status!=0) {
      O2SCL_ADD_ERR_RET("Function iterate() failed in msolve().",
                        gsl_efailed);
    }

    return gsl_success;
  }

#endif

  public:
      
  gsl_mroot_hybrids() {
    shrink_step=true;
    dim=0;
    ajac=&def_jac;
    def_jac.epsrel=GSL_SQRT_DBL_EPSILON;
    int_scaling=true;
    jac_given=false;
    set_called=false;
  }

  virtual ~gsl_mroot_hybrids() {}

  /** \brief If true, iterate() will shrink the step-size automatically if
      the function returns a non-zero value (default true)

      The original GSL behavior can be obtained by setting 
      this to \c false.
  */
  bool shrink_step;

  /// If true, use the internal scaling method (default true)
  bool int_scaling;

  /// Default automatic Jacobian object
  simple_jacobian<param_t,func_t,vec_t,mat_t,
  alloc_vec_t,alloc_t> def_jac;

  /// Set the automatic Jacobian object
  virtual int set_jacobian(jacobian<param_t,func_t,vec_t,mat_t> &j) {
    ajac=&j;
    return 0;
  }
    
  /** 
      \brief The value of the function at the present iteration
        
      \comment
      We need this to be public so that the user can see if 
      iterate() is working
      \endcomment
        
  */
  alloc_vec_t f;

  /// The present solution
  alloc_vec_t x;

  /** 
      \brief Perform an iteration
        
      At the end of the iteration, the current value of the solution 
      is stored in \ref x.
  */
  int iterate() {
        
    if (!set_called) {
      O2SCL_ERR2_RET("Function set() not called or most recent call ",
                     "failed in gsl_mroot_hybrids::iterate().",
                     gsl_efailed);
    }

    const double fnorm=state.fnorm;

    this->last_conv=0;
 
    gsl_matrix *q=state.q;
    gsl_matrix *r=state.r;
    gsl_vector *tau=state.tau;
    gsl_vector *diag=state.diag;
    gsl_vector *qtf=state.qtf;
    gsl_vector *df=state.df;
    gsl_vector *qtdf=state.qtdf;
    gsl_vector *rdx=state.rdx;
    gsl_vector *w=state.w;
    gsl_vector *v=state.v;
 
    double prered,actred;
    double pnorm,fnorm1,fnorm1p;
    double ratio;
    double p1=0.1,p5=0.5,p001=0.001,p0001=0.0001;
      
    /* Compute qtf=Q^T f */
      
    compute_qtf_tvec(dim,q,f,qtf);
  
    /* Compute dogleg step */
      
    dogleg(dim,r,qtf,diag,state.delta,state.newton,state.gradient,dx);
      
    /* Take a trial step */
      
    compute_trial_step_tvec(dim,x,dx,x_trial);
    pnorm=scaled_enorm_tvec(dim,diag,dx);
    if (state.iter==1) {
      if (pnorm < state.delta) {
        state.delta=pnorm;
      }
    } 

    /* Evaluate function at x + p */
      
    int status;
    
    if (shrink_step==false) {
          
      // Evaluate the function at the new point, exit if it fails
        
      status=(*fnewp)(dim,x_trial,f_trial,*params);
          
      if (status != gsl_success) {
        this->last_conv=gsl_ebadfunc;
        O2SCL_CONV_ADD_ERR_RET
          ("Returned non-zero value in gsl_mroot_hybrids::iterate().",
           o2scl::gsl_ebadfunc,this->err_nonconv);
      }
    
    } else {
          
      // Evaluate the function at the new point, try to recover
      // if it fails
          
      status=(*fnewp)(dim,x_trial,f_trial,*params);
          
      int bit=0;
      while(status!=0 && bit<20) {
        for(size_t ib=0;ib<dim;ib++) {
          x_trial[ib]=(x_trial[ib]+x[ib])/2.0;
        }
        status=(*fnewp)(dim,x_trial,f_trial,*params);
        bit++;
      }
    
      // Exit if we weren't able to find a new good point

      if (status != gsl_success) {
        this->last_conv=gsl_ebadfunc;
        O2SCL_CONV_ADD_ERR_RET
          ("No suitable step found in gsl_mroot_hybrids::iterate().",
           o2scl::gsl_ebadfunc,this->err_nonconv);
      }
    
    }

    /* Set df=f_trial-f */

    compute_df_tvec(dim,f_trial,f,df);
      
    /* Compute the scaled actual reduction */
  
    fnorm1=enorm_tvec(f_trial);
    actred=compute_actual_reduction(fnorm,fnorm1);
      
    /* Compute rdx=R dx */
  
    compute_rdx(dim,r,dx,rdx);
  
    /* Compute the scaled predicted reduction phi1p=|Q^T f + R dx| */
  
    fnorm1p=enorm_sum(qtf,rdx);
    prered=compute_predicted_reduction(fnorm,fnorm1p);

    /* Compute the ratio of the actual to predicted reduction */
  
    if (prered > 0) {
      ratio=actred / prered;
    } else {
      ratio=0;
    }
  
    /* Update the step bound */
  
    if (ratio < p1)     {
      state.ncsuc=0;
      state.ncfail++;
      state.delta *= p5;
    } else {
      state.ncfail=0;
      state.ncsuc++;
      
      if (ratio >= p5 || state.ncsuc > 1) {
        state.delta=GSL_MAX (state.delta,pnorm / p5);
      }
      if (fabs (ratio - 1) <= p1) {
        state.delta=pnorm / p5;
      }
    }
  
    /* Test for successful iteration */

    if (ratio >= p0001) {
      for(size_t i=0;i<dim;i++) {
        x[i]=x_trial[i];
        f[i]=f_trial[i];
      }
      state.fnorm=fnorm1;
      state.iter++;
    }
      
    /* Determine the progress of the iteration */
  
    state.nslow1++;
    if (actred >= p001) {
      state.nslow1=0;
    }
  
    if (actred >= p1) {
      state.nslow2=0;
    }
    
    if (state.ncfail==2) {

      int rx;
          
      if (jac_given) rx=(*jac)(dim,x,f,state.J,*params);
      else rx=(*ajac)(dim,x,f,state.J,*params);

      if (rx!=0) {
        this->last_conv=gsl_efailed;
        O2SCL_CONV_ADD_ERR_RET
          ("Jacobian failed in gsl_mroot_hybrids::iterate().",
           gsl_efailed,this->err_nonconv);
      }
      
      state.nslow2++;
      
      if (state.iter==1) {
        if (int_scaling) {
          compute_diag_tvec(dim,state.J,diag);
        }
        state.delta=compute_delta(dim,diag,x);
      } else {
        if (int_scaling) {
          update_diag_tvec(dim,state.J,diag);
        }
      }
      
      /* Factorize J into QR decomposition */
      ovector *otau=(ovector *)tau;
      omatrix *oq=(omatrix *)q;
      omatrix *oor=(omatrix *)r;
      o2scl_linalg::QR_decomp(dim,dim,state.J,*otau);
      o2scl_linalg::QR_unpack(dim,dim,state.J,*otau,*oq,*oor);
        
      return gsl_success;
    }
  
    /* Compute qtdf=Q^T df,w=(Q^T df - R dx)/|dx|, v=D^2 dx/|dx| */
  
    compute_qtf(q,df,qtdf);
    compute_wv(dim,qtdf,rdx,dx,diag,pnorm,w,v);
  
    /* Rank-1 update of the jacobian Q'R'=Q(R + w v^T) */
  
    gsl_linalg_QR_update(q,r,w,v);
      
    /* No progress as measured by jacobian evaluations */

    if (state.nslow2==5) {
      this->last_conv=gsl_enoprogj;
      O2SCL_CONV_RET
        ("No progress in Jacobian in gsl_mroot_hybrids::iterate().",
         gsl_enoprogj,this->err_nonconv);
    }
  
    /* No progress as measured by function evaluations */

    if (state.nslow1==10) {
      this->last_conv=gsl_enoprog;
      O2SCL_CONV_RET
        ("No progress in function in gsl_mroot_hybrids::iterate().",
         gsl_enoprog,this->err_nonconv);
    }
      
    return gsl_success;
  }
      
  /// Allocate the memory
  int allocate(size_t n) {
    int status;

    if (dim>0) free();
  
    ao.allocate(x,n);
    ao.allocate(dx,n);
    ao.allocate(f,n);
    ao.allocate(x_trial,n);
    ao.allocate(f_trial,n);
    status=state.allocate(n);

    for(size_t i=0;i<n;i++) {
      x[i]=0.0;
      dx[i]=0.0;
      f[i]=0.0;
    }
      
    dim=n;

    return o2scl::gsl_success;
  }
    
  /// Free the allocated memory
  int free() {
    if (dim>0) {
      state.free();
      ao.free(x_trial);
      ao.free(f_trial);
      ao.free(dx);
      ao.free(x);
      ao.free(f);
      dim=0;
    }
    return o2scl::gsl_success;
  }
      
  /// Return the type,\c "gsl_mroot_hybrids".
  virtual const char *type() { return "gsl_mroot_hybrids"; }
    
  /** \brief Solve \c func with derivatives \c dfunc using \c x as 
      an initial guess, returning \c x.

  */
  virtual int msolve_de(size_t nn, vec_t &xx, param_t &pa, func_t &ufunc,
                        jfunc_t &dfunc) {
    int status;
  
    status=set_de(nn,xx,ufunc,dfunc,pa);
    if (status!=0) {
      O2SCL_ADD_ERR_RET("Function set() failed in msolve().",gsl_efailed);
    }
      
    return solve_set(nn,xx,pa,ufunc);
  }

  /// Solve \c ufunc using \c xx as an initial guess, returning \c xx.
  virtual int msolve(size_t nn, vec_t &xx, param_t &pa, func_t &ufunc) {
      
    int status;

    status=set(nn,xx,ufunc,pa);
    if (status!=0) {
      O2SCL_ADD_ERR_RET("Function set() failed in msolve().",gsl_efailed);
    }
      
    status=solve_set(nn,xx,pa,ufunc);
    for(size_t i=0;i<nn;i++) xx[i]=x[i];

    return status;
  }
      
  /** \brief Set the function, the parameters, and the initial guess 
   */
  int set(size_t nn, vec_t &ax, func_t &ufunc, param_t &pa) {

    int status;
  
    if (nn!=dim) { 
      status=allocate(nn);
      if (status!=0) {
        O2SCL_ADD_ERR_RET("Function allocate() failed in msolve().",
                          gsl_efailed);
      }
    }
      
    fnewp=&ufunc;

    // Specify function for automatic Jacobian
    ajac->set_function(ufunc);
      
    if (dim==0) {
      O2SCL_ERR_RET("Memory has not been allocated. Use allocate().",
                    o2scl::gsl_ebadlen);
    }
      
    params=&pa;
      
    // Copy the user-specified solution
    for(size_t i=0;i<dim;i++) x[i]=ax[i];
        
    gsl_matrix *q=state.q;
    gsl_matrix *r=state.r;
    gsl_vector *tau=state.tau;
    gsl_vector *diag=state.diag;
      
    status=ufunc(dim,ax,f,*params);
    if (status!=0) {
      O2SCL_ADD_ERR_RET("Function returned non-zero in set().",
                        gsl_ebadfunc);
    }
      
    if (jac_given) status=(*jac)(dim,ax,f,state.J,*params);
    else status=(*ajac)(dim,ax,f,state.J,*params);

    if (status!=0) {
      O2SCL_ADD_ERR_RET("Jacobian failed in set().",gsl_efailed);
    }

    state.iter=1;
    state.fnorm=enorm_tvec(f);
    state.ncfail=0;
    state.ncsuc=0;
    state.nslow1=0;
    state.nslow2=0;
      
    for(size_t i=0;i<dim;i++) dx[i]=0.0;
      
    /* Store column norms in diag */
  
    if (int_scaling) {
      compute_diag_tvec(dim,state.J,diag);
    } else {
      gsl_vector_set_all(diag,1.0);
    }
        
    /* Set delta to factor |D x| or to factor if |D x| is zero */
        
    state.delta=compute_delta(dim,diag,x);
  
    /* Factorize J into QR decomposition */

    ovector *otau=(ovector *)tau;
    omatrix *oq=(omatrix *)q;
    omatrix *oor=(omatrix *)r;
      
    status=o2scl_linalg::QR_decomp(dim,dim,state.J,*otau);
    if (status!=0) {
      O2SCL_ADD_ERR_RET("QR decomposition failed in set().",gsl_efailed);
    }
    
    status=o2scl_linalg::QR_unpack(dim,dim,state.J,*otau,*oq,*oor);
    if (status!=0) {
      O2SCL_ADD_ERR_RET("QR unpacking failed in set().",gsl_efailed);
    }
      
    set_called=true;
    jac_given=false;

    return gsl_success;
  }

  /** \brief Set the function, the Jacobian, the parameters,
      and the initial guess 
  */
  int set_de(size_t nn, vec_t &ax, func_t &ufunc, jfunc_t &dfunc, 
             param_t &pa) {

    fnewp=&ufunc;
    jac=&dfunc;

    // Make sure set() uses the right Jacobian
    jac_given=true;
      
    int ret=set(nn,ax,ufunc,pa);

    // Reset jac_given since set() will set it back to false
    jac_given=true;
    set_called=true;
    
    return ret;
  }
  
  };
  
#ifndef DOXYGENP
}
#endif

#endif

---

Documentation generated with Doxygen and provided under the GNU Free Documentation License. See License Information for details.

Project hosting provided by , O2scl Sourceforge Project Page