source: trunk/cppad/local/sqrt_op.hpp @ 3301

Last change on this file since 3301 was 3301, checked in by bradbell, 6 years ago

merge in multiple forward direcitons from branches/forward_dir

  • Property svn:keywords set to Id
File size: 5.0 KB
Line 
1/* $Id: sqrt_op.hpp 3301 2014-05-24 05:20:21Z bradbell $ */
2# ifndef CPPAD_SQRT_OP_INCLUDED
3# define CPPAD_SQRT_OP_INCLUDED
4
5/* --------------------------------------------------------------------------
6CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-14 Bradley M. Bell
7
8CppAD is distributed under multiple licenses. This distribution is under
9the terms of the
10                    Eclipse Public License Version 1.0.
11
12A copy of this license is included in the COPYING file of this distribution.
13Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
14-------------------------------------------------------------------------- */
15
16
17namespace CppAD { // BEGIN_CPPAD_NAMESPACE
18/*!
19\file sqrt_op.hpp
20Forward and reverse mode calculations for z = sqrt(x).
21*/
22
23
24/*!
25Compute forward mode Taylor coefficient for result of op = SqrtOp.
26
27The C++ source code corresponding to this operation is
28\verbatim
29        z = sqrt(x)
30\endverbatim
31
32\copydetails forward_unary1_op
33*/
34template <class Base>
35inline void forward_sqrt_op(
36        size_t p           ,
37        size_t q           ,
38        size_t i_z         ,
39        size_t i_x         ,
40        size_t cap_order   , 
41        Base*  taylor      )
42{       
43        // check assumptions
44        CPPAD_ASSERT_UNKNOWN( NumArg(SqrtOp) == 1 );
45        CPPAD_ASSERT_UNKNOWN( NumRes(SqrtOp) == 1 );
46        CPPAD_ASSERT_UNKNOWN( i_x < i_z );
47        CPPAD_ASSERT_UNKNOWN( q < cap_order );
48        CPPAD_ASSERT_UNKNOWN( p <= q );
49
50        // Taylor coefficients corresponding to argument and result
51        Base* x = taylor + i_x * cap_order;
52        Base* z = taylor + i_z * cap_order;
53
54        size_t k;
55        if( p == 0 )
56        {       z[0] = sqrt( x[0] );
57                p++;
58        }
59        for(size_t j = p; j <= q; j++)
60        {
61                CPPAD_ASSERT_KNOWN(
62                        x[0] != Base(0),
63                        "Forward: attempt to take derivatve of square root of zero"
64                )
65                z[j] = Base(0);
66                for(k = 1; k < j; k++)
67                        z[j] -= Base(k) * z[k] * z[j-k];
68                z[j] /= Base(j);
69                z[j] += x[j] / Base(2);
70                z[j] /= z[0];
71        }
72}
73
74/*!
75Multiple direction forward mode Taylor coefficient for op = SqrtOp.
76
77The C++ source code corresponding to this operation is
78\verbatim
79        z = sqrt(x)
80\endverbatim
81
82\copydetails forward_unary1_op_dir
83*/
84template <class Base>
85inline void forward_sqrt_op_dir(
86        size_t q           ,
87        size_t r           ,
88        size_t i_z         ,
89        size_t i_x         ,
90        size_t cap_order   , 
91        Base*  taylor      )
92{       
93        // check assumptions
94        CPPAD_ASSERT_UNKNOWN( NumArg(SqrtOp) == 1 );
95        CPPAD_ASSERT_UNKNOWN( NumRes(SqrtOp) == 1 );
96        CPPAD_ASSERT_UNKNOWN( i_x < i_z );
97        CPPAD_ASSERT_UNKNOWN( 0 < q );
98        CPPAD_ASSERT_UNKNOWN( q < cap_order );
99
100        // Taylor coefficients corresponding to argument and result
101        size_t num_taylor_per_var = (cap_order-1) * r + 1;
102        Base* z = taylor + i_z * num_taylor_per_var;
103        Base* x = taylor + i_x * num_taylor_per_var;
104        CPPAD_ASSERT_KNOWN(
105                x[0] != Base(0),
106                "Forward: attempt to take derivatve of square root of zero"
107        )
108
109        size_t m = (q-1) * r + 1;
110        for(size_t ell = 0; ell < r; ell++)
111        {       z[m+ell] = Base(0);
112                for(size_t k = 1; k < q; k++)
113                        z[m+ell] -= Base(k) * z[(k-1)*r+1+ell] * z[(q-k-1)*r+1+ell];
114                z[m+ell] /= Base(q);
115                z[m+ell] += x[m+ell] / Base(2);
116                z[m+ell] /= z[0];
117        }       
118}
119
120/*!
121Compute zero order forward mode Taylor coefficient for result of op = SqrtOp.
122
123The C++ source code corresponding to this operation is
124\verbatim
125        z = sqrt(x)
126\endverbatim
127
128\copydetails forward_unary1_op_0
129*/
130template <class Base>
131inline void forward_sqrt_op_0(
132        size_t i_z         ,
133        size_t i_x         ,
134        size_t cap_order   , 
135        Base*  taylor      )
136{
137        // check assumptions
138        CPPAD_ASSERT_UNKNOWN( NumArg(SqrtOp) == 1 );
139        CPPAD_ASSERT_UNKNOWN( NumRes(SqrtOp) == 1 );
140        CPPAD_ASSERT_UNKNOWN( i_x < i_z );
141        CPPAD_ASSERT_UNKNOWN( 0 < cap_order );
142
143        // Taylor coefficients corresponding to argument and result
144        Base* x = taylor + i_x * cap_order;
145        Base* z = taylor + i_z * cap_order;
146
147        z[0] = sqrt( x[0] );
148}
149/*!
150Compute reverse mode partial derivatives for result of op = SqrtOp.
151
152The C++ source code corresponding to this operation is
153\verbatim
154        z = sqrt(x)
155\endverbatim
156
157\copydetails reverse_unary1_op
158*/
159
160template <class Base>
161inline void reverse_sqrt_op(
162        size_t      d            ,
163        size_t      i_z          ,
164        size_t      i_x          ,
165        size_t      cap_order    , 
166        const Base* taylor       ,
167        size_t      nc_partial   ,
168        Base*       partial      )
169{
170        // check assumptions
171        CPPAD_ASSERT_UNKNOWN( NumArg(SqrtOp) == 1 );
172        CPPAD_ASSERT_UNKNOWN( NumRes(SqrtOp) == 1 );
173        CPPAD_ASSERT_UNKNOWN( i_x < i_z );
174        CPPAD_ASSERT_UNKNOWN( d < cap_order );
175        CPPAD_ASSERT_UNKNOWN( d < nc_partial );
176
177        // Taylor coefficients and partials corresponding to argument
178        Base* px       = partial + i_x * nc_partial;
179
180        // Taylor coefficients and partials corresponding to result
181        const Base* z  = taylor  + i_z * cap_order;
182        Base* pz       = partial + i_z * nc_partial;
183
184        CPPAD_ASSERT_KNOWN(
185                z[0] != Base(0),
186                "Reverse: attempt to take derivatve of square root of zero"
187        )
188
189        // number of indices to access
190        size_t j = d;
191        size_t k;
192        while(j)
193        {
194
195                // scale partial w.r.t. z[j]
196                pz[j]   /= z[0];
197
198                pz[0]   -= pz[j] * z[j];
199                px[j]   += pz[j] / Base(2);
200                for(k = 1; k < j; k++)
201                        pz[k]   -= pz[j] * z[j-k];
202                --j;
203        }
204        px[0] += pz[0] / (Base(2) * z[0]);
205}
206
207} // END_CPPAD_NAMESPACE
208# endif
Note: See TracBrowser for help on using the repository browser.