mroot_hybrids.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 /* multiroots/hybrid.c
24  * multiroots/dogleg.c
25  *
26  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Brian Gough
27  *
28  * This program is free software; you can redistribute it and/or modify
29  * it under the terms of the GNU General Public License as published by
30  * the Free Software Foundation; either version 3 of the License, or (at
31  * your option) any later version.
32  *
33  * This program is distributed in the hope that it will be useful, but
34  * WITHOUT ANY WARRANTY; without even the implied warranty of
35  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
36  * General Public License for more details.
37  *
38  * You should have received a copy of the GNU General Public License
39  * along with this program; if not, write to the Free Software
40  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
41  * 02110-1301, USA.
42  */
43 /** \file mroot_hybrids.h
44  \brief File defining \ref o2scl::mroot_hybrids and
45  specializations
46 */
47 #ifndef O2SCL_MROOT_HYBRIDS_H
48 #define O2SCL_MROOT_HYBRIDS_H
49 
50 #include <string>
51 
52 #include <gsl/gsl_multiroots.h>
53 #include <gsl/gsl_linalg.h>
54 
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>
59 
60 #include <o2scl/vector.h>
61 #include <o2scl/mroot.h>
62 #include <o2scl/jacobian.h>
63 #include <o2scl/qr.h>
64 // For matrix_out() below
65 #include <o2scl/columnify.h>
66 
67 #ifndef DOXYGEN_NO_O2NS
68 namespace o2scl {
69 #endif
70 
71  /** \brief Base functions for \ref mroot_hybrids
72  */
74 
75  public:
76 
79 
80  /// Compute the actual reduction
81  double compute_actual_reduction(double fnorm0, double fnorm1) {
82  double actred;
83  if (fnorm1<fnorm0) {
84  double u=fnorm1/fnorm0;
85  actred=1-u*u;
86  } else {
87  actred=-1;
88  }
89  return actred;
90  }
91 
92  /// Compute the predicted reduction phi1p=|Q^T f + R dx|
93  double compute_predicted_reduction(double fnorm0, double fnorm1) {
94  double prered;
95  if (fnorm1<fnorm0) {
96  double u=fnorm1/fnorm0;
97  prered=1-u*u;
98  } else {
99  prered=0;
100  }
101  return prered;
102  }
103 
104  /// Compute \f$ R \cdot g \f$ where \c g is the \gradient
105  template<class vec2_t, class mat_t>
106  void compute_Rg(size_t N, const mat_t &r2,
107  const ubvector &gradient2, vec2_t &Rg) {
108 
109  for (size_t i=0;i<N;i++) {
110  double sum=0.0;
111  for (size_t j=i;j<N;j++) {
112  sum+=r2(i,j)*gradient2[j];
113  }
114  Rg[i]=sum;
115  }
116 
117  return;
118  }
119 
120  /// Compute \c w and \c v
121  template<class vec2_t>
122  void compute_wv(size_t n, const ubvector &qtdf2,
123  const ubvector &rdx2, const vec2_t &dx2,
124  const ubvector &diag2, double pnorm,
125  ubvector &w2, ubvector &v2) {
126 
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;
131  }
132 
133  return;
134  }
135 
136  /// Compute \f$ R \cdot \mathrm{dx} \f$
137  template<class vec2_t, class mat_t>
138  void compute_rdx(size_t N, const mat_t &r2,
139  const vec2_t &dx2, ubvector &rdx2) {
140 
141  for (size_t i=0;i<N;i++) {
142  double sum=0.0;
143  for (size_t j=i;j<N;j++) {
144  sum+=r2(i,j)*dx2[j];
145  }
146  rdx2[i]=sum;
147  }
148 
149  return;
150  }
151 
152  /** \brief Compute the norm of the vector \f$ \vec{v} \f$ defined
153  by \f$ v_i = d_i ff_i \f$
154  */
155  template<class vec2_t, class vec3_t>
156  double scaled_enorm(size_t n, const vec2_t &d, const vec3_t &ff) {
157  double e2=0.0;
158  for (size_t i=0;i<n;i++) {
159  double u=d[i]*ff[i];
160  e2+=u*u;
161  }
162  return sqrt(e2);
163  }
164 
165  /// Compute delta
166  template<class vec2_t>
167  double compute_delta(size_t n, ubvector &diag2, vec2_t &x2) {
168  double Dx=scaled_enorm(n,diag2,x2);
169  double factor=100.0;
170  return (Dx > 0) ? factor*Dx : factor;
171  }
172 
173  /** \brief Compute the Euclidean norm of \c ff
174 
175  \future Replace this with \c dnrm2 from \ref cblas_base.h
176  */
177  template<class vec2_t> double enorm(size_t N, const vec2_t &ff) {
178  double e2=0.0;
179  for (size_t i=0;i<N;i++) {
180  double fi=ff[i];
181  e2+=fi*fi;
182  }
183  return sqrt(e2);
184  }
185 
186  /// Compute the Euclidean norm of the sum of \c a and \c b
187  double enorm_sum(size_t n, const ubvector &a,
188  const ubvector &b) {
189  double e2=0.0;
190  for (size_t i=0;i<n;i++) {
191  double u=a[i]+b[i];
192  e2+=u*u;
193  }
194  return sqrt(e2);
195  }
196 
197  /** \brief Compute a trial step and put the result in \c xx_trial
198 
199  \future Replace with daxpy.
200  */
201  template<class vec2_t>
202  int compute_trial_step(size_t N, vec2_t &xl, vec2_t &dxl,
203  vec2_t &xx_trial) {
204  for (size_t i=0;i<N;i++) {
205  xx_trial[i]=xl[i]+dxl[i];
206  }
207  return 0;
208  }
209 
210  /** \brief Compute the change in the function value
211 
212  \future Replace with daxpy?
213  */
214  template<class vec2_t>
215  int compute_df(size_t n, const vec2_t &ff_trial,
216  const vec2_t &fl, ubvector &dfl) {
217  for (size_t i=0;i<n;i++) {
218  dfl[i]=ff_trial[i]-fl[i];
219  }
220 
221  return 0;
222  }
223 
224  /// Compute \c diag, the norm of the columns of the Jacobian
225  template<class mat2_t>
226  void compute_diag(size_t n, const mat2_t &J2, ubvector &diag2) {
227  for (size_t j=0;j<n;j++) {
228  double sum=0;
229  for (size_t i=0;i<n;i++) {
230  const double Jij=J2(i,j);
231  sum+=Jij*Jij;
232  }
233  if (sum == 0) {
234  sum=1.0;
235  }
236  diag2[j]=sqrt(sum);
237  }
238 
239  return;
240  }
241 
242  /** \brief Compute \f$ Q^{T} f \f$
243 
244  \future This is just right-multiplication, so we could
245  use the \o2 cblas routines instead.
246  */
247  template<class vec2_t, class vec3_t, class vec4_t>
248  void compute_qtf(size_t N, const vec2_t &q2,
249  const vec3_t &f2, vec4_t &qtf2) {
250  for (size_t j=0;j<N;j++) {
251  double sum=0.0;
252  for (size_t i=0;i<N;i++) {
253  sum+=q2(i,j)*f2[i];
254  }
255  qtf2[j]=sum;
256  }
257 
258  return;
259  }
260 
261  /// Update \c diag
262  template<class mat2_t>
263  void update_diag(size_t n, const mat2_t &J2, ubvector &diag2) {
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++) {
267  double Jij=J2(i,j);
268  sum+=Jij*Jij;
269  }
270  if (sum == 0) {
271  sum=1.0;
272  }
273  cnorm=sqrt(sum);
274  diagj=diag2[j];
275  if (cnorm > diagj) {
276  diag2[j]=cnorm;
277  }
278  }
279 
280  return;
281  }
282 
283  /** \brief Form appropriate convex combination of the Gauss-Newton
284  direction and the scaled gradient direction.
285 
286  Using the Gauss-Newton direction given in \c newton (a vector of
287  size N), and the gradient direction in \c gradient (also a
288  vector of size N), this computes
289  \f[
290  \mathrm{pp}=\alpha \mathrm{newton}+\beta \mathrm{gradient}
291  \f]
292  */
293  template<class vec2_t>
294  void scaled_addition(size_t N, double alpha, ubvector &newton2,
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];
298  }
299 
300  return;
301  }
302 
303  /// Compute the Gauss-Newton direction
304  template<class mat_t>
305  int newton_direction(const size_t N, const mat_t &r2,
306  const ubvector &qtf2, ubvector &p) {
307  size_t i;
308  int status=0;
309 
310  vector_copy(N,qtf2,p);
311  o2scl_cblas::dtrsv(o2scl_cblas::o2cblas_RowMajor,
312  o2scl_cblas::o2cblas_Upper,
313  o2scl_cblas::o2cblas_NoTrans,
314  o2scl_cblas::o2cblas_NonUnit,N,N,r2,p);
315 
316  for (i=0;i<N;i++) {
317  p[i]*=-1.0;
318  }
319 
320  return status;
321  }
322 
323  /// Compute the gradient direction
324  template<class mat_t>
325  void gradient_direction(const size_t M, const size_t N,
326  const mat_t &r2, const ubvector &qtf2,
327  const ubvector &diag2, ubvector &g) {
328  for (size_t j=0;j<M;j++) {
329  double sum=0.0;
330  for (size_t i=0;i<N;i++) {
331  sum+=r2(i,j)*qtf2[i];
332  }
333  g[j]=-sum/diag2[j];
334  }
335 
336  return;
337  }
338 
339  /// Compute the point at which the gradient is minimized
340  void minimum_step(const size_t N, double gnorm,
341  const ubvector &diag2, ubvector &g) {
342  for (size_t i=0;i<N;i++) {
343  g[i]=(g[i]/gnorm)/diag2[i];
344  }
345 
346  return;
347  }
348 
349  /** \brief Take a dogleg step
350 
351  Given the QR decomposition of an \c n by \c n matrix "A",
352  a diagonal matrix \c diag, a vector \c b, and a positive
353  number \c delta, this function determines the convex
354  combination \c x of the Gauss-Newton and scaled gradient
355  directions which minimizes \f$ A x-b \f$ in the least
356  squares sense, subject to the restriction that the
357  Euclidean norm of \f$ d x \f$ is smaller than the
358  value of \c delta.
359  */
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,
364  vec2_t &p) {
365 
366  double qnorm, gnorm, sgnorm, bnorm, temp;
367 
368  // Compute the Gauss-Newton direction
369  newton_direction(n,r2,qtf2,newton2);
370  qnorm=scaled_enorm(n,diag2,newton2);
371 
372  // Test whether the gauss-newton direction is acceptable
373  if (qnorm <= delta2) {
374  for(size_t i=0;i<n;i++) p[i]=newton2[i];
375  return success;
376  }
377 
378  // Compute the gradient direction
379  gradient_direction(n,n,r2,qtf2,diag2,gradient2);
380  gnorm=enorm(n,gradient2);
381 
382  // The special case that the scaled gradient is zero
383  if (gnorm == 0) {
384  double alpha=delta2/qnorm;
385  double beta=0;
386  scaled_addition(n,alpha,newton2,beta,gradient2,p);
387  return success;
388  }
389 
390  // -----------------------------------------------------
391  // The next set of statements calculates the point along the
392  // scaled gradient at which the quadratic is minimized.
393 
394  minimum_step(n,gnorm,diag2,gradient2);
395 
396  // Use p as temporary space to compute Rg
397  compute_Rg(n,r2,gradient2,p);
398 
399  temp=enorm(n,p);
400  sgnorm=(gnorm/temp)/temp;
401 
402  // -----------------------------------------------------
403  // Test whether the scaled gradient direction is acceptable
404 
405  if (sgnorm>delta2) {
406 
407  double alpha=0;
408  double beta=delta2;
409 
410  scaled_addition(n,alpha,newton2,beta,gradient2,p);
411 
412  return success;
413  }
414 
415  // If not, calculate the point along the dogleg step
416  // at which the quadratic is minimized
417  bnorm=enorm(n,qtf2);
418 
419  double bg=bnorm/gnorm;
420  double bq=bnorm/qnorm;
421  double dq=delta2/qnorm;
422  double dq2=dq*dq;
423  double sd=sgnorm/delta2;
424  double sd2=sd*sd;
425 
426  double t1=bg*bq*sd;
427  double u=t1-dq;
428  double t2=t1-dq*sd2+sqrt(u*u+(1-dq2)*(1-sd2));
429 
430  double alpha=dq*(1-sd2)/t2;
431  double beta=(1-alpha)*sgnorm;
432 
433  // Form the appropriate convex combination of the gauss-newton
434  // direction and the scaled gradient direction.
435 
436  scaled_addition(n,alpha,newton2,beta,gradient2,p);
437 
438  return success;
439  }
440 
441  };
442 
443  /** \brief Multidimensional root-finding algorithm using
444  Powell's Hybrid method (GSL)
445 
446  This is a recasted version of the GSL routines which use a
447  modified version of Powell's Hybrid method as implemented in the
448  HYBRJ algorithm in MINPACK (\ref Garbow80). Both the scaled and
449  unscaled options are available by setting \ref int_scaling (the
450  scaled version is the default). If derivatives are not provided,
451  they will be computed automatically. This class provides the
452  GSL-like interface using allocate(), set() (or set_de() in case
453  where derivatives are available), iterate(), and the
454  higher-level interfaces, msolve() and msolve_de(), which perform
455  the solution and the memory allocation automatically. Some
456  additional checking is performed in case the user calls the
457  functions out of order (i.e. set() without allocate()).
458 
459  The functions msolve() and msolve_de() use the condition \f$
460  \sum_i |f_i|<\f$ \ref mroot::tol_rel to determine if the solver has
461  succeeded.
462 
463  See the \ref multisolve_subsect section of the User's guide for
464  general information about \o2 solvers.
465  There is an example for the usage of the multidimensional solver
466  classes given in <tt>examples/ex_mroot.cpp</tt>, see the \ref
467  ex_mroot_sect .
468 
469  \note The \ref set() and \ref set_de() functions store a pointer
470  to the function object and the user must ensure that the object
471  is still valid for a later call to \ref iterate().
472 
473  The original GSL algorithm has been modified to shrink the
474  stepsize if a proposed step causes the function to return a
475  non-zero value. This allows the routine to automatically try to
476  avoid regions where the function is not defined. The algorithm
477  is also modified to check that it is not sending non-finite
478  values to the user-specified function. To return to the default
479  GSL behavior, set \ref shrink_step and \ref extra_finite_check
480  to false.
481 
482  The default method for numerically computing the Jacobian is
483  from \ref jacobian_gsl. This default is identical to the GSL
484  approach, except that the default value of \ref
485  jacobian_gsl::epsmin is non-zero. See \ref jacobian_gsl
486  for more details.
487 
488  By default convergence failures result in calling the exception
489  handler, but this can be turned off by setting \ref
490  mroot::err_nonconv to false. If \ref mroot::err_nonconv is
491  false, then the functions \ref iterate(), \ref msolve() and \ref
492  msolve_de() will return a non-zero value if convergence fails.
493  Note that the default Jacobian object, \ref def_jac also has a
494  data member \ref jacobian_gsl::err_nonconv which separately
495  handles the case where the one row of the Jacobian is all zero.
496 
497  \future Is all the setting of vectors and matrices to zero really
498  necessary? Do they need to be executed even if memory hasn't
499  been recently allocated?
500 
501  \future Convert more ubvectors to vec_t.
502 
503  \future Some more of the element-wise vector manipulation could be
504  converted to BLAS routines.
505 
506  \future It's kind of strange that set() sets jac_given to false
507  and set_de() has to reset it to true. Can this be simplified?
508 
509  \future Many of these minpack functions could be put in their
510  own "minpack_tools" class, or possibly moved to be
511  linear algebra routines instead.
512 
513  \future There are still some numbers in here which the user
514  could have control over, for example, the <tt>nslow2</tt>
515  threshold which indicates failure.
516  */
517  template<
518  class func_t=mm_funct11,
521  class jfunc_t=jac_funct11 > class mroot_hybrids :
522  public mroot<func_t,vec_t,jfunc_t>, mroot_hybrids_base {
523 
524 #ifndef DOXYGEN_INTERNAL
525 
526  protected:
527 
528  /// Number of iterations
529  int iter;
530  /// Compute the number of failures
531  size_t ncfail;
532  /// Compute the number of successes
533  size_t ncsuc;
534  /// The number of times the actual reduction is less than 0.001
535  size_t nslow1;
536  /// The number of times the actual reduction is less than 0.1
537  size_t nslow2;
538  /// The norm of the current function value
539  double fnorm;
540  /// The limit of the Nuclidean norm
541  double delta;
542  /// Jacobian
543  mat_t J;
544  /// Q matrix from QR decomposition
545  mat_t q;
546  /// R matrix from QR decomposition
547  mat_t r;
548  /// The diagonal elements
550  /// The value of \f$ Q^T f \f$
552  /// The Newton direction
554  /// The gradient direction
556  /// The change in the function value
558  /// The value of \f$ Q^T \cdot \mathrm{df} \f$
560  /// The value of \f$ R \cdot \mathrm{dx} \f$
562  /// The value of \f$ w=(Q^T df-R dx)/|dx| \f$
564  /// The value of \f$ v=D^2 dx/|dx| \f$
566 
567  /// The user-specified Jacobian
568  jfunc_t *jac;
569 
570  /// The automatic Jacobian
572 
573  /// The value of the derivative
574  vec_t dx;
575 
576  /// Trial root
577  vec_t x_trial;
578 
579  /// Trial function value
580  vec_t f_trial;
581 
582  /// The number of equations and unknowns
583  size_t dim;
584 
585  /// True if the jacobian has been given
586  bool jac_given;
587 
588  /// The user-specified function
589  func_t *fnewp;
590 
591  /// True if "set" has been called
593 
594  /// Finish the solution after set() or set_de() has been called
595  virtual int solve_set(size_t nn, vec_t &xx, func_t &ufunc) {
596 
597  int status;
598  iter=0;
599 
600  do {
601  iter++;
602 
603  if (iterate()!=0) {
604  O2SCL_CONV2_RET("Function iterate() failed in mroot_hybrids::",
605  "solve_set().",exc_efailed,this->err_nonconv);
606  }
607 
608  // ------------------------------------------------------
609  // The equivalent of the statement:
610  //
611  // status=gsl_multiroot_test_residual(f,this->tol_rel);
612 
613  double resid=0.0;
614  for(size_t i=0;i<nn;i++) {
615  resid+=fabs(f[i]);
616  }
617  if (resid<this->tol_rel) status=success;
618  else status=gsl_continue;
619 
620  // ------------------------------------------------------
621 
622  if (this->verbose>0) {
623  this->print_iter(nn,x,f,iter,resid,this->tol_rel,
624  "mroot_hybrids");
625  }
626 
627  } while (status==gsl_continue && iter<this->ntrial);
628 
629  for(size_t i=0;i<nn;i++) {
630  xx[i]=x[i];
631  }
632 
633  this->last_ntrial=iter;
634 
635  if (((int)iter)>=this->ntrial) {
636  O2SCL_CONV2_RET("Function mroot_hybrids::msolve() ",
637  "exceeded max. number of iterations.",
638  exc_emaxiter,this->err_nonconv);
639  }
640 
641  return success;
642  }
643 
644 #endif
645 
646  public:
647 
648  mroot_hybrids() {
649  shrink_step=true;
650  dim=0;
651  ajac=&def_jac;
652  //def_jac.set_epsrel(sqrt(std::numeric_limits<double>::epsilon()));
653  int_scaling=true;
654  jac_given=false;
655  set_called=false;
656  extra_finite_check=true;
657  }
658 
659  virtual ~mroot_hybrids() {
660  }
661 
662  /** \brief If true, iterate() will shrink the step-size automatically if
663  the function returns a non-zero value (default true)
664 
665  The original GSL behavior can be obtained by setting
666  this to \c false.
667  */
669 
670  /** \brief If true, double check that the input function values are
671  finite (default true)
672  */
674 
675  /// If true, use the internal scaling method (default true)
677 
678  /// Default automatic Jacobian object
680 
681  /// Set the automatic Jacobian object
683  ajac=&j;
684  return 0;
685  }
686 
687  /** \brief The value of the function at the present iteration
688 
689  \comment
690  We need this to be public so that the user can see if
691  iterate() is working
692  \endcomment
693  */
694  vec_t f;
695 
696  /// The present solution
697  vec_t x;
698 
699  /** \brief Perform an iteration
700 
701  At the end of the iteration, the current value of the solution
702  is stored in \ref x.
703  */
704  int iterate() {
705 
706  if (!set_called) {
707  O2SCL_ERR2("Function set() not called or most recent call ",
708  "failed in mroot_hybrids::iterate().",
709  exc_efailed);
710  }
711 
712  double prered, actred;
713  double pnorm, fnorm1, fnorm1p;
714  double ratio;
715  double p1=0.1, p5=0.5, p001=0.001, p0001=0.0001;
716 
717  /* Compute qtf=Q^T f */
718 
719  compute_qtf(dim,q,f,qtf);
720 
721  /* Compute dogleg step */
722 
723  dogleg(dim,r,qtf,diag,delta,newton,gradient,dx);
724 
725  /* Take a trial step */
726 
727  compute_trial_step(dim,x,dx,x_trial);
728  pnorm=scaled_enorm(dim,diag,dx);
729 
730  if (iter==1) {
731  if (pnorm<delta) {
732  delta=pnorm;
733  }
734  }
735 
736  if (extra_finite_check) {
737  // AWS, 11/14/13: This is not strictly necessary, because if
738  // x_trial is not finite the iteration will fail to converge
739  // anyway, but it is disconcerting to the user to have
740  // non-finite values sent to the user-specified function for
741  // no apparent reason. In reality what appears to be happening
742  // is that pnorm was previously zero (because of a vanishingly
743  // small step size), and the call to compute_wv() below then
744  // leads to non-finite values. On the other hand, checking all
745  // the vector values is time consuming, so I perform this
746  // check only if extra_finite_check is true.
747  for(size_t ik=0;ik<dim;ik++) {
748  if (!std::isfinite(x_trial[ik])) {
749  O2SCL_CONV2_RET("Iteration lead to non-finite values in ",
750  "mroot_hybrids::iterate().",exc_efailed,
751  this->err_nonconv);
752  }
753  }
754  }
755 
756  /* Evaluate function at x+p */
757 
758  int status;
759 
760  if (shrink_step==false) {
761 
762  // Evaluate the function at the new point, exit if it fails
763 
764  status=(*fnewp)(dim,x_trial,f_trial);
765 
766  if (status != success) {
767  std::string str="Function returned non-zero value ("+itos(status)+
768  ") in mroot_hybrids::iterate().";
769  O2SCL_CONV_RET(str.c_str(),o2scl::exc_ebadfunc,this->err_nonconv);
770  }
771 
772  } else {
773 
774  // Evaluate the function at the new point, try to recover
775  // if it fails
776 
777  status=(*fnewp)(dim,x_trial,f_trial);
778 
779  int bit=0;
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;
783  }
784  status=(*fnewp)(dim,x_trial,f_trial);
785  bit++;
786  }
787 
788  // Exit if we weren't able to find a new good point
789 
790  if (status != success) {
791  std::string str="No suitable step found, function returned ("+
792  itos(status)+") in mroot_hybrids::iterate().";
793  O2SCL_CONV_RET(str.c_str(),o2scl::exc_ebadfunc,this->err_nonconv);
794  }
795 
796  }
797 
798  /* Set df=f_trial-f */
799 
800  compute_df(dim,f_trial,f,df);
801 
802  /* Compute the scaled actual reduction */
803 
804  fnorm1=enorm(dim,f_trial);
805  actred=compute_actual_reduction(fnorm,fnorm1);
806 
807  /* Compute rdx=R dx */
808 
809  compute_rdx(dim,r,dx,rdx);
810 
811  /* Compute the scaled predicted reduction phi1p=|Q^T f+R dx| */
812 
813  fnorm1p=enorm_sum(dim,qtf,rdx);
814  prered=compute_predicted_reduction(fnorm,fnorm1p);
815 
816  /* Compute the ratio of the actual to predicted reduction */
817 
818  if (prered > 0) {
819  ratio=actred/prered;
820  } else {
821  ratio=0;
822  }
823 
824  /* Update the step bound */
825 
826  if (ratio<p1) {
827  ncsuc=0;
828  ncfail++;
829  delta*=p5;
830  } else {
831  ncfail=0;
832  ncsuc++;
833 
834  if (ratio >= p5 || ncsuc > 1) {
835  delta=GSL_MAX(delta,pnorm/p5);
836  }
837  if (fabs (ratio-1) <= p1) {
838  delta=pnorm/p5;
839  }
840  }
841 
842  /* Test for successful iteration */
843 
844  if (ratio >= p0001) {
845  for(size_t i=0;i<dim;i++) {
846  x[i]=x_trial[i];
847  f[i]=f_trial[i];
848  }
849  fnorm=fnorm1;
850  iter++;
851  }
852 
853  /* Determine the progress of the iteration */
854 
855  nslow1++;
856  if (actred >= p001) {
857  nslow1=0;
858  }
859 
860  if (actred >= p1) {
861  nslow2=0;
862  }
863 
864  if (ncfail==2) {
865 
866  int jac_ret;
867 
868  if (jac_given) jac_ret=(*jac)(dim,x,dim,f,J);
869  else jac_ret=(*ajac)(dim,x,dim,f,J);
870 
871  if (jac_ret!=0) {
872  std::string str="Jacobian failed and returned ("+
873  itos(jac_ret)+") in mroot_hybrids::iterate().";
874  O2SCL_CONV_RET(str.c_str(),exc_efailed,this->err_nonconv);
875  }
876 
877  nslow2++;
878 
879  if (iter==1) {
880  if (int_scaling) {
881  compute_diag(dim,J,diag);
882  }
883  delta=compute_delta(dim,diag,x);
884  } else {
885  if (int_scaling) {
886  update_diag(dim,J,diag);
887  }
888  }
889 
890  o2scl_linalg::QR_decomp_unpack(dim,dim,this->J,this->q,this->r);
891 
892  return success;
893  }
894 
895  /* Compute qtdf=Q^T df, w=(Q^T df-R dx)/|dx|, v=D^2 dx/|dx| */
896 
897  compute_qtf(dim,q,df,qtdf);
898  compute_wv(dim,qtdf,rdx,dx,diag,pnorm,w,v);
899 
900  /* Rank-1 update of the jacobian Q'R'=Q(R+w v^T) */
901 
902  o2scl_linalg::QR_update(dim,dim,q,r,w,v);
903 
904  /* No progress as measured by jacobian evaluations */
905 
906  if (nslow2==5) {
907  O2SCL_CONV2_RET("No progress in Jacobian in mroot_hybrids::",
908  "iterate().",exc_enoprogj,this->err_nonconv);
909  }
910 
911  /* No progress as measured by function evaluations */
912 
913  if (nslow1==10) {
914  O2SCL_CONV2_RET("No progress in function in mroot_hybrids::",
915  "iterate().",exc_enoprog,this->err_nonconv);
916  }
917 
918  return success;
919  }
920 
921  /// Allocate the memory
922  void allocate(size_t n) {
923 
924  q.resize(n,n);
925  r.resize(n,n);
926  diag.resize(n);
927  qtf.resize(n);
928  newton.resize(n);
929  gradient.resize(n);
930  df.resize(n);
931  qtdf.resize(n);
932  rdx.resize(n);
933  w.resize(n);
934  v.resize(n);
935 
936  for(size_t ii=0;ii<n;ii++) {
937  for(size_t jj=0;jj<n;jj++) {
938  q(ii,jj)=0.0;
939  }
940  }
941  for(size_t ii=0;ii<n;ii++) {
942  for(size_t jj=0;jj<n;jj++) {
943  r(ii,jj)=0.0;
944  }
945  }
946 
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;
956 
957  J.resize(n,n);
958 
959  x.resize(n);
960  dx.resize(n);
961  f.resize(n);
962  x_trial.resize(n);
963  f_trial.resize(n);
964 
965  for(size_t i=0;i<n;i++) {
966  x[i]=0.0;
967  dx[i]=0.0;
968  f[i]=0.0;
969  }
970 
971  dim=n;
972 
973  return;
974  }
975 
976  /// Return the type,\c "mroot_hybrids".
977  virtual const char *type() { return "mroot_hybrids"; }
978 
979  /** \brief Solve \c func with derivatives \c dfunc using \c x as
980  an initial guess, returning \c x.
981 
982  */
983  virtual int msolve_de(size_t nn, vec_t &xx, func_t &ufunc,
984  jfunc_t &dfunc) {
985 
986  int ret=set_de(nn,xx,ufunc,dfunc);
987  if (ret!=success) {
988  O2SCL_CONV2_RET("Function set_de() failed in mroot_hybrids::",
989  "msolve_de().",exc_efailed,this->err_nonconv);
990  }
991 
992  return solve_set(nn,xx,ufunc);
993  }
994 
995  /// Solve \c ufunc using \c xx as an initial guess, returning \c xx.
996  virtual int msolve(size_t nn, vec_t &xx, func_t &ufunc) {
997 
998  int ret=set(nn,xx,ufunc);
999  if (ret!=success) {
1000  O2SCL_CONV2_RET("Function set() failed in mroot_hybrids::",
1001  "msolve().",exc_efailed,this->err_nonconv);
1002  }
1003  int status=solve_set(nn,xx,ufunc);
1004  for(size_t i=0;i<nn;i++) xx[i]=x[i];
1005 
1006  return status;
1007  }
1008 
1009  /** \brief Set the function, the parameters, and the initial guess
1010  */
1011  int set(size_t nn, vec_t &ax, func_t &ufunc) {
1012 
1013  int status;
1014 
1015  if (nn!=dim) {
1016  allocate(nn);
1017  }
1018 
1019  fnewp=&ufunc;
1020 
1021  // Specify function for automatic Jacobian
1022  ajac->set_function(ufunc);
1023 
1024  if (dim==0) {
1025  O2SCL_ERR2("No memory allocated in ",
1026  "mroot_hybrids::set().",o2scl::exc_ebadlen);
1027  }
1028 
1029  // Copy the user-specified solution
1030  for(size_t i=0;i<dim;i++) x[i]=ax[i];
1031 
1032  status=ufunc(dim,ax,f);
1033  if (status!=0) {
1034  O2SCL_ERR2("Function returned non-zero in ",
1035  "mroot_hybrids::set().",exc_ebadfunc);
1036  }
1037 
1038  if (jac_given) status=(*jac)(dim,ax,dim,f,J);
1039  else status=(*ajac)(dim,ax,dim,f,J);
1040 
1041  if (status!=0) {
1042  O2SCL_CONV2_RET("Jacobian failed in ",
1043  "mroot_hybrids::set().",exc_efailed,this->err_nonconv);
1044  }
1045 
1046  iter=1;
1047  fnorm=enorm(dim,f);
1048  ncfail=0;
1049  ncsuc=0;
1050  nslow1=0;
1051  nslow2=0;
1052 
1053  for(size_t i=0;i<dim;i++) dx[i]=0.0;
1054 
1055  /* Store column norms in diag */
1056 
1057  if (int_scaling) {
1058  compute_diag(dim,J,diag);
1059  } else {
1060  for(size_t ii=0;ii<diag.size();ii++) diag[ii]=0.0;
1061  }
1062 
1063  /* Set delta to factor |D x| or to factor if |D x| is zero */
1064 
1065  delta=compute_delta(dim,diag,x);
1066 
1067  /* Factorize J into QR decomposition */
1068  o2scl_linalg::QR_decomp_unpack(dim,dim,this->J,this->q,this->r);
1069  set_called=true;
1070  jac_given=false;
1071 
1072  return 0;
1073  }
1074 
1075  /** \brief Set the function, the Jacobian, the parameters,
1076  and the initial guess
1077  */
1078  int set_de(size_t nn, vec_t &ax, func_t &ufunc, jfunc_t &dfunc) {
1079 
1080  fnewp=&ufunc;
1081  jac=&dfunc;
1082 
1083  // Make sure set() uses the right Jacobian
1084  jac_given=true;
1085 
1086  int ret=set(nn,ax,ufunc);
1087 
1088  // Reset jac_given since set() will set it back to false
1089  jac_given=true;
1090  set_called=true;
1091 
1092  return ret;
1093  }
1094 
1095 #ifndef DOXYGEN_INTERNAL
1096 
1097  private:
1098 
1102  (const mroot_hybrids<func_t,vec_t,mat_t,jfunc_t>&);
1103 
1104 #endif
1105 
1106  };
1107 
1108 #ifndef DOXYGEN_NO_O2NS
1109 }
1110 #endif
1111 
1112 #if defined (O2SCL_COND_FLAG) || defined (DOXYGEN)
1113 #if defined (O2SCL_ARMA) || defined (DOXYGEN)
1114 #include <armadillo>
1115 namespace o2scl {
1116  /** \brief A version of \ref mroot_hybrids which uses
1117  Armadillo for the QR decomposition
1118 
1119  \note This class template is only defined if Armadillo
1120  was enabled when \o2 was installed.
1121  */
1122  template<class func_t, class vec_t, class mat_t, class jfunc_t>
1124  public mroot_hybrids<func_t,vec_t,mat_t,jfunc_t> {
1125 
1126  virtual void qr_decomp_unpack() {
1127  arma::qr_econ(this->q,this->r,this->J);
1128  return;
1129  }
1130 
1131  };
1132 }
1133 #endif
1134 #if defined (O2SCL_HAVE_EIGEN) || defined (DOXYGEN)
1135 #include <eigen3/Eigen/Dense>
1136 namespace o2scl {
1137  /** \brief A version of \ref mroot_hybrids
1138  which uses Eigen for the QR decomposition
1139 
1140  \note This class template is only defined if Eigen
1141  was enabled when \o2 was installed.
1142  */
1143  template<class func_t, class vec_t, class mat_t, class jfunc_t>
1145  public mroot_hybrids<func_t,vec_t,mat_t,jfunc_t> {
1146 
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>();
1151  return;
1152  }
1153 
1154  };
1155 }
1156 #endif
1157 #else
1158 #include <o2scl/mroot_special.h>
1159 #endif
1160 
1161 #endif
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=&#39;position: relative; top: 0.3em; font-size: 0.8em&#39;>2</span>scl O$_2$scl names...
Definition: anneal.h:42
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 .
Definition: cblas_base.h:308
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.
Definition: jacobian.h:78
vec_t x_trial
Trial root.
ubvector v
The value of .
ubvector w
The value of .
#define O2SCL_CONV_RET(d, n, b)
Set a "convergence" error and return the error value.
Definition: err_hnd.h:292
jacobian_gsl< func_t, vec_t, mat_t > def_jac
Default automatic Jacobian object.
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.
mat_t J
Jacobian.
void compute_Rg(size_t N, const mat_t &r2, const ubvector &gradient2, vec2_t &Rg)
Compute where g is the gradient gradient .
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
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|.
Definition: mroot_hybrids.h:93
exceeded max number of iterations
Definition: err_hnd.h:73
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.
Definition: err_hnd.h:298
Multidimensional root-finding [abstract base].
Definition: mroot.h:66
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
Definition: err_hnd.h:51
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 .
generic failure
Definition: err_hnd.h:61
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.
Definition: qr_base.h:231
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&#39;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.
Definition: vector.h:127
iteration is not making progress toward solution
Definition: err_hnd.h:105
matrix, vector lengths are not conformant
Definition: err_hnd.h:89
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.
Definition: err_hnd.h:281
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
Definition: err_hnd.h:107
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
Base functions for mroot_hybrids.
Definition: mroot_hybrids.h:73
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&#39; = A + u v^T,.
Definition: qr_base.h:165
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
Definition: err_hnd.h:69
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]
Definition: inte_qng_gsl.h:66
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.
Definition: mroot_hybrids.h:81
Success.
Definition: err_hnd.h:47
void compute_qtf(size_t N, const vec2_t &q2, const vec3_t &f2, vec4_t &qtf2)
Compute .
Simple automatic Jacobian.
Definition: jacobian.h:144
Base for providing a numerical jacobian [abstract base].
Definition: jacobian.h:64
ubvector rdx
The value of .

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