\(\newcommand{\W}[1]{ \; #1 \; }\) \(\newcommand{\R}[1]{ {\rm #1} }\) \(\newcommand{\B}[1]{ {\bf #1} }\) \(\newcommand{\D}[2]{ \frac{\partial #1}{\partial #2} }\) \(\newcommand{\DD}[3]{ \frac{\partial^2 #1}{\partial #2 \partial #3} }\) \(\newcommand{\Dpow}[2]{ \frac{\partial^{#1}}{\partial {#2}^{#1}} }\) \(\newcommand{\dpow}[2]{ \frac{ {\rm d}^{#1}}{{\rm d}\, {#2}^{#1}} }\)

multi_newton_worker

View page source

Do One Thread’s Work for Multi-Threaded Newton Method

Syntax

multi_newton_worker ()

Purpose

This function finds all the zeros in the interval [ low , up ] .

low

This is the value of the multi_newton_common information

low = work_all_ [ thread_num ] ->xlow

up

This is the value of the multi_newton_common information

up = work_all_ [ thread_num ] ->xup

thread_num

This is the number for the current thread; see thread_num .

Source

namespace {
void multi_newton_worker(void)
{
   // Split [xlow, xup] into num_sub intervales and
   // look for one zero in each sub-interval.
   size_t thread_num    = thread_alloc::thread_num();
   size_t num_threads   = std::max(num_threads_, size_t(1));
   bool   ok            = thread_num < num_threads;
   size_t num_sub       = work_all_[thread_num]->num_sub;
   double xlow          = work_all_[thread_num]->xlow;
   double xup           = work_all_[thread_num]->xup;
   vector<double>& x    = work_all_[thread_num]->x;

   // check arguments
   ok &= max_itr_ > 0;
   ok &= num_sub > 0;
   ok &= xlow < xup;
   ok &= x.size() == 0;

   // check for special case where there is nothing for this thread to do
   if( num_sub == 0 )
   {  work_all_[thread_num]->ok = ok;
      return;
   }

   // check for a zero on each sub-interval
   size_t i;
   double xlast = xlow - 2.0 * sub_length_; // over sub_length_ away from x_low
   double flast = 2.0 * epsilon_;           // any value > epsilon_ would do
   for(i = 0; i < num_sub; i++)
   {
      // note that when i == 0, xlow_i == xlow (exactly)
      double xlow_i = xlow + double(i) * sub_length_;

      // note that when i == num_sub - 1, xup_i = xup (exactly)
      double xup_i  = xup  - double(num_sub - i - 1) * sub_length_;

      // initial point for Newton iterations
      double xcur = (xup_i + xlow_i) / 2.;

      // Newton iterations
      bool more_itr = true;
      size_t itr    = 0;
      // initialize these values to avoid MSC C++ warning
      double fcur=0.0, dfcur=0.0;
      while( more_itr )
      {  fun_(xcur, fcur, dfcur);

         // check end of iterations
         if( fabs(fcur) <= epsilon_ )
            more_itr = false;
         if( (xcur == xlow_i ) & (fcur * dfcur > 0.) )
            more_itr = false;
         if( (xcur == xup_i)   & (fcur * dfcur < 0.) )
            more_itr = false;

         // next Newton iterate
         if( more_itr )
         {  xcur = xcur - fcur / dfcur;
            // keep in bounds
            xcur = std::max(xcur, xlow_i);
            xcur = std::min(xcur, xup_i);

            more_itr = ++itr < max_itr_;
         }
      }
      if( fabs( fcur ) <= epsilon_ )
      {  // check for case where xcur is lower bound for this
         // sub-interval and upper bound for previous sub-interval
         if( fabs(xcur - xlast) >= sub_length_ )
         {  x.push_back( xcur );
            xlast = xcur;
            flast = fcur;
         }
         else if( fabs(fcur) < fabs(flast) )
         {  x[ x.size() - 1] = xcur;
            xlast            = xcur;
            flast            = fcur;
         }
      }
   }
   work_all_[thread_num]->ok = ok;
}
}