Prev Next multi_newton.hpp Headings

OpenMP Multi-Threading Newton's Method Source Code


# 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


Input File: openmp/multi_newton.hpp