multi_atomic_three_user¶

Defines a atomic_three Operation that Computes Square Root¶

Syntax¶

atomic_user a_square_root

a_square_root ( au , ay )

Purpose¶

This atomic function operation computes a square root using Newton’s method. It is meant to be very inefficient in order to demonstrate timing results.

au¶

This argument has prototype

const ADvector & au

where ADvector is a simple vector class with elements of type AD<double> . The size of au is three.

num_itr¶

We use the notation

num_itr = size_t ( Integer ( au [0] ) )

for the number of Newton iterations in the computation of the square root function. The component au [0] must be a Parameter .

y_initial¶

We use the notation

y_initial = au [1]

for the initial value of the Newton iterate.

y_squared¶

We use the notation

y_squared = au [2]

for the value we are taking the square root of.

ay¶

This argument has prototype

ADvector & ay

The size of ay is one and ay [0] is the square root of y_squared .

Limitations¶

Only zero order forward mode is implements for the atomic_user class.

Source¶

// includes used by all source code in multi_atomic_three.cpp file
# include <cppad/cppad.hpp>
# include "multi_atomic_three.hpp"
# include "team_thread.hpp"
//
namespace {
using CppAD::thread_alloc;                // multi-threading memory allocator
using CppAD::vector;                      // uses thread_alloc
typedef CppAD::ad_type_enum ad_type_enum; // constant, dynamic or variable

class atomic_user : public CppAD::atomic_three<double> {
public:
   // ctor
   atomic_user(void)
   : CppAD::atomic_three<double>("atomic_square_root")
   { }
private:
   // for_type
   bool for_type(
      const vector<double>&        parameter_u ,
      const vector<ad_type_enum>&  type_u      ,
      vector<ad_type_enum>&        type_y      ) override
   {  bool ok = parameter_u.size() == 3;
      ok     &= type_u.size() == 3;
      ok     &= type_y.size() == 1;
      if( ! ok )
         return false;
      ok     &= type_u[0] < CppAD::variable_enum;
      if( ! ok )
         return false;
      type_y[0] = std::max( type_u[0], type_u[1] );
      type_y[0] = std::max( type_y[0], type_u[2] );
      //
      return true;
   }
   // forward
   bool forward(
      const vector<double>&        parameter_u ,
      const vector<ad_type_enum>&  type_u      ,
      size_t                       need_y      ,
      size_t                       order_low   ,
      size_t                       order_up    ,
      const vector<double>&        taylor_u    ,
      vector<double>&              taylor_y    ) override
   {
# ifndef NDEBUG
      size_t n = taylor_u.size() / (order_up + 1);
      size_t m = taylor_y.size() / (order_up + 1);
      assert( n == 3 );
      assert( m == 1 );
# endif
      // only implementing zero order forward for this example
      if( order_up != 0 )
         return false;

      // extract components of argument vector
      size_t num_itr    = size_t( taylor_u[0] );
      double y_initial  = taylor_u[1];
      double y_squared  = taylor_u[2];

      // Use Newton's method to solve f(y) = y^2 = y_squared
      double y_itr = y_initial;
      for(size_t itr = 0; itr < num_itr; itr++)
      {  // solve (y - y_itr) * f'(y_itr) = y_squared - y_itr^2
         double fp_itr = 2.0 * y_itr;
         y_itr         = y_itr + (y_squared - y_itr * y_itr) / fp_itr;
      }

      // return the Newton approximation for f(y) = y_squared
      taylor_y[0] = y_itr;
      return true;
   }
};
}