24 #ifndef MROOT_BROYDEN_H 25 #define MROOT_BROYDEN_H 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> 36 #include <o2scl/mroot.h> 37 #include <o2scl/permutation.h> 40 #ifndef DOXYGEN_NO_O2NS 60 typedef boost::numeric::ublas::matrix_column<ubmatrix> ubmatrix_column;
127 for(
size_t i=0;i<lu.size1();i++) {
128 for(
size_t j=0;j<lu.size2();j++) {
135 for(
size_t i=0;i<H.size1();i++) {
136 for(
size_t j=0;j<H.size2();j++) {
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;
155 def_jac.set_epsrel(sqrt(std::numeric_limits<double>::epsilon()));
189 double enorm(
size_t nvar,
const vec_t &ff) {
192 for (i=0;i<nvar;i++) {
204 void set(func_t &func,
size_t nvar, vec_t &x, vec_t &f, vec_t &dx) {
221 (*ajac)(nvar,x,nvar,f,
lu);
223 o2scl_linalg::LU_decomp<>(nvar,
lu,
perm,signum);
224 o2scl_linalg::LU_invert<ubmatrix,ubmatrix,ubmatrix_column>
227 for (i=0;i<nvar;i++) {
228 for (j=0;j<nvar;j++) {
233 for(
size_t i=0;i<dx.size();i++) dx[i]=0.0;
249 double phi0, phi1, t, lambda;
255 for (i=0;i<nvar;i++) {
257 for (j=0;j<nvar;j++) {
272 for (i=0;i<nvar;i++) {
273 x_trial[i]=x[i]+t*p[i];
277 int status=func(nvar,x_trial,fnew);
283 phi1=
enorm(nvar,fnew);
287 if (phi1>phi0 && iter<10 && t>0.1) {
290 double theta=phi1/phi0;
291 t*=(sqrt(1.0+6.0*theta)-1.0)/(3.0*theta);
295 if (new_step==
false && phi1>phi0) {
300 (*ajac)(nvar,x,nvar,f,
lu);
302 for (i=0;i<nvar;i++) {
303 for (j=0;j<nvar;j++) {
309 o2scl_linalg::LU_invert<ubmatrix,ubmatrix,ubmatrix_column>
315 for (i=0;i<nvar;i++) {
316 x_trial[i]=x[i]+t*p[i];
320 int status=func(nvar,x_trial,fnew);
326 phi1=
enorm(nvar,fnew);
332 for (i=0;i<nvar;i++) {
337 for (i=0;i<nvar;i++) {
339 for (j=0;j<nvar;j++) {
347 for (i=0;i<nvar;i++) {
351 O2SCL_ERR2(
"Approximation to Jacobian has collapsed in ",
356 for (i=0;i<nvar;i++) {
361 for (i=0;i<nvar;i++) {
363 for (j=0;j<nvar;j++) {
370 for (i=0;i<nvar;i++) {
372 for (j=0;j<nvar;j++) {
374 H(i,j)-=vi*wj/lambda;
383 for (i=0;i<nvar;i++) {
393 virtual int msolve(
size_t n, vec_t &x, func_t &func) {
413 for(
size_t i=0;i<n;i++) {
414 resid+=fabs(f_int[i]);
422 this->
print_iter(n,x,f_int,iter,resid,this->tol_rel,
430 if (iter>=this->ntrial) {
432 "exceeded max. number of iterations.",
444 #ifndef DOXYGEN_INTERNAL 451 (
const mroot_broyden<func_t,vec_t,mat_t,jfunc_t>&);
457 #ifndef DOXYGEN_NO_O2NS 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.
size_t user_nvar
Number of variables.
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)
size_t mem_size
Size of memory allocated.
virtual int set_function(func_t &f)
Set the function to compute the Jacobian of.
void resize(size_t dim)
Resize.
std::function< int(size_t, const boost::numeric::ublas::vector< double > &, boost::numeric::ublas::vector< double > &) > mm_funct
Array of multi-dimensional functions typedef.
vec_t dx_int
Stepsize vector.
vec_t * user_f
Function values.
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)
A class for representing permutations.
exceeded max number of iterations
jacobian< func_t, vec_t, mat_t > * ajac
Jacobian.
int LU_decomp(const size_t N, mat_t &A, o2scl::permutation &p, int &signum)
Compute the LU decomposition of the matrix A.
#define O2SCL_CONV2_RET(d, d2, n, b)
Set an error and return the error value, two-string version.
Multidimensional root-finding [abstract base].
iteration has not converged
virtual int msolve(size_t n, vec_t &x, func_t &func)
Desc.
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.
int ntrial
Maximum number of iterations (default 100)
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.
bool err_nonconv
If true, call the error handler if msolve() or msolve_de() does not converge (default true) ...
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
int last_ntrial
The number of iterations for in the most recent minimization.
vec_t f_int
Function value vector.
permutation perm
Permutation object for the LU decomposition.
problem with user-supplied function
Multidimensional root-finding using Broyden's method (GSL)
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.
double tol_rel
The maximum value of the functions for success (default 1.0e-8)
ubmatrix lu
LU decomposition.
Simple automatic Jacobian.
Base for providing a numerical jacobian [abstract base].