Class mroot_hybrids (o2scl)

O2scl : Class List

template<class func_t = mm_funct, class vec_t = boost::numeric::ublas::vector<double>, class mat_t = boost::numeric::ublas::matrix<double>, class jfunc_t = jac_funct>
class mroot_hybrids : public o2scl::mroot<mm_funct, boost::numeric::ublas::vector<double>, jac_funct>, private o2scl::mroot_hybrids_base

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 [Garbow80].

Both the scaled and unscaled options are available by setting 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 the 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 \( \sum_i |f_i|<\) mroot::tol_rel to determine if the solver has succeeded.

See the Multi-dimensional solvers section of the User’s guide for general information about the O2scl solvers. There is an example for the usage of the multidimensional solver classes given in examples/ex_mroot.cpp, see the Multi-dimensional solver example.

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. The algorithm is also modified to check that it is not sending non-finite values to the user-specified function. To return to the default GSL behavior, set shrink_step and extra_finite_check to false.

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

By default convergence failures result in calling the exception handler, but this can be turned off by setting mroot::err_nonconv to false. If mroot::err_nonconv is false, then the functions iterate(), msolve() and msolve_de() will return a non-zero value if convergence fails. Note that the default Jacobian object, def_jac also has a data member jacobian_gsl::err_nonconv which separately handles the case where the one row of the Jacobian is all zero.

Todo

class mroot_hybrids

Future:

• Is all the setting of vectors and matrices to zero really necessary? Do they need to be executed even if memory hasn’t been recently allocated?

• Convert more ubvectors to vec_t.

• Some more of the element-wise vector manipulation could be converted to BLAS routines.

• It’s kind of strange that set() sets jac_given to false and set_de() has to reset it to true. Can this be simplified?

• Many of these minpack functions could be put in their own “minpack_tools” class, or possibly moved to be linear algebra routines instead.

• There are still some numbers in here which the user could have control over, for example, the nslow2 threshold which indicates failure.

Note

The set() and set_de() functions store a pointer to the function object and the user must ensure that the object is still valid for a later call to iterate().

Public Functions

inline mroot_hybrids()

inline virtual ~mroot_hybrids()

inline virtual int set_jacobian(jacobian<func_t, vec_t, mat_t> &j)

Set the automatic Jacobian object.

inline int iterate()

Perform an iteration.

At the end of the iteration, the current value of the solution is stored in x.

inline void allocate(size_t n)

Allocate the memory.

inline virtual const char *type()

Return the type,"mroot_hybrids".

inline virtual int msolve_de(size_t nn, vec_t &xx, func_t &ufunc, jfunc_t &dfunc)

Solve func with derivatives dfunc using x as an initial guess, returning x.

inline virtual int msolve(size_t nn, vec_t &xx, func_t &ufunc)

Solve ufunc using xx as an initial guess, returning xx.

inline int set(size_t nn, vec_t &ax, func_t &ufunc)

Set the function, the parameters, and the initial guess.

inline int set_de(size_t nn, vec_t &ax, func_t &ufunc, jfunc_t &dfunc)

Set the function, the Jacobian, the parameters, and the initial guess.

Public Members

bool shrink_step

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 false.

bool extra_finite_check

If true, double check that the input function values are finite (default true)

bool int_scaling

If true, use the internal scaling method (default true)

jacobian_gsl<func_t, vec_t, mat_t> def_jac

Default automatic Jacobian object.

vec_t f

The value of the function at the present iteration.

vec_t x

The present solution.

Protected Functions

inline virtual int solve_set(size_t nn, vec_t &xx, func_t &ufunc)

Finish the solution after set() or set_de() has been called.

Protected Attributes

int iter

Number of iterations.

size_t ncfail

Compute the number of failures.

size_t ncsuc

Compute the number of successes.

size_t nslow1

The number of times the actual reduction is less than 0.001.

size_t nslow2

The number of times the actual reduction is less than 0.1.

double fnorm

The norm of the current function value.

double delta

The limit of the Nuclidean norm.

mat_t J

Jacobian.

mat_t q

Q matrix from QR decomposition.

mat_t r

R matrix from QR decomposition.

ubvector diag

The diagonal elements.

ubvector qtf

The value of \( Q^T f \).

ubvector newton

The Newton direction.

ubvector gradient

The gradient direction.

ubvector df

The change in the function value.

ubvector qtdf

The value of \( Q^T \cdot \mathrm{df} \).

ubvector rdx

The value of \( R \cdot \mathrm{dx} \).

ubvector w

The value of \( w=(Q^T df-R dx)/|dx| \).

ubvector v

The value of \( v=D^2 dx/|dx| \).

jfunc_t *jac

The user-specified Jacobian.

jacobian<func_t, vec_t, mat_t> *ajac

The automatic Jacobian.

vec_t dx

The value of the derivative.

vec_t x_trial

Trial root.

vec_t f_trial

Trial function value.

size_t dim

The number of equations and unknowns.

bool jac_given

True if the jacobian has been given.

func_t *fnewp

The user-specified function.

bool set_called

True if “set” has been called.

Private Functions

mroot_hybrids(const mroot_hybrids<func_t, vec_t, mat_t, jfunc_t>&)

mroot_hybrids<func_t, vec_t, mat_t, jfunc_t> &operator=(const mroot_hybrids<func_t, vec_t, mat_t, jfunc_t>&)