Go to the documentation of this file.
47 #ifndef O2SCL_MROOT_HYBRIDS_H
48 #define O2SCL_MROOT_HYBRIDS_H
52 #include <gsl/gsl_multiroots.h>
53 #include <gsl/gsl_linalg.h>
55 #include <boost/numeric/ublas/vector.hpp>
56 #include <boost/numeric/ublas/vector_proxy.hpp>
57 #include <boost/numeric/ublas/matrix.hpp>
58 #include <boost/numeric/ublas/matrix_proxy.hpp>
60 #include <o2scl/vector.h>
61 #include <o2scl/mroot.h>
62 #include <o2scl/jacobian.h>
65 #include <o2scl/columnify.h>
67 #ifndef DOXYGEN_NO_O2NS
84 double u=fnorm1/fnorm0;
96 double u=fnorm1/fnorm0;
105 template<
class vec2_t,
class mat_t>
107 const ubvector &gradient2, vec2_t &Rg) {
109 for (
size_t i=0;i<N;i++) {
111 for (
size_t j=i;j<N;j++) {
112 sum+=r2(i,j)*gradient2[j];
121 template<
class vec2_t>
123 const ubvector &rdx2,
const vec2_t &dx2,
124 const ubvector &diag2,
double pnorm,
127 for (
size_t i=0;i<n;i++) {
128 double diagi=diag2[i];
129 w2[i]=(qtdf2[i]-rdx2[i])/pnorm;
130 v2[i]=diagi*diagi*dx2[i]/pnorm;
137 template<
class vec2_t,
class mat_t>
139 const vec2_t &dx2,
ubvector &rdx2) {
141 for (
size_t i=0;i<N;i++) {
143 for (
size_t j=i;j<N;j++) {
155 template<
class vec2_t,
class vec3_t>
158 for (
size_t i=0;i<n;i++) {
166 template<
class vec2_t>
170 return (Dx > 0) ? factor*Dx : factor;
177 template<
class vec2_t>
double enorm(
size_t N,
const vec2_t &ff) {
179 for (
size_t i=0;i<N;i++) {
190 for (
size_t i=0;i<n;i++) {
201 template<
class vec2_t>
204 for (
size_t i=0;i<N;i++) {
205 xx_trial[i]=xl[i]+dxl[i];
214 template<
class vec2_t>
217 for (
size_t i=0;i<n;i++) {
218 dfl[i]=ff_trial[i]-fl[i];
225 template<
class mat2_t>
227 for (
size_t j=0;j<n;j++) {
229 for (
size_t i=0;i<n;i++) {
230 const double Jij=J2(i,j);
247 template<
class vec2_t,
class vec3_t,
class vec4_t>
249 const vec3_t &f2, vec4_t &qtf2) {
250 for (
size_t j=0;j<N;j++) {
252 for (
size_t i=0;i<N;i++) {
262 template<
class mat2_t>
264 for (
size_t j=0;j<n;j++) {
265 double cnorm, diagj, sum=0.0;
266 for (
size_t i=0;i<n;i++) {
293 template<
class vec2_t>
295 double beta,
ubvector &gradient2, vec2_t &pp) {
296 for (
size_t i=0;i<N;i++) {
297 pp[i]=alpha*newton2[i]+beta*gradient2[i];
304 template<
class mat_t>
312 o2scl_cblas::o2cblas_Upper,
313 o2scl_cblas::o2cblas_NoTrans,
314 o2scl_cblas::o2cblas_NonUnit,N,N,r2,p);
324 template<
class mat_t>
326 const mat_t &r2,
const ubvector &qtf2,
328 for (
size_t j=0;j<M;j++) {
330 for (
size_t i=0;i<N;i++) {
331 sum+=r2(i,j)*qtf2[i];
342 for (
size_t i=0;i<N;i++) {
343 g[i]=(g[i]/gnorm)/diag2[i];
360 template<
class vec2_t,
class mat_t>
362 const ubvector &diag2,
double delta2,
366 double qnorm, gnorm, sgnorm, bnorm, temp;
373 if (qnorm <= delta2) {
374 for(
size_t i=0;i<n;i++) p[i]=newton2[i];
380 gnorm=
enorm(n,gradient2);
384 double alpha=delta2/qnorm;
400 sgnorm=(gnorm/temp)/temp;
419 double bg=bnorm/gnorm;
420 double bq=bnorm/qnorm;
421 double dq=delta2/qnorm;
423 double sd=sgnorm/delta2;
428 double t2=t1-dq*sd2+sqrt(u*u+(1-dq2)*(1-sd2));
430 double alpha=dq*(1-sd2)/t2;
431 double beta=(1-alpha)*sgnorm;
524 #ifndef DOXYGEN_INTERNAL
595 virtual int solve_set(
size_t nn, vec_t &xx, func_t &ufunc) {
614 for(
size_t i=0;i<nn;i++) {
629 for(
size_t i=0;i<nn;i++) {
635 if (((
int)iter)>=this->ntrial) {
637 "exceeded max. number of iterations.",
659 virtual ~mroot_hybrids() {
707 O2SCL_ERR2(
"Function set() not called or most recent call ",
708 "failed in mroot_hybrids::iterate().",
712 double prered, actred;
713 double pnorm, fnorm1, fnorm1p;
715 double p1=0.1, p5=0.5, p001=0.001, p0001=0.0001;
747 for(
size_t ik=0;ik<
dim;ik++) {
748 if (!std::isfinite(
x_trial[ik])) {
767 std::string str=
"Function returned non-zero value ("+
itos(status)+
768 ") in mroot_hybrids::iterate().";
780 while(status!=0 && bit<20) {
781 for(
size_t ib=0;ib<
dim;ib++) {
791 std::string str=
"No suitable step found, function returned ("+
792 itos(status)+
") in mroot_hybrids::iterate().";
834 if (ratio >= p5 ||
ncsuc > 1) {
837 if (fabs (ratio-1) <= p1) {
844 if (ratio >= p0001) {
845 for(
size_t i=0;i<
dim;i++) {
856 if (actred >= p001) {
872 std::string str=
"Jacobian failed and returned ("+
873 itos(jac_ret)+
") in mroot_hybrids::iterate().";
936 for(
size_t ii=0;ii<n;ii++) {
937 for(
size_t jj=0;jj<n;jj++) {
941 for(
size_t ii=0;ii<n;ii++) {
942 for(
size_t jj=0;jj<n;jj++) {
947 for(
size_t ii=0;ii<
diag.size();ii++)
diag[ii]=0.0;
948 for(
size_t ii=0;ii<
qtf.size();ii++)
qtf[ii]=0.0;
951 for(
size_t ii=0;ii<
df.size();ii++)
df[ii]=0.0;
952 for(
size_t ii=0;ii<
qtdf.size();ii++)
qtdf[ii]=0.0;
953 for(
size_t ii=0;ii<
rdx.size();ii++)
rdx[ii]=0.0;
954 for(
size_t ii=0;ii<
w.size();ii++)
w[ii]=0.0;
955 for(
size_t ii=0;ii<
v.size();ii++)
v[ii]=0.0;
965 for(
size_t i=0;i<n;i++) {
977 virtual const char *
type() {
return "mroot_hybrids"; }
983 virtual int msolve_de(
size_t nn, vec_t &xx, func_t &ufunc,
986 int ret=
set_de(nn,xx,ufunc,dfunc);
996 virtual int msolve(
size_t nn, vec_t &xx, func_t &ufunc) {
998 int ret=
set(nn,xx,ufunc);
1004 for(
size_t i=0;i<nn;i++) xx[i]=
x[i];
1011 int set(
size_t nn, vec_t &ax, func_t &ufunc) {
1022 ajac->set_function(ufunc);
1030 for(
size_t i=0;i<
dim;i++)
x[i]=ax[i];
1032 status=ufunc(
dim,ax,
f);
1038 if (jac_given) status=(*jac)(
dim,ax,
dim,
f,
J);
1039 else status=(*ajac)(
dim,ax,
dim,
f,
J);
1053 for(
size_t i=0;i<
dim;i++)
dx[i]=0.0;
1060 for(
size_t ii=0;ii<
diag.size();ii++)
diag[ii]=0.0;
1078 int set_de(
size_t nn, vec_t &ax, func_t &ufunc, jfunc_t &dfunc) {
1086 int ret=
set(nn,ax,ufunc);
1095 #ifndef DOXYGEN_INTERNAL
1108 #ifndef DOXYGEN_NO_O2NS
1112 #if defined (O2SCL_COND_FLAG) || defined (DOXYGEN)
1113 #if defined (O2SCL_ARMA) || defined (DOXYGEN)
1114 #include <armadillo>
1122 template<
class func_t,
class vec_t,
class mat_t,
class jfunc_t>
1126 virtual void qr_decomp_unpack() {
1127 arma::qr_econ(this->
q,this->
r,this->
J);
1134 #if defined (O2SCL_HAVE_EIGEN) || defined (DOXYGEN)
1135 #include <eigen3/Eigen/Dense>
1143 template<
class func_t,
class vec_t,
class mat_t,
class jfunc_t>
1147 virtual void qr_decomp_unpack() {
1148 Eigen::HouseholderQR<Eigen::MatrixXd> hqr(this->
J);
1149 this->
q=hqr.householderQ();
1150 this->
r=hqr.matrixQR().triangularView<Eigen::Upper>();
1158 #include <o2scl/mroot_special.h>
int dogleg(size_t n, const mat_t &r2, const ubvector &qtf2, const ubvector &diag2, double delta2, ubvector &newton2, ubvector &gradient2, vec2_t &p)
Take a dogleg step.
double scaled_enorm(size_t n, const vec2_t &d, const vec3_t &ff)
Compute the norm of the vector defined by .
bool shrink_step
If true, iterate() will shrink the step-size automatically if the function returns a non-zero value (...
double enorm_sum(size_t n, const ubvector &a, const ubvector &b)
Compute the Euclidean norm of the sum of a and b.
#define O2SCL_CONV2_RET(d, d2, n, b)
Set an error and return the error value, two-string version.
ubvector qtdf
The value of .
virtual int set_jacobian(jacobian< func_t, vec_t, mat_t > &j)
Set the automatic Jacobian object.
int iter
Number of iterations.
Class for automatically computing gradients [abstract base].
vec_t x
The present solution.
int compute_trial_step(size_t N, vec2_t &xl, vec2_t &dxl, vec2_t &xx_trial)
Compute a trial step and put the result in xx_trial.
@ exc_efailed
generic failure
size_t nslow1
The number of times the actual reduction is less than 0.001.
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
void compute_rdx(size_t N, const mat_t &r2, const vec2_t &dx2, ubvector &rdx2)
Compute .
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
void compute_Rg(size_t N, const mat_t &r2, const ubvector &gradient2, vec2_t &Rg)
Compute where g is the gradient gradient .
size_t ncfail
Compute the number of failures.
double compute_actual_reduction(double fnorm0, double fnorm1)
Compute the actual reduction.
bool set_called
True if "set" has been called.
A version of mroot_hybrids which uses Armadillo for the QR decomposition.
jacobian< func_t, vec_t, mat_t > * ajac
The automatic Jacobian.
void minimum_step(const size_t N, double gnorm, const ubvector &diag2, ubvector &g)
Compute the point at which the gradient is minimized.
ubvector diag
The diagonal elements.
int iterate()
Perform an iteration.
void vector_copy(const vec_t &src, vec2_t &dest)
Simple vector copy.
void update_diag(size_t n, const mat2_t &J2, ubvector &diag2)
Update diag.
@ exc_emaxiter
exceeded max number of iterations
void compute_wv(size_t n, const ubvector &qtdf2, const ubvector &rdx2, const vec2_t &dx2, const ubvector &diag2, double pnorm, ubvector &w2, ubvector &v2)
Compute w and v.
double delta
The limit of the Nuclidean norm.
int last_ntrial
The number of iterations for in the most recent minimization.
double tol_rel
The maximum value of the functions for success (default 1.0e-8)
vec_t f_trial
Trial function value.
int newton_direction(const size_t N, const mat_t &r2, const ubvector &qtf2, ubvector &p)
Compute the Gauss-Newton direction.
static const double x2[5]
size_t ncsuc
Compute the number of successes.
A version of mroot_hybrids which uses Eigen for the QR decomposition.
Base functions for mroot_hybrids.
vec_t f
The value of the function at the present iteration.
int ntrial
Maximum number of iterations (default 100)
virtual int msolve(size_t nn, vec_t &xx, func_t &ufunc)
Solve ufunc using xx as an initial guess, returning xx.
ubvector newton
The Newton direction.
void allocate(size_t n)
Allocate the memory.
mat_t q
Q matrix from QR decomposition.
double fnorm
The norm of the current function value.
#define O2SCL_CONV_RET(d, n, b)
Set a "convergence" error and return the error value.
void QR_decomp_unpack(const size_t M, const size_t N, mat_t &A, mat2_t &Q, mat3_t &R)
Compute the unpacked QR decomposition of matrix A.
int compute_df(size_t n, const vec2_t &ff_trial, const vec2_t &fl, ubvector &dfl)
Compute the change in the function value.
double compute_predicted_reduction(double fnorm0, double fnorm1)
Compute the predicted reduction phi1p=|Q^T f + R dx|.
@ gsl_continue
iteration has not converged
Multidimensional root-finding algorithm using Powell's Hybrid method (GSL)
bool extra_finite_check
If true, double check that the input function values are finite (default true)
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.
Base for providing a numerical jacobian [abstract base].
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.
int verbose
Output control (default 0)
vec_t dx
The value of the derivative.
Multidimensional root-finding [abstract base].
@ exc_enoprog
iteration is not making progress toward solution
jfunc_t * jac
The user-specified Jacobian.
ubvector df
The change in the function value.
ubvector qtf
The value of .
Simple automatic Jacobian.
void compute_diag(size_t n, const mat2_t &J2, ubvector &diag2)
Compute diag, the norm of the columns of the Jacobian.
void compute_qtf(size_t N, const vec2_t &q2, const vec3_t &f2, vec4_t &qtf2)
Compute .
@ exc_ebadlen
matrix, vector lengths are not conformant
std::string itos(int x)
Convert an integer to a string.
std::function< int(size_t, const boost::numeric::ublas::vector< double > &, boost::numeric::ublas::vector< double > &) > mm_funct
Array of multi-dimensional functions typedef in src/base/mm_funct.h.
@ exc_enoprogj
jacobian jacobian evaluations are not improving the solution
bool err_nonconv
If true, call the error handler if msolve() or msolve_de() does not converge (default true)
void dtrsv(const enum o2cblas_order order, const enum o2cblas_uplo Uplo, const enum o2cblas_transpose TransA, const enum o2cblas_diag Diag, const size_t M, const size_t N, const mat_t &A, vec_t &X)
Compute .
void scaled_addition(size_t N, double alpha, ubvector &newton2, double beta, ubvector &gradient2, vec2_t &pp)
Form appropriate convex combination of the Gauss-Newton direction and the scaled gradient direction.
bool int_scaling
If true, use the internal scaling method (default true)
ubvector gradient
The gradient direction.
std::function< int(size_t, boost::numeric::ublas::vector< double > &, size_t, boost::numeric::ublas::vector< double > &, boost::numeric::ublas::matrix< double > &) > jac_funct
Jacobian function (not necessarily square) in src/root/jacobian.h.
double compute_delta(size_t n, ubvector &diag2, vec2_t &x2)
Compute delta.
void QR_update(size_t M, size_t N, mat1_t &Q, mat2_t &R, vec1_t &w, vec2_t &v)
Update a QR factorisation for A= Q R, A' = A + u v^T,.
mat_t r
R matrix from QR decomposition.
ubvector rdx
The value of .
int set(size_t nn, vec_t &ax, func_t &ufunc)
Set the function, the parameters, and the initial guess.
int print_iter(size_t n, const vec2_t &x, const vec3_t &y, int iter, double value=0.0, double limit=0.0, std::string comment="")
Print out iteration information.
jacobian_gsl< func_t, vec_t, mat_t > def_jac
Default automatic Jacobian object.
bool jac_given
True if the jacobian has been given.
double enorm(size_t N, const vec2_t &ff)
Compute the Euclidean norm of ff.
virtual const char * type()
Return the type,"mroot_hybrids".
size_t nslow2
The number of times the actual reduction is less than 0.1.
size_t dim
The number of equations and unknowns.
void gradient_direction(const size_t M, const size_t N, const mat_t &r2, const ubvector &qtf2, const ubvector &diag2, ubvector &g)
Compute the gradient direction.
virtual int solve_set(size_t nn, vec_t &xx, func_t &ufunc)
Finish the solution after set() or set_de() has been called.
func_t * fnewp
The user-specified function.
@ exc_ebadfunc
problem with user-supplied function
Documentation generated with Doxygen. Provided under the
GNU Free Documentation License (see License Information).