Prev Next

Multi-Threaded Newton's Method Main Program

Syntax
multi_newton n_thread repeat n_zero n_grid n_sum

Purpose
Runs a timing test of the multi_newton routine. This routine uses Newton's method to determine if there is a zero of a function on each of a set of sub-intervals. CppAD is used to calculate the derivatives required by Newton's method. OpenMP is used to parallelize the calculation on the different sub-intervals.

n_thread
If the argument n_thread is equal to automatic, dynamic thread adjustment is used. Otherwise, n_thread must be a positive number specifying the number of OpenMP threads to use.

repeat
If the argument repeat is equal to automatic, the number of times to repeat the calculation of the number of zeros in total interval is automatically determined. In this case, the rate of execution of the total solution is reported.

If the argument repeat is not equal to automatic, it must be a positive integer. In this case repeat determination of the number of times the calculation of the zeros in the total interval is repeated. The rate of execution is not reported (it is assumed that the program execution time is being calculated some other way).

n_zero
The argument n_zero is the actual number of zeros that there should be in the test function,  \sin(x) . It must be an integer greater than one. The total interval searched for zeros is  [ 0 , (n\_zero - 1) \pi ] .

n_grid
The argument n_grid specifies the number of sub-intervals to divide the total interval into. It must an integer greater than zero (should probably be greater than two times n_zero).

n_sum
The actual function that is used is  \[
     f(x) = \frac{1}{n\_sum} \sum_{i=1}^{n\_sum} \sin (x)
\] 
where n_sum is a positive integer. The larger the value of n_sum, the more computation is required to calculate the function.

Subroutines
multi_newton Multi-Threaded Newton's Method Routine
multi_newton.hpp OpenMP Multi-Threading Newton's Method Source Code

Example Source


# include <cppad/cppad.hpp>
# include <cmath>
# include <cstring>
# include "multi_newton.hpp"

# ifdef _OPENMP
# include <omp.h>
# endif


namespace { // empty namespace
     size_t n_sum;  // larger values make fun(x) take longer to calculate
        size_t n_zero; // number of zeros of fun(x) in the total interval
}

// A slow version of the sine function
CppAD::AD<double> fun(const CppAD::AD<double> &x)
{    CppAD::AD<double> sum = 0.;
     size_t i;
     for(i = 0; i < n_sum; i++)
          sum += sin(x);

     return sum / double(n_sum);
}

void test_once(CppAD::vector<double> &xout, size_t n_grid)
{    assert( n_zero > 1 );
     double pi      = 4. * std::atan(1.); 
     double xlow    = 0.;
     double xup     = (n_zero - 1) * pi;
     double epsilon = 1e-6;
     size_t max_itr = 20;

     multi_newton(
          xout    ,
          fun     ,
          n_grid  ,
          xlow    ,
          xup     ,
          epsilon ,
          max_itr
     );
     return;
}

void test_repeat(size_t size, size_t repeat)
{    size_t i;
     CppAD::vector<double> xout;
     for(i = 0; i < repeat; i++)
          test_once(xout, size);
     return;
}

int main(int argc, char *argv[])
{
     using std::cout;
     using std::endl;
     using CppAD::vector;

     char *usage = "multi_newton n_thread repeat n_zero n_grid n_sum";
     if( argc != 6 )
     {    std::cerr << usage << endl;
          exit(1);
     }
     argv++;

     // n_thread command line argument
     int n_thread;
     if( std::strcmp(*argv, "automatic") == 0 )
          n_thread = 0;
     else n_thread = std::atoi(*argv);
     argv++;

     // repeat command line argument
     size_t repeat;
     if( std::strcmp(*argv, "automatic") == 0 )
          repeat = 0;
     else
     {    assert( std::atoi(*argv) > 0 );
          repeat = std::atoi(*argv);
     }
     argv++;

     // n_zero command line argument 
     assert( std::atoi(*argv) > 1 );
     n_zero = std::atoi(*argv++);

     // n_grid command line argument
     assert( std::atoi(*argv) > 0 );
     size_t n_grid = std::atoi(*argv++);
       
     // n_sum command line argument 
     assert( std::atoi(*argv) > 0 );
     n_sum = std::atoi(*argv++);

     // minimum time for test (repeat until this much time)
     double time_min = 1.;

# ifdef _OPENMP
     if( n_thread > 0 )
     {    omp_set_dynamic(0);            // off dynamic thread adjust
          omp_set_num_threads(n_thread); // set the number of threads 
     }
     // now determine the maximum number of threads
     n_thread = omp_get_max_threads();
     assert( n_thread > 0 );
     
     // No tapes are currently active,
     // so we can inform CppAD of the maximum number of threads
     CppAD::AD<double>::omp_max_thread(size_t(n_thread));

     // inform the user of the maximum number of threads
     cout << "OpenMP: version = "         << _OPENMP;
     cout << ", max number of threads = " << n_thread << endl;
# else
     cout << "_OPENMP is not defined, ";
     cout << "running in single tread mode" << endl;
     n_thread = 1;
# endif
     // initialize flag
     bool ok = true;

     // sub-block so xout gets deallocated before call to CppADTrackCount
     {    // Correctness check
          vector<double> xout;
          test_once(xout, n_grid);
          double epsilon = 1e-6;
          double pi      = 4. * std::atan(1.);
          ok            &= (xout.size() == n_zero);
          size_t i       = 0;
          while( ok & (i < n_zero) )
          {    ok &= std::fabs( xout[i] - pi * i) <= 2 * epsilon;
               ++i;
          }
     }
     if( repeat > 0 )
     {    // run the calculation the requested number of time
          test_repeat(n_grid, repeat);
     }
     else
     {    // actually time the calculation    

          // size of the one test case
          vector<size_t> size_vec(1);
          size_vec[0] = n_grid;

          // run the test case
          vector<size_t> rate_vec =
          CppAD::speed_test(test_repeat, size_vec, time_min);

          // report results
          cout << "n_grid           = " << size_vec[0] << endl;
          cout << "repeats per sec  = " << rate_vec[0] << endl;
     }
     // check all the threads for a CppAD memory leak
     if( CppADTrackCount() != 0 )
     {    ok = false;
          cout << "Error: memory leak detected" << endl;
     }
     if( ok )
          cout << "Correctness Test Passed" << endl;
     else cout << "Correctness Test Failed" << endl;

     return static_cast<int>( ! ok );
}


Input File: openmp/multi_newton.cpp