![]() |
Prev | Next | multi_newton_worker |
multi_newton_worker()
[ low , up ]
.
low = work_all_[thread_num]->xlow
up = work_all_[thread_num]->xup
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; } }