mroot_cern.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2017, Andrew W. Steiner
5 
6  This file is part of O2scl.
7 
8  O2scl is free software; you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version.
12 
13  O2scl is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with O2scl. If not, see <http://www.gnu.org/licenses/>.
20 
21  -------------------------------------------------------------------
22 */
23 /** \file mroot_cern.h
24  \brief File defining \ref o2scl::mroot_cern
25 */
26 #ifndef O2SCL_MROOT_CERN_H
27 #define O2SCL_MROOT_CERN_H
28 
29 #include <string>
30 
31 #include <boost/numeric/ublas/vector.hpp>
32 #include <boost/numeric/ublas/matrix.hpp>
33 
34 #include <o2scl/vector.h>
35 #include <o2scl/mroot.h>
36 
37 namespace o2scl {
38 
39  /** \brief Multi-dimensional mroot-finding routine (CERNLIB)
40 
41  If \f$ x_i \f$ denotes the current iteration, and \f$
42  x^{\prime}_i \f$ denotes the previous iteration, then the
43  calculation is terminated if either of the following tests is
44  successful
45  \f[
46  1:\quad \mathrm{max} | f_i(x) | \leq \mathrm{tol\_rel}
47  \f]
48  \f[
49  2:\quad \mathrm{max} |x_i-x^{\prime}_i| \leq
50  \mathrm{tol\_abs} \times \mathrm{max} | x_i |
51  \f]
52 
53  This routine treats the functions specified as a \ref mm_funct11
54  object slightly differently than \ref o2scl::mroot_hybrids. First
55  the equations should be numbered (as much as is possible) in
56  order of increasing nonlinearity. Also, instead of calculating
57  all of the equations on each function call, only the equation
58  specified by the \c size_t parameter needs to be calculated. If
59  the equations are specified as
60  \f{eqnarray*}
61  &0=f_0(x_0,x_1,...,x_{n-1})& \\
62  &0=f_1(x_0,x_1,...,x_{n-1})& \\
63  &...& \\
64  &0=f_{n-1}(x_0,x_1,...,x_{n-1})& \\
65  \f}
66  then when the \c size_t argument is given as \c i, then
67  only the function \f$ f_i \f$ needs to be calculated.
68 
69  \warning This code has not been checked to ensure that it cannot
70  fail to solve the equations without calling the error handler
71  and returning a non-zero value. Until then, the solution may
72  need to be checked explicitly by the caller.
73 
74  See the \ref multisolve_subsect section of the User's guide for
75  general information about \o2 solvers.
76  There is an example for the usage of the multidimensional solver
77  classes given in <tt>examples/ex_mroot.cpp</tt>, see \ref
78  ex_mroot_sect .
79 
80  \future Modify this so it handles functions which return
81  non-zero values.
82  \future Move some of the memory allocation out of msolve()
83  \future Give the user access to the number of function
84  calls
85  \future Rename nier6, nier7, and nier8 to something sensible.
86  \future It may be that the \o2 native Householder transformations
87  should be used here instead of the inline version given here.
88 
89  Based on the CERNLIB routines RSNLEQ and DSNLEQ, which was
90  based on \ref More79 and \ref More80 and is documented at
91  http://wwwasdoc.web.cern.ch/wwwasdoc/shortwrupsdir/c201/top.html
92  */
93  template<class func_t=mm_funct11,
95  class jfunc_t=jac_funct11> class mroot_cern :
96  public mroot<func_t,vec_t,jfunc_t> {
97 
98 #ifndef DOXYGEN_INTERNAL
99 
100  protected:
101 
102  /// Desc
104 
105 #endif
106 
107  public:
108 
109  mroot_cern() {
110  info=0;
111  eps=0.1490116119384766e-07;
112  scale=10.0;
113  maxf=0;
114 
115  int tmp_mpt[289]=
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,
129  48,48,48,48,48,48};
130  // The for loop is just a convenient way of using
131  // aggregate initialization
132  for(size_t i=0;i<289;i++) mpt[i]=tmp_mpt[i];
133  }
134 
135  /** \brief Get the value of \c INFO from the last call to msolve()
136 
137  The value of info is assigned according to the following list.
138  The values 1-8 are the standard behavior from CERNLIB.
139  0 - The function solve() has not been called.
140  1 - Test 1 was successful. \n
141  2 - Test 2 was successful. \n
142  3 - Both tests were successful. \n
143  4 - Number of iterations is greater than mroot_cern_root::maxf. \n
144  5 - Approximate (finite difference) Jacobian matrix is
145  singular. \n
146  6 - Iterations are not making good progress. \n
147  7 - Iterations are diverging. \n
148  8 - Iterations are converging, but either mroot_cern_root::tol_abs
149  is too small or the Jacobian is nearly singular
150  or the variables are badly scaled. \n
151  9 - Either root::tol_rel or root::tol_abs is not greater than zero
152  or the specified number of variables is \f$ \leq 0\f$.
153 
154  The return values returned by msolve() corresponding
155  to the values of \c INFO above are
156  1 - \ref success
157  2 - \ref success
158  3 - \ref success
159  4 - \ref exc_emaxiter
160  5 - \ref exc_esing
161  6 - \ref exc_enoprog
162  7 - \ref exc_erunaway
163  8 - \ref exc_efailed
164  9 - \ref exc_einval
165  */
166  int get_info() { return info; }
167 
168  /** \brief Get the a string corresponding to the integer returned
169  by \ref mroot_cern::get_info().
170  */
171  std::string get_info_string() {
172  if (info==0) {
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.";
192  }
193  }
194 
195  /** \brief Maximum number of function evaluations
196 
197  If \f$ \mathrm{maxf}\leq 0 \f$ , then \f$ 50(\mathrm{nv}+3) \f$
198  (which is the CERNLIB default) is used. The default value of
199  \c maxf is zero which then implies the default from CERNLIB.
200  */
201  int maxf;
202 
203  /// Return the type, \c "mroot_cern".
204  virtual const char *type() { return "mroot_cern"; }
205 
206  /** \brief The original scale parameter from CERNLIB (default 10.0)
207  */
208  double scale;
209 
210  /** \brief The smallest floating point number
211  (default \f$ \sim 1.49012 \times 10^{-8} \f$ )
212 
213  The original prescription from CERNLIB for \c eps is
214  given below:
215  \verbatim
216  #if !defined(CERNLIB_DOUBLE)
217  PARAMETER (EPS = 0.84293 69702 17878 97282 52636 392E-07)
218  #endif
219  #if defined(CERNLIB_IBM)
220  PARAMETER (EPS = 0.14901 16119 38476 562D-07)
221  #endif
222  #if defined(CERNLIB_VAX)
223  PARAMETER (EPS = 0.37252 90298 46191 40625D-08)
224  #endif
225  #if (defined(CERNLIB_UNIX))&&(defined(CERNLIB_DOUBLE))
226  PARAMETER (EPS = 0.14901 16119 38476 600D-07)
227  #endif
228  \endverbatim
229  */
230  double eps;
231 
232  /// Solve \c func using \c x as an initial guess, returning \c x.
233  virtual int msolve(size_t nvar, vec_t &x, func_t &func) {
234 
235  int mopt=0, i, j, k, it=0;
236  double fky;
237 
238  int lmaxf;
239  if (maxf<=0) lmaxf=50*(nvar+3);
240  else lmaxf=maxf;
241 
242  info=0;
243 
244  if (nvar<=0 || this->tol_rel<=0.0 || this->tol_abs<=0.0) {
245  info=9;
246  std::string str="Invalid value of tol_rel ("+dtos(this->tol_rel)+
247  "), tol_abs ("+dtos(this->tol_abs)+"), or nvar ("+itos(nvar)+
248  " in mroot_cern::msolve().";
249  O2SCL_ERR(str.c_str(),exc_einval);
250  }
251 
252  // Find optimal value of mopt for iterative refinement
253 
254  if (nvar<=288) mopt=mpt[nvar-1];
255  else {
256  bool done=false;
257  double h=0.0;
258  for(i=49;i<=((int)nvar) && done==false;i++) {
259  double temp=log(((double)i)+1.0)/((double)(nvar+2*i+1));
260  if (temp<h) {
261  mopt=i-1;
262  done=true;
263  }
264  if (!done) h=temp;
265  }
266  }
267 
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;
270  bool set=false;
271 
272  for(i=0;i<((int)nvar);i++) {
273  if (xnorm<fabs(x[i])) {
274  xnorm=fabs(x[i]);
275  set=true;
276  }
277  }
278  double delta=scale*xnorm;
279  if (set==false) delta=scale;
280 
281  w.resize(nvar,nvar);
282 
283  vec_t f(nvar), w0(nvar), w1(nvar), w2(nvar);
284 
285  bool solve_done=false;
286  while (solve_done==false) {
287  bool bskip=false;
288 
289  int nsing=nvar;
290  double fnorm1=fnorm;
291  double difit1=difit;
292  fnorm=0.0;
293 
294  // Compute step H for the divided difference which approximates
295  // the K-th row of the Jacobian matrix
296 
297  double h=eps*xnorm;
298  if (h==0.0) h=eps;
299  for(j=0;j<((int)nvar);j++) {
300  for(i=0;i<((int)nvar);i++) {
301  w(j,i)=0.0;
302  }
303  w(j,j)=h;
304  w1[j]=x[j];
305  }
306 
307  // Enter a subiteration
308 
309  for(k=0;k<((int)nvar);k++) {
310  iflag=k;
311 
312  func(iflag,w1,f);
313 
314  fky=f[k];
315  nfcall++;
316  numf=(int)(((double)nfcall)/((double)nvar));
317  if (fnorm<fabs(fky)) fnorm=fabs(fky);
318 
319  // Compute the K-th row of the Jacobian matrix
320 
321  for(j=k;j<((int)nvar);j++) {
322  for(i=0;i<((int)nvar);i++) {
323  w2[i]=w1[i]+w(j,i);
324  }
325 
326  func(iflag,w2,f);
327 
328  double fkz=f[k];
329  nfcall++;
330  numf=(int)(((double)nfcall)/((double)nvar));
331  w0[j]=fkz-fky;
332  }
333 
334  f[k]=fky;
335 
336  // Compute the Householder transformation to reduce the K-th row
337  // of the Jacobian matrix to a multiple of the K-th unit vector
338 
339  double eta=0.0;
340  for(i=k;i<((int)nvar);i++) if (eta<fabs(w0[i])) eta=fabs(w0[i]);
341 
342  if (eta!=0.0) {
343  nsing--;
344  double sknorm=0.0;
345  for(i=k;i<((int)nvar);i++) {
346  w0[i]/=eta;
347  sknorm+=w0[i]*w0[i];
348  }
349  sknorm=sqrt(sknorm);
350  if (w0[k]<0.0) sknorm=-sknorm;
351  w0[k]+=sknorm;
352 
353  // Apply the transformation
354 
355  for(i=0;i<((int)nvar);i++) {
356  w2[i]=0.0;
357  }
358  for(j=k;j<((int)nvar);j++) {
359  for(i=0;i<((int)nvar);i++) {
360  w2[i]+=w0[j]*w(j,i);
361  }
362  }
363  for(j=k;j<((int)nvar);j++) {
364  double temp=w0[j]/(sknorm*w0[k]);
365  for(i=0;i<((int)nvar);i++) {
366  w(j,i)-=temp*w2[i];
367  }
368  }
369 
370  // Compute the subiterate
371 
372  w0[k]=sknorm*eta;
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++) {
377  w1[i]+=temp2*w(k,i);
378  }
379  }
380  }
381 
382  // Compute the norms of the iterate and correction vector
383 
384  xnorm=0.0;
385  difit=0.0;
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]);
389  x[i]=w1[i];
390  }
391 
392  // Update the bound on the correction vector
393 
394  if(delta<scale*xnorm) delta=scale*xnorm;
395 
396  // Determine the progress of the iteration
397 
398  bool lcv=(fnorm<fnorm1 && difit<difit1 && nsing==0);
399  nier6++;
400  nier7++;
401  nier8++;
402  if (lcv) nier6=0;
403  if (fnorm<fnorm1 || difit<difit1) nier7=0;
404  if (difit>eps*xnorm) nier8=0;
405 
406  // Print iteration information
407 
408  if (this->verbose>0) {
409  this->print_iter(nvar,x,f,++it,fnorm,this->tol_rel,"mroot_cern");
410  }
411 
412  // Tests for convergence
413 
414  if (fnorm<=this->tol_rel) info=1;
415  if (difit<=this->tol_abs*xnorm && lcv) info=2;
416  if (fnorm<=this->tol_rel && info==2) info=3;
417  if (info!=0) {
418  return 0;
419  }
420 
421  // Tests for termination
422 
423  if (numf>=lmaxf) {
424  info=4;
425  O2SCL_CONV_RET("Too many iterations in mroot_cern::msolve().",
426  exc_emaxiter,this->err_nonconv);
427  }
428  if (nsing==((int)nvar)) {
429  info=5;
430  O2SCL_CONV_RET("Jacobian matrix singular in mroot_cern::msolve().",
431  exc_esing,this->err_nonconv);
432  }
433  if (nier6==5) {
434  info=6;
435  O2SCL_CONV_RET("No progress in mroot_cern::msolve().",
436  exc_enoprog,this->err_nonconv);
437  }
438  if (nier7==3) {
439  info=7;
440  O2SCL_CONV_RET("Iterations diverging in mroot_cern::msolve().",
441  exc_erunaway,this->err_nonconv);
442  }
443  if (nier8==4) {
444  info=8;
445  std::string s="Variable tol_abs too small, J singular, or bad ";
446  s+="scaling in mroot_cern::msolve().";
447  O2SCL_CONV_RET(s.c_str(),exc_efailed,this->err_nonconv);
448  }
449 
450  // Exit if necessary
451 
452  if (info!=0) {
453  O2SCL_ERR("Unspecified error in mroot_cern::msolve().",
454  exc_efailed);
455  }
456 
457  if (!((!lcv) || difit>0.05*xnorm)) {
458  // 8/20/08: Could this just be rewritten?
459  // if (lcv && difit<=0.05*xnorm)
460 
461  // Iterative refinement (if the iteration is converging)
462 
463  for(int m=2;m<=mopt && bskip==false;m++) {
464  fnorm1=fnorm;
465  fnorm=0.0;
466  for(k=0;k<((int)nvar) && bskip==false;k++) {
467  iflag=k;
468 
469  func(iflag,w1,f);
470 
471  fky=f[k];
472  nfcall++;
473  numf=(int)(((double)nfcall)/((double)nvar));
474 
475  if (fnorm<fabs(fky)) fnorm=fabs(fky);
476 
477  // Iterative refinement is terminated if it does not give a
478  // reduction on residuals
479 
480  if (fnorm>=fnorm1) {
481  fnorm=fnorm1;
482  bskip=true;
483  }
484 
485  if (!bskip) {
486  double temp3=fky/w0[k];
487 
488  for(i=0;i<((int)nvar);i++) {
489  w1[i]+=temp3*w(k,i);
490  }
491  }
492  }
493 
494  if (!bskip) {
495 
496  // Compute the norms of the iterate and correction vector
497 
498  xnorm=0.0;
499  difit=0.0;
500 
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]);
504  x[i]=w1[i];
505  }
506 
507  // Stopping criteria for iterative refinement
508 
509  if (fnorm<=this->tol_rel) info=1;
510  if (difit<=xnorm*this->tol_abs) info=2;
511  if (fnorm<=this->tol_rel && info==2) info=3;
512  if (numf>=lmaxf && info==0) {
513  info=4;
514  O2SCL_CONV_RET("Too many iterations in mroot_cern::msolve().",
515  exc_emaxiter,this->err_nonconv);
516  }
517 
518  if (info!=0) {
519  return 0;
520  }
521  }
522  }
523  }
524  }
525 
526  return 0;
527  }
528 
529 #ifndef DOXYGEN_INTERNAL
530 
531  protected:
532 
533  /// Internal storage for the value of \c info
534  int info;
535 
536  /// Store the number of function evaluations
537  int mpt[289];
538 
539 #endif
540 
541  };
542 
543 
544 }
545 
546 #endif
547 
iterative process is out of control
Definition: err_hnd.h:71
The main O<span style=&#39;position: relative; top: 0.3em; font-size: 0.8em&#39;>2</span>scl O$_2$scl names...
Definition: anneal.h:42
int verbose
Output control (default 0)
Definition: mroot.h:88
int get_info()
Get the value of INFO from the last call to msolve()
Definition: mroot_cern.h:166
#define O2SCL_CONV_RET(d, n, b)
Set a "convergence" error and return the error value.
Definition: err_hnd.h:292
double scale
The original scale parameter from CERNLIB (default 10.0)
Definition: mroot_cern.h:208
invalid argument supplied by user
Definition: err_hnd.h:59
std::function< int(size_t, const boost::numeric::ublas::vector< double > &, boost::numeric::ublas::vector< double > &) > mm_funct11
Array of multi-dimensional functions typedef.
Definition: mm_funct.h:43
apparent singularity detected
Definition: err_hnd.h:93
exceeded max number of iterations
Definition: err_hnd.h:73
int mpt[289]
Store the number of function evaluations.
Definition: mroot_cern.h:537
Multidimensional root-finding [abstract base].
Definition: mroot.h:66
generic failure
Definition: err_hnd.h:61
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.
Definition: mroot.h:133
virtual const char * type()
Return the type, "mroot_cern".
Definition: mroot_cern.h:204
iteration is not making progress toward solution
Definition: err_hnd.h:105
double eps
The smallest floating point number (default )
Definition: mroot_cern.h:230
bool err_nonconv
If true, call the error handler if msolve() or msolve_de() does not converge (default true) ...
Definition: mroot.h:99
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(). ...
Definition: mroot_cern.h:171
int maxf
Maximum number of function evaluations.
Definition: mroot_cern.h:201
virtual int msolve(size_t nvar, vec_t &x, func_t &func)
Solve func using x as an initial guess, returning x.
Definition: mroot_cern.h:233
std::function< int(size_t, boost::numeric::ublas::vector< double > &, size_t, boost::numeric::ublas::vector< double > &, boost::numeric::ublas::matrix< double > &) > jac_funct11
Jacobian function (not necessarily square)
Definition: jacobian.h:44
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
Multi-dimensional mroot-finding routine (CERNLIB)
Definition: mroot_cern.h:95
double tol_abs
The minimum allowable stepsize (default 1.0e-12)
Definition: mroot.h:85
double tol_rel
The maximum value of the functions for success (default 1.0e-8)
Definition: mroot.h:82
boost::numeric::ublas::matrix< double > w
Desc.
Definition: mroot_cern.h:103
std::string itos(int x)
Convert an integer to a string.
int info
Internal storage for the value of info.
Definition: mroot_cern.h:534

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).