26 #ifndef O2SCL_MROOT_CERN_H 27 #define O2SCL_MROOT_CERN_H 31 #include <boost/numeric/ublas/vector.hpp> 32 #include <boost/numeric/ublas/matrix.hpp> 34 #include <o2scl/vector.h> 35 #include <o2scl/mroot.h> 96 public mroot<func_t,vec_t,jfunc_t> {
98 #ifndef DOXYGEN_INTERNAL 111 eps=0.1490116119384766e-07;
116 {0,1,2,3,3,3,4,4,4,5,5,5,5,6,6,6,6,7,7,7,7,8,8,8,8,9,9,9,9,9,10,
117 10,10,10,10,11,11,11,11,11,12,12,12,12,12,13,13,13,13,13,14,14,
118 14,14,14,15,15,15,15,15,15,16,16,16,16,16,16,17,17,17,17,17,18,
119 18,18,18,18,18,19,19,19,19,19,19,20,20,20,20,20,20,21,21,21,21,
120 21,21,21,22,22,22,22,22,22,23,23,23,23,23,23,24,24,24,24,24,24,
121 24,25,25,25,25,25,25,26,26,26,26,26,26,26,27,27,27,27,27,27,28,
122 28,28,28,28,28,28,29,29,29,29,29,29,29,30,30,30,30,30,30,30,31,
123 31,31,31,31,31,31,32,32,32,32,32,32,32,33,33,33,33,33,33,33,34,
124 34,34,34,34,34,34,35,35,35,35,35,35,35,36,36,36,36,36,36,36,37,
125 37,37,37,37,37,37,37,38,38,38,38,38,38,38,39,39,39,39,39,39,39,
126 40,40,40,40,40,40,40,40,41,41,41,41,41,41,41,42,42,42,42,42,42,
127 42,42,43,43,43,43,43,43,43,44,44,44,44,44,44,44,44,45,45,45,45,
128 45,45,45,45,46,46,46,46,46,46,46,47,47,47,47,47,47,47,47,48,48,
132 for(
size_t i=0;i<289;i++)
mpt[i]=tmp_mpt[i];
173 return "The function solve() has not been called.";
174 }
else if (
info==1) {
175 return "Test 1 was successful.";
176 }
else if (
info==2) {
177 return "Test 2 was successful.";
178 }
else if (
info==3) {
179 return "Both tests were successful.";
180 }
else if (
info==4) {
181 return "Number of iterations is greater than maxf.";
182 }
else if (
info==5) {
183 return "Approximate Jacobian matrix is singular.";
184 }
else if (
info==6) {
185 return "Iterations are not making good progress.";
186 }
else if (
info==7) {
187 return "Iterations are diverging.";
188 }
else if (
info==8) {
189 return "Either tol_abs is too small or Jacobian is nearly singular.";
190 }
else if (
info==9) {
191 return "Either tol_rel, tol_abs, or the number of vars is not positive.";
204 virtual const char *
type() {
return "mroot_cern"; }
233 virtual int msolve(
size_t nvar, vec_t &x, func_t &func) {
235 int mopt=0, i, j, k, it=0;
239 if (maxf<=0) lmaxf=50*(nvar+3);
244 if (nvar<=0 || this->tol_rel<=0.0 || this->
tol_abs<=0.0) {
246 std::string str=
"Invalid value of tol_rel ("+
dtos(this->
tol_rel)+
248 " in mroot_cern::msolve().";
254 if (nvar<=288) mopt=
mpt[nvar-1];
258 for(i=49;i<=((int)nvar) && done==
false;i++) {
259 double temp=log(((
double)i)+1.0)/((double)(nvar+2*i+1));
268 int iflag=0, numf=0, nfcall=0, nier6=-1, nier7=-1, nier8=0;
269 double fnorm=0.0, difit=0.0, xnorm=0.0;
272 for(i=0;i<((int)nvar);i++) {
273 if (xnorm<fabs(x[i])) {
278 double delta=scale*xnorm;
279 if (
set==
false) delta=
scale;
283 vec_t f(nvar), w0(nvar), w1(nvar), w2(nvar);
285 bool solve_done=
false;
286 while (solve_done==
false) {
299 for(j=0;j<((int)nvar);j++) {
300 for(i=0;i<((int)nvar);i++) {
309 for(k=0;k<((int)nvar);k++) {
316 numf=(int)(((
double)nfcall)/((double)nvar));
317 if (fnorm<fabs(fky)) fnorm=fabs(fky);
321 for(j=k;j<((int)nvar);j++) {
322 for(i=0;i<((int)nvar);i++) {
330 numf=(int)(((
double)nfcall)/((double)nvar));
340 for(i=k;i<((int)nvar);i++)
if (eta<fabs(w0[i])) eta=fabs(w0[i]);
345 for(i=k;i<((int)nvar);i++) {
350 if (w0[k]<0.0) sknorm=-sknorm;
355 for(i=0;i<((int)nvar);i++) {
358 for(j=k;j<((int)nvar);j++) {
359 for(i=0;i<((int)nvar);i++) {
363 for(j=k;j<((int)nvar);j++) {
364 double temp=w0[j]/(sknorm*w0[k]);
365 for(i=0;i<((int)nvar);i++) {
373 double temp2=fky/w0[k];
374 if (h*fabs(temp2)>delta)
375 temp2=(temp2>=0.0) ? fabs(delta/h) : -fabs(delta/h);
376 for(i=0;i<((int)nvar);i++) {
386 for(i=0;i<((int)nvar);i++) {
387 if (xnorm<fabs(w1[i])) xnorm=fabs(w1[i]);
388 if (difit<fabs(x[i]-w1[i])) difit=fabs(x[i]-w1[i]);
394 if(delta<scale*xnorm) delta=scale*xnorm;
398 bool lcv=(fnorm<fnorm1 && difit<difit1 && nsing==0);
403 if (fnorm<fnorm1 || difit<difit1) nier7=0;
404 if (difit>eps*xnorm) nier8=0;
416 if (fnorm<=this->tol_rel &&
info==2)
info=3;
428 if (nsing==((
int)nvar)) {
430 O2SCL_CONV_RET(
"Jacobian matrix singular in mroot_cern::msolve().",
445 std::string s=
"Variable tol_abs too small, J singular, or bad ";
446 s+=
"scaling in mroot_cern::msolve().";
453 O2SCL_ERR(
"Unspecified error in mroot_cern::msolve().",
457 if (!((!lcv) || difit>0.05*xnorm)) {
463 for(
int m=2;m<=mopt && bskip==
false;m++) {
466 for(k=0;k<((int)nvar) && bskip==
false;k++) {
473 numf=(int)(((
double)nfcall)/((double)nvar));
475 if (fnorm<fabs(fky)) fnorm=fabs(fky);
486 double temp3=fky/w0[k];
488 for(i=0;i<((int)nvar);i++) {
501 for(i=0;i<((int)nvar);i++) {
502 if (xnorm<fabs(w1[i])) xnorm=fabs(w1[i]);
503 if (difit<fabs(x[i]-w1[i])) difit=fabs(x[i]-w1[i]);
509 if (fnorm<=this->tol_rel)
info=1;
511 if (fnorm<=this->tol_rel &&
info==2)
info=3;
512 if (numf>=lmaxf &&
info==0) {
529 #ifndef DOXYGEN_INTERNAL iterative process is out of control
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
int verbose
Output control (default 0)
int get_info()
Get the value of INFO from the last call to msolve()
#define O2SCL_CONV_RET(d, n, b)
Set a "convergence" error and return the error value.
std::function< int(size_t, const boost::numeric::ublas::vector< double > &, boost::numeric::ublas::vector< double > &) > mm_funct
Array of multi-dimensional functions typedef.
double scale
The original scale parameter from CERNLIB (default 10.0)
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)
invalid argument supplied by user
apparent singularity detected
exceeded max number of iterations
int mpt[289]
Store the number of function evaluations.
Multidimensional root-finding [abstract base].
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.
virtual const char * type()
Return the type, "mroot_cern".
iteration is not making progress toward solution
double eps
The smallest floating point number (default )
bool err_nonconv
If true, call the error handler if msolve() or msolve_de() does not converge (default true) ...
std::string dtos(double x, int prec=6, bool auto_prec=false)
Convert a double to a string.
std::string get_info_string()
Get the a string corresponding to the integer returned by mroot_cern::get_info(). ...
int maxf
Maximum number of function evaluations.
virtual int msolve(size_t nvar, vec_t &x, func_t &func)
Solve func using x as an initial guess, returning x.
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Multi-dimensional mroot-finding routine (CERNLIB)
double tol_abs
The minimum allowable stepsize (default 1.0e-12)
double tol_rel
The maximum value of the functions for success (default 1.0e-8)
boost::numeric::ublas::matrix< double > w
Desc.
std::string itos(int x)
Convert an integer to a string.
int info
Internal storage for the value of info.