![]() |
Prev | Next | multi_newton.hpp | Headings |
# include <cppad/cppad.hpp>
# include <cassert>
# ifdef _OPENMP
# include <omp.h>
# endif
namespace { // BEGIN CppAD namespace
template <class Fun>
void one_newton(double &fcur, double &xcur, Fun &fun,
double xlow, double xin, double xup, double epsilon, size_t max_itr)
{ using CppAD::AD;
using CppAD::vector;
using CppAD::abs;
// domain space vector
size_t n = 1;
vector< AD<double> > X(n);
// range space vector
size_t m = 1;
vector< AD<double> > Y(m);
// domain and range differentials
vector<double> dx(n), dy(m);
size_t itr;
xcur = xin;
for(itr = 0; itr < max_itr; itr++)
{ // domain space vector
X[0] = xcur;
CppAD::Independent(X);
// range space vector
Y[0] = fun(X[0]);
// F : X -> Y
CppAD::ADFun<double> F(X, Y);
// fcur = F(xcur)
fcur = Value(Y[0]);
// evaluate dfcur = F'(xcur)
dx[0] = 1;
dy = F.Forward(1, dx);
double dfcur = dy[0];
// check end of iterations
if( abs(fcur) <= epsilon )
return;
if( (xcur == xlow) & (fcur * dfcur > 0.) )
return;
if( (xcur == xup) & (fcur * dfcur < 0.) )
return;
if( dfcur == 0. )
return;
// next Newton iterate
double delta_x = - fcur / dfcur;
if( xlow - xcur >= delta_x )
xcur = xlow;
else if( xup - xcur <= delta_x )
xcur = xup;
else xcur = xcur + delta_x;
}
return;
}
template <class Fun>
void multi_newton(
CppAD::vector<double> &xout ,
Fun &fun ,
size_t n_grid ,
double xlow ,
double xup ,
double epsilon ,
size_t max_itr )
{ using CppAD::AD;
using CppAD::vector;
using CppAD::abs;
// check argument values
assert( xlow < xup );
assert( n_grid > 0 );
// OpenMP uses integers in place of size_t
int i, n = int(n_grid);
// set up grid
vector<double> grid(n_grid + 1);
vector<double> fcur(n_grid), xcur(n_grid), xmid(n_grid);
double dx = (xup - xlow) / double(n_grid);
for(i = 0; size_t(i) < n_grid; i++)
{ grid[i] = xlow + i * dx;
xmid[i] = xlow + (i + .5) * dx;
}
grid[n_grid] = xup;
# ifdef _OPENMP
# pragma omp parallel for
# endif
for(i = 0; i < n; i++)
{ one_newton(
fcur[i] ,
xcur[i] ,
fun ,
grid[i] ,
xmid[i] ,
grid[i+1] ,
epsilon ,
max_itr
);
}
// end omp parallel for
// remove duplicates and points that are not solutions
double xlast = xlow;
size_t ilast = 0;
size_t n_zero = 0;
for(i = 0; size_t(i) < n_grid; i++)
{ if( abs( fcur[i] ) <= epsilon )
{ if( n_zero == 0 )
{ xcur[n_zero++] = xlast = xcur[i];
ilast = i;
}
else if( fabs( xcur[i] - xlast ) > dx )
{ xcur[n_zero++] = xlast = xcur[i];
ilast = i;
}
else if( fabs( fcur[i] ) < fabs( fcur[ilast] ) )
{ xcur[n_zero - 1] = xlast = xcur[i];
ilast = i;
}
}
}
// resize output vector and set its values
xout.resize(n_zero);
for(i = 0; size_t(i) < n_zero; i++)
xout[i] = xcur[i];
return;
}
} // END CppAD namespace