Defines a User Atomic Operation that Computes Square Root

Syntax
atomic_user a_square_root  a_square_root(au, ay)

Purpose
This user atomic 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.cpp file
# include "multi_atomic.hpp"
//
namespace {

class atomic_user : public CppAD::atomic_base<double> {
public:
// ctor
atomic_user(void)
{ }
private:
// forward mode routine called by CppAD
virtual bool forward(
size_t                   p   ,
size_t                   q   ,
const vector<bool>&      vu  ,
vector<bool>&            vy  ,
const vector<double>&    tu  ,
vector<double>&          ty  )
{
# ifndef NDEBUG
size_t n = tu.size() / (q + 1);
size_t m = ty.size() / (q + 1);
assert( n == 3 );
assert( m == 1 );
# endif
// only implementing zero order forward for this example
if( q != 0 )
return false;

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

// check for setting variable information
if( vu.size() > 0 )
{     if( vu[0] )
return false;
vy[0] = vu[1] || vu[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
ty[0] = y_itr;
return true;
}
};
}