mroot_broyden.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2011-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 
24 #ifndef MROOT_BROYDEN_H
25 #define MROOT_BROYDEN_H
26 
27 /** \file mroot_broyden.h
28  \brief File defining \ref o2scl::mroot_broyden
29 */
30 
31 #include <boost/numeric/ublas/vector.hpp>
32 #include <boost/numeric/ublas/vector_proxy.hpp>
33 #include <boost/numeric/ublas/matrix.hpp>
34 #include <boost/numeric/ublas/matrix_proxy.hpp>
35 
36 #include <o2scl/mroot.h>
37 #include <o2scl/permutation.h>
38 #include <o2scl/lu.h>
39 
40 #ifndef DOXYGEN_NO_O2NS
41 namespace o2scl {
42 #endif
43 
44  /** \brief Multidimensional root-finding using Broyden's method (GSL)
45 
46  Experimental.
47 
48  See \ref Broyden65.
49  */
50  template<class func_t=mm_funct11,
53  class jfunc_t=jac_funct11>
54  class mroot_broyden : public mroot<func_t,vec_t,jfunc_t> {
55 
56  public:
57 
60  typedef boost::numeric::ublas::matrix_column<ubmatrix> ubmatrix_column;
61 
62  protected:
63 
64  /// Desc
65  ubmatrix H;
66 
67  /// LU decomposition
68  ubmatrix lu;
69 
70  /// Permutation object for the LU decomposition
72 
73  /// Desc
74  ubvector v;
75 
76  /// Desc
77  ubvector w;
78 
79  /// Desc
80  ubvector y;
81 
82  /// Desc
83  ubvector p;
84 
85  /// Desc
86  ubvector fnew;
87 
88  /// Desc
89  ubvector x_trial;
90 
91  /// Desc
92  double phi;
93 
94  /// Stepsize vector
95  vec_t dx_int;
96 
97  /// Function value vector
98  vec_t f_int;
99 
100  /// A pointer to the user-specified function
101  func_t *user_func;
102 
103  /// Function values
104  vec_t *user_f;
105 
106  /// Initial guess and current solution
107  vec_t *user_x;
108 
109  /// Initial and current step
110  vec_t *user_dx;
111 
112  /// Number of variables
113  size_t user_nvar;
114 
115  /// Size of memory allocated
116  size_t mem_size;
117 
118  /// Jacobian
120 
121  /** \brief Clear allocated vectors and matrices
122 
123  This function is called by set() before each solve.
124  */
125  void clear() {
126  if (mem_size>0) {
127  for(size_t i=0;i<lu.size1();i++) {
128  for(size_t j=0;j<lu.size2();j++) {
129  lu(i,j)=0.0;
130  }
131  }
132  for(size_t i=0;i<mem_size;i++) {
133  perm[i]=0;
134  }
135  for(size_t i=0;i<H.size1();i++) {
136  for(size_t j=0;j<H.size2();j++) {
137  H(i,j)=0.0;
138  }
139  }
140  for(size_t i=0;i<v.size();i++) v[i]=0.0;
141  for(size_t i=0;i<w.size();i++) w[i]=0.0;
142  for(size_t i=0;i<y.size();i++) y[i]=0.0;
143  for(size_t i=0;i<fnew.size();i++) fnew[i]=0.0;
144  for(size_t i=0;i<x_trial.size();i++) x_trial[i]=0.0;
145  for(size_t i=0;i<p.size();i++) p[i]=0.0;
146  }
147  return;
148  }
149 
150  public:
151 
152  mroot_broyden() {
153  mem_size=0;
154  ajac=&def_jac;
155  def_jac.set_epsrel(sqrt(std::numeric_limits<double>::epsilon()));
156  }
157 
158  /// Default Jacobian object
160 
161  virtual ~mroot_broyden() {
162  }
163 
164  /// Allocate memory
165  void allocate(size_t n) {
166 
167  if (n!=mem_size) {
168 
169  lu.resize(n,n);
170  perm.resize(n);
171  H.resize(n,n);
172  v.resize(n);
173  w.resize(n);
174  y.resize(n);
175  fnew.resize(n);
176  x_trial.resize(n);
177  p.resize(n);
178  dx_int.resize(n);
179  f_int.resize(n);
180 
181  mem_size=n;
182 
183  }
184 
185  return;
186  }
187 
188  /// Euclidean norm
189  double enorm(size_t nvar, const vec_t &ff) {
190  double e2=0;
191  size_t i;
192  for (i=0;i<nvar;i++) {
193  double fi=ff[i];
194  e2+=fi*fi;
195  }
196  return sqrt(e2);
197  }
198 
199  /** \brief Set the function, initial guess, and provide vectors to
200  store function values and stepsize
201 
202  The initial values of \c f and \c dx are ignored.
203  */
204  void set(func_t &func, size_t nvar, vec_t &x, vec_t &f, vec_t &dx) {
205 
206  if (nvar!=mem_size) allocate(nvar);
207  clear();
208 
209  user_func=&func;
210  user_nvar=nvar;
211  user_x=&x;
212  user_f=&f;
213  user_dx=&dx;
214 
215  size_t i, j;
216  int signum=0;
217 
218  func(nvar,x,f);
219 
220  ajac->set_function(func);
221  (*ajac)(nvar,x,nvar,f,lu);
222 
223  o2scl_linalg::LU_decomp<>(nvar,lu,perm,signum);
224  o2scl_linalg::LU_invert<ubmatrix,ubmatrix,ubmatrix_column>
225  (nvar,lu,perm,H);
226 
227  for (i=0;i<nvar;i++) {
228  for (j=0;j<nvar;j++) {
229  H(i,j)=-H(i,j);
230  }
231  }
232 
233  for(size_t i=0;i<dx.size();i++) dx[i]=0.0;
234 
235  phi=enorm(nvar,f);
236 
237  return;
238  }
239 
240  /// Perform an iteration
241  int iterate() {
242 
243  func_t &func=*user_func;
244  size_t nvar=user_nvar;
245  vec_t &x=*user_x;
246  vec_t &f=*user_f;
247  vec_t &dx=*user_dx;
248 
249  double phi0, phi1, t, lambda;
250 
251  size_t i, j, iter;
252 
253  /* p=H f */
254 
255  for (i=0;i<nvar;i++) {
256  double sum=0;
257  for (j=0;j<nvar;j++) {
258  sum+=H(i,j)*f[j];
259  }
260  p[i]=sum;
261  }
262 
263  t=1.0;
264  iter=0;
265 
266  phi0=phi;
267 
268  bool new_step=true;
269  while (new_step) {
270  new_step=false;
271 
272  for (i=0;i<nvar;i++) {
273  x_trial[i]=x[i]+t*p[i];
274  }
275 
276  {
277  int status=func(nvar,x_trial,fnew);
278  if (status != success) {
279  return exc_ebadfunc;
280  }
281  }
282 
283  phi1=enorm(nvar,fnew);
284 
285  iter++;
286 
287  if (phi1>phi0 && iter<10 && t>0.1) {
288 
289  // [GSL] Full step goes uphill, take a reduced step instead.
290  double theta=phi1/phi0;
291  t*=(sqrt(1.0+6.0*theta)-1.0)/(3.0*theta);
292  new_step=true;
293  }
294 
295  if (new_step==false && phi1>phi0) {
296 
297  // [GSL] Need to recompute Jacobian
298  int signum=0;
299 
300  (*ajac)(nvar,x,nvar,f,lu);
301 
302  for (i=0;i<nvar;i++) {
303  for (j=0;j<nvar;j++) {
304  lu(i,j)=-lu(i,j);
305  }
306  }
307 
308  o2scl_linalg::LU_decomp(nvar,lu,perm,signum);
309  o2scl_linalg::LU_invert<ubmatrix,ubmatrix,ubmatrix_column>
310  (nvar,lu,perm,H);
311  o2scl_linalg::LU_solve(nvar,lu,perm,f,p);
312 
313  t=1.0;
314 
315  for (i=0;i<nvar;i++) {
316  x_trial[i]=x[i]+t*p[i];
317  }
318 
319  {
320  int status=func(nvar,x_trial,fnew);
321  if (status != success) {
322  return exc_ebadfunc;
323  }
324  }
325 
326  phi1=enorm(nvar,fnew);
327  }
328 
329  }
330 
331  /* y=f'-f */
332  for (i=0;i<nvar;i++) {
333  y[i]=fnew[i]-f[i];
334  }
335 
336  /* v=H y */
337  for (i=0;i<nvar;i++) {
338  double sum=0.0;
339  for (j=0;j<nvar;j++) {
340  sum+=H(i,j)*y[j];
341  }
342  v[i]=sum;
343  }
344 
345  /* lambda=p dot v */
346  lambda=0.0;
347  for (i=0;i<nvar;i++) {
348  lambda+=p[i]*v[i];
349  }
350  if (lambda==0) {
351  O2SCL_ERR2("Approximation to Jacobian has collapsed in ",
352  "mroot_broyden::iterate().",exc_ezerodiv);
353  }
354 
355  /* v'=v+t*p */
356  for (i=0;i<nvar;i++) {
357  v[i]=v[i]+t*p[i];
358  }
359 
360  /* w^T=p^T H */
361  for (i=0;i<nvar;i++) {
362  double sum=0;
363  for (j=0;j<nvar;j++) {
364  sum+=H(j,i)*p[j];
365  }
366  w[i]=sum;
367  }
368 
369  /* Hij -> Hij-(vi wj/lambda) */
370  for (i=0;i<nvar;i++) {
371  double vi=v[i];
372  for (j=0;j<nvar;j++) {
373  double wj=w[j];
374  H(i,j)-=vi*wj/lambda;
375  }
376  }
377 
378  /* copy fnew into f */
379  vector_copy(nvar,fnew,f);
380 
381  /* copy x_trial into x */
382  vector_copy(nvar,x_trial,x);
383  for (i=0;i<nvar;i++) {
384  dx[i]=t*p[i];
385  }
386 
387  phi=phi1;
388 
389  return success;
390  }
391 
392  /// Desc
393  virtual int msolve(size_t n, vec_t &x, func_t &func) {
394 
395  int status;
396 
397  set(func,n,x,f_int,dx_int);
398 
399  int iter=0;
400 
401  do {
402  iter++;
403 
404  status=iterate();
405  if (status) break;
406 
407  // ------------------------------------------------------
408  // The equivalent of the statement:
409  //
410  // status=gsl_multiroot_test_residual(f,this->tol_rel);
411 
412  double resid=0.0;
413  for(size_t i=0;i<n;i++) {
414  resid+=fabs(f_int[i]);
415  }
416  if (resid<this->tol_rel) status=success;
417  else status=gsl_continue;
418 
419  // ------------------------------------------------------
420 
421  if (this->verbose>0) {
422  this->print_iter(n,x,f_int,iter,resid,this->tol_rel,
423  "mroot_broyden");
424  }
425 
426  } while (status==gsl_continue && iter<this->ntrial);
427 
428  this->last_ntrial=iter;
429 
430  if (iter>=this->ntrial) {
431  O2SCL_CONV2_RET("Function mroot_broyden::msolve() ",
432  "exceeded max. number of iterations.",
433  exc_emaxiter,this->err_nonconv);
434  }
435 
436  if (status!=0) {
437  O2SCL_ERR2("Function iterate() failed in ",
438  "mroot_broyden::solve_set().",exc_efailed);
439  }
440 
441  return success;
442  }
443 
444 #ifndef DOXYGEN_INTERNAL
445 
446  private:
447 
451  (const mroot_broyden<func_t,vec_t,mat_t,jfunc_t>&);
452 
453 #endif
454 
455  };
456 
457 #ifndef DOXYGEN_NO_O2NS
458 }
459 #endif
460 
461 #endif
int iterate()
Perform an iteration.
void allocate(size_t n)
Allocate memory.
void clear()
Clear allocated vectors and matrices.
jacobian_gsl< func_t, vec_t, mat_t > def_jac
Default Jacobian object.
vec_t * user_x
Initial guess and current solution.
func_t * user_func
A pointer to the user-specified function.
tried to divide by zero
Definition: err_hnd.h:75
size_t user_nvar
Number of variables.
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
size_t mem_size
Size of memory allocated.
virtual int set_function(func_t &f)
Set the function to compute the Jacobian of.
Definition: jacobian.h:78
void resize(size_t dim)
Resize.
Definition: permutation.h:242
vec_t dx_int
Stepsize vector.
Definition: mroot_broyden.h:95
vec_t * user_f
Function values.
ubvector p
Desc.
Definition: mroot_broyden.h:83
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
A class for representing permutations.
Definition: permutation.h:70
exceeded max number of iterations
Definition: err_hnd.h:73
jacobian< func_t, vec_t, mat_t > * ajac
Jacobian.
ubvector v
Desc.
Definition: mroot_broyden.h:74
int LU_decomp(const size_t N, mat_t &A, o2scl::permutation &p, int &signum)
Compute the LU decomposition of the matrix A.
Definition: lu_base.h:86
#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
iteration has not converged
Definition: err_hnd.h:51
virtual int msolve(size_t n, vec_t &x, func_t &func)
Desc.
ubvector x_trial
Desc.
Definition: mroot_broyden.h:89
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
int ntrial
Maximum number of iterations (default 100)
Definition: mroot.h:91
vec_t * user_dx
Initial and current step.
double enorm(size_t nvar, const vec_t &ff)
Euclidean norm.
void vector_copy(const vec_t &src, vec2_t &dest)
Simple vector copy.
Definition: vector.h:127
bool err_nonconv
If true, call the error handler if msolve() or msolve_de() does not converge (default true) ...
Definition: mroot.h:99
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
int last_ntrial
The number of iterations for in the most recent minimization.
Definition: mroot.h:94
vec_t f_int
Function value vector.
Definition: mroot_broyden.h:98
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
ubvector y
Desc.
Definition: mroot_broyden.h:80
permutation perm
Permutation object for the LU decomposition.
Definition: mroot_broyden.h:71
problem with user-supplied function
Definition: err_hnd.h:69
Multidimensional root-finding using Broyden&#39;s method (GSL)
Definition: mroot_broyden.h:54
int LU_solve(const size_t N, const mat_t &LU, const o2scl::permutation &p, const vec_t &b, vec2_t &x)
Solve a linear system after LU decomposition.
Definition: lu_base.h:248
double tol_rel
The maximum value of the functions for success (default 1.0e-8)
Definition: mroot.h:82
ubmatrix lu
LU decomposition.
Definition: mroot_broyden.h:68
ubvector w
Desc.
Definition: mroot_broyden.h:77
Success.
Definition: err_hnd.h:47
ubmatrix H
Desc.
Definition: mroot_broyden.h:65
Simple automatic Jacobian.
Definition: jacobian.h:144
Base for providing a numerical jacobian [abstract base].
Definition: jacobian.h:64
ubvector fnew
Desc.
Definition: mroot_broyden.h:86

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