multi_atomic_two_user¶

Defines a atomic_two 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_two.cpp file
# include <cppad/cppad.hpp>
# include "multi_atomic_two.hpp"
# include "team_thread.hpp"
//
namespace {
using CppAD::thread_alloc; // fast multi-threading memory allocator
using CppAD::vector;       // uses thread_alloc

class atomic_user : public CppAD::atomic_base<double> {
public:
    // ctor
    atomic_user(void)
    : CppAD::atomic_base<double>("atomic_square_root")
    { }
private:
    // forward mode routine called by CppAD
    bool forward(
        size_t                   p   ,
        size_t                   q   ,
        const vector<bool>&      vu  ,
        vector<bool>&            vy  ,
        const vector<double>&    tu  ,
        vector<double>&          ty  ) override
    {
# 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;
    }
};
}