source: branches/cache/cppad/local/sqrt_op.hpp @ 3324

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

merge trunk changes into cache

  • Property svn:keywords set to Id
File size: 4.8 KB
Line 
1/* $Id: sqrt_op.hpp 3324 2014-09-12 12:14:53Z 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( q < cap_order );
47        CPPAD_ASSERT_UNKNOWN( p <= q );
48
49        // Taylor coefficients corresponding to argument and result
50        Base* x = taylor + i_x * cap_order;
51        Base* z = taylor + i_z * cap_order;
52
53        size_t k;
54        if( p == 0 )
55        {       z[0] = sqrt( x[0] );
56                p++;
57        }
58        for(size_t j = p; j <= q; j++)
59        {
60                CPPAD_ASSERT_KNOWN(
61                        x[0] != Base(0),
62                        "Forward: attempt to take derivatve of square root of zero"
63                )
64                z[j] = Base(0);
65                for(k = 1; k < j; k++)
66                        z[j] -= Base(k) * z[k] * z[j-k];
67                z[j] /= Base(j);
68                z[j] += x[j] / Base(2);
69                z[j] /= z[0];
70        }
71}
72
73/*!
74Multiple direction forward mode Taylor coefficient for op = SqrtOp.
75
76The C++ source code corresponding to this operation is
77\verbatim
78        z = sqrt(x)
79\endverbatim
80
81\copydetails forward_unary1_op_dir
82*/
83template <class Base>
84inline void forward_sqrt_op_dir(
85        size_t q           ,
86        size_t r           ,
87        size_t i_z         ,
88        size_t i_x         ,
89        size_t cap_order   , 
90        Base*  taylor      )
91{       
92        // check assumptions
93        CPPAD_ASSERT_UNKNOWN( NumArg(SqrtOp) == 1 );
94        CPPAD_ASSERT_UNKNOWN( NumRes(SqrtOp) == 1 );
95        CPPAD_ASSERT_UNKNOWN( 0 < q );
96        CPPAD_ASSERT_UNKNOWN( q < cap_order );
97
98        // Taylor coefficients corresponding to argument and result
99        size_t num_taylor_per_var = (cap_order-1) * r + 1;
100        Base* z = taylor + i_z * num_taylor_per_var;
101        Base* x = taylor + i_x * num_taylor_per_var;
102        CPPAD_ASSERT_KNOWN(
103                x[0] != Base(0),
104                "Forward: attempt to take derivatve of square root of zero"
105        )
106
107        size_t m = (q-1) * r + 1;
108        for(size_t ell = 0; ell < r; ell++)
109        {       z[m+ell] = Base(0);
110                for(size_t k = 1; k < q; k++)
111                        z[m+ell] -= Base(k) * z[(k-1)*r+1+ell] * z[(q-k-1)*r+1+ell];
112                z[m+ell] /= Base(q);
113                z[m+ell] += x[m+ell] / Base(2);
114                z[m+ell] /= z[0];
115        }       
116}
117
118/*!
119Compute zero order forward mode Taylor coefficient for result of op = SqrtOp.
120
121The C++ source code corresponding to this operation is
122\verbatim
123        z = sqrt(x)
124\endverbatim
125
126\copydetails forward_unary1_op_0
127*/
128template <class Base>
129inline void forward_sqrt_op_0(
130        size_t i_z         ,
131        size_t i_x         ,
132        size_t cap_order   , 
133        Base*  taylor      )
134{
135        // check assumptions
136        CPPAD_ASSERT_UNKNOWN( NumArg(SqrtOp) == 1 );
137        CPPAD_ASSERT_UNKNOWN( NumRes(SqrtOp) == 1 );
138        CPPAD_ASSERT_UNKNOWN( 0 < cap_order );
139
140        // Taylor coefficients corresponding to argument and result
141        Base* x = taylor + i_x * cap_order;
142        Base* z = taylor + i_z * cap_order;
143
144        z[0] = sqrt( x[0] );
145}
146/*!
147Compute reverse mode partial derivatives for result of op = SqrtOp.
148
149The C++ source code corresponding to this operation is
150\verbatim
151        z = sqrt(x)
152\endverbatim
153
154\copydetails reverse_unary1_op
155*/
156
157template <class Base>
158inline void reverse_sqrt_op(
159        size_t      d            ,
160        size_t      i_z          ,
161        size_t      i_x          ,
162        size_t      cap_order    , 
163        const Base* taylor       ,
164        size_t      nc_partial   ,
165        Base*       partial      )
166{
167        // check assumptions
168        CPPAD_ASSERT_UNKNOWN( NumArg(SqrtOp) == 1 );
169        CPPAD_ASSERT_UNKNOWN( NumRes(SqrtOp) == 1 );
170        CPPAD_ASSERT_UNKNOWN( d < cap_order );
171        CPPAD_ASSERT_UNKNOWN( d < nc_partial );
172
173        // Taylor coefficients and partials corresponding to argument
174        Base* px       = partial + i_x * nc_partial;
175
176        // Taylor coefficients and partials corresponding to result
177        const Base* z  = taylor  + i_z * cap_order;
178        Base* pz       = partial + i_z * nc_partial;
179
180        CPPAD_ASSERT_KNOWN(
181                z[0] != Base(0),
182                "Reverse: attempt to take derivatve of square root of zero"
183        )
184
185        // number of indices to access
186        size_t j = d;
187        size_t k;
188        while(j)
189        {
190
191                // scale partial w.r.t. z[j]
192                pz[j]   /= z[0];
193
194                pz[0]   -= pz[j] * z[j];
195                px[j]   += pz[j] / Base(2);
196                for(k = 1; k < j; k++)
197                        pz[k]   -= pz[j] * z[j-k];
198                --j;
199        }
200        px[0] += pz[0] / (Base(2) * z[0]);
201}
202
203} // END_CPPAD_NAMESPACE
204# endif
Note: See TracBrowser for help on using the repository browser.