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,
125 ubvector &w2, ubvector &v2) {
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>
216 const vec2_t &fl, ubvector &dfl) {
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>
306 const ubvector &qtf2, ubvector &p) {
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,
327 const ubvector &diag2, ubvector &g) {
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];
341 const ubvector &diag2, ubvector &g) {
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>
361 int dogleg(
size_t n,
const mat_t &r2,
const ubvector &qtf2,
362 const ubvector &diag2,
double delta2,
363 ubvector &newton2, ubvector &gradient2,
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++) {
617 if (resid<this->tol_rel) status=
success;
622 if (this->verbose>0) {
623 this->print_iter(nn,x,f,iter,resid,this->tol_rel,
629 for(
size_t i=0;i<nn;i++) {
633 this->last_ntrial=iter;
635 if (((
int)iter)>=this->ntrial) {
637 "exceeded max. number of iterations.",
656 extra_finite_check=
true;
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;
723 dogleg(dim,r,qtf,diag,delta,newton,gradient,dx);
736 if (extra_finite_check) {
747 for(
size_t ik=0;ik<dim;ik++) {
748 if (!std::isfinite(x_trial[ik])) {
760 if (shrink_step==
false) {
764 status=(*fnewp)(dim,x_trial,f_trial);
767 std::string str=
"Function returned non-zero value ("+
itos(status)+
768 ") in mroot_hybrids::iterate().";
777 status=(*fnewp)(dim,x_trial,f_trial);
780 while(status!=0 && bit<20) {
781 for(
size_t ib=0;ib<dim;ib++) {
782 x_trial[ib]=(x_trial[ib]+x[ib])/2.0;
784 status=(*fnewp)(dim,x_trial,f_trial);
791 std::string str=
"No suitable step found, function returned ("+
792 itos(status)+
") in mroot_hybrids::iterate().";
804 fnorm1=
enorm(dim,f_trial);
834 if (ratio >= p5 || ncsuc > 1) {
835 delta=GSL_MAX(delta,pnorm/p5);
837 if (fabs (ratio-1) <= p1) {
844 if (ratio >= p0001) {
845 for(
size_t i=0;i<dim;i++) {
856 if (actred >= p001) {
868 if (jac_given) jac_ret=(*jac)(dim,x,dim,f,J);
869 else jac_ret=(*ajac)(dim,x,dim,f,J);
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;
949 for(
size_t ii=0;ii<newton.size();ii++) newton[ii]=0.0;
950 for(
size_t ii=0;ii<gradient.size();ii++) gradient[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);
992 return solve_set(nn,xx,ufunc);
996 virtual int msolve(
size_t nn, vec_t &xx, func_t &ufunc) {
998 int ret=
set(nn,xx,ufunc);
1003 int status=solve_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) {
1030 for(
size_t i=0;i<dim;i++) x[i]=ax[i];
1032 status=ufunc(dim,ax,f);
1035 "mroot_hybrids::set().",
exc_ebadfunc,this->err_nonconv);
1038 if (jac_given) status=(*jac)(dim,ax,dim,f,J);
1039 else status=(*ajac)(dim,ax,dim,f,J);
1043 "mroot_hybrids::set().",
exc_efailed,this->err_nonconv);
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 1102 (
const mroot_hybrids<func_t,vec_t,mat_t,jfunc_t>&);
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 iter
Number of iterations.
double enorm(size_t N, const vec2_t &ff)
Compute the Euclidean norm of ff.
bool shrink_step
If true, iterate() will shrink the step-size automatically if the function returns a non-zero value (...
double compute_delta(size_t n, ubvector &diag2, vec2_t &x2)
Compute delta.
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
vec_t dx
The value of the derivative.
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 .
ubvector qtdf
The value of .
size_t ncsuc
Compute the number of successes.
virtual int set_function(func_t &f)
Set the function to compute the Jacobian of.
#define O2SCL_CONV_RET(d, n, b)
Set a "convergence" error and return the error value.
jacobian_gsl< func_t, vec_t, mat_t > def_jac
Default automatic Jacobian object.
std::function< int(size_t, const boost::numeric::ublas::vector< double > &, boost::numeric::ublas::vector< double > &) > mm_funct
Array of multi-dimensional functions typedef.
size_t ncfail
Compute the number of failures.
void allocate(size_t n)
Allocate the memory.
double fnorm
The norm of the current function value.
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)
void compute_Rg(size_t N, const mat_t &r2, const ubvector &gradient2, vec2_t &Rg)
Compute where g is the gradient gradient .
int compute_df(size_t n, const vec2_t &ff_trial, const vec2_t &fl, ubvector &dfl)
Compute the change in the function value.
ubvector diag
The diagonal elements.
double compute_predicted_reduction(double fnorm0, double fnorm1)
Compute the predicted reduction phi1p=|Q^T f + R dx|.
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.
void update_diag(size_t n, const mat2_t &J2, ubvector &diag2)
Update diag.
#define O2SCL_CONV2_RET(d, d2, n, b)
Set an error and return the error value, two-string version.
Multidimensional root-finding [abstract base].
A version of mroot_hybrids which uses Eigen for the QR decomposition.
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.
jacobian< func_t, vec_t, mat_t > * ajac
The automatic Jacobian.
iteration has not converged
size_t nslow1
The number of times the actual reduction is less than 0.001.
double enorm_sum(size_t n, const ubvector &a, const ubvector &b)
Compute the Euclidean norm of the sum of a and b.
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.
void compute_rdx(size_t N, const mat_t &r2, const vec2_t &dx2, ubvector &rdx2)
Compute .
size_t nslow2
The number of times the actual reduction is less than 0.1.
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.
double delta
The limit of the Nuclidean norm.
vec_t f
The value of the function at the present iteration.
void minimum_step(const size_t N, double gnorm, const ubvector &diag2, ubvector &g)
Compute the point at which the gradient is minimized.
mat_t q
Q matrix from QR decomposition.
virtual int set_jacobian(jacobian< func_t, vec_t, mat_t > &j)
Set the automatic Jacobian object.
bool int_scaling
If true, use the internal scaling method (default true)
Multidimensional root-finding algorithm using Powell's Hybrid method (GSL)
void compute_diag(size_t n, const mat2_t &J2, ubvector &diag2)
Compute diag, the norm of the columns of the Jacobian.
void vector_copy(const vec_t &src, vec2_t &dest)
Simple vector copy.
iteration is not making progress toward solution
matrix, vector lengths are not conformant
vec_t x
The present solution.
A version of mroot_hybrids which uses Armadillo for the QR decomposition.
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
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...
int iterate()
Perform an iteration.
jacobian jacobian evaluations are not improving the solution
Base functions for mroot_hybrids.
jfunc_t * jac
The user-specified Jacobian.
func_t * fnewp
The user-specified function.
ubvector df
The change in the function value.
ubvector gradient
The gradient direction.
ubvector qtf
The value of .
bool extra_finite_check
If true, double check that the input function values are finite (default true)
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,.
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.
ubvector newton
The Newton direction.
virtual const char * type()
Return the type,"mroot_hybrids".
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.
problem with user-supplied function
virtual int msolve(size_t nn, vec_t &xx, func_t &ufunc)
Solve ufunc using xx as an initial guess, returning xx.
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.
mat_t r
R matrix from QR decomposition.
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.
bool set_called
True if "set" has been called.
virtual int solve_set(size_t nn, vec_t &xx, func_t &ufunc)
Finish the solution after set() or set_de() has been called.
bool jac_given
True if the jacobian has been given.
size_t dim
The number of equations and unknowns.
std::string itos(int x)
Convert an integer to a string.
static const double x2[5]
double scaled_enorm(size_t n, const vec2_t &d, const vec3_t &ff)
Compute the norm of the vector defined by .
double compute_actual_reduction(double fnorm0, double fnorm1)
Compute the actual reduction.
void compute_qtf(size_t N, const vec2_t &q2, const vec3_t &f2, vec4_t &qtf2)
Compute .
Simple automatic Jacobian.
Base for providing a numerical jacobian [abstract base].
ubvector rdx
The value of .