source: trunk/cppad/local/atan_op.hpp @ 3320

Last change on this file since 3320 was 3320, checked in by bradbell, 6 years ago
  1. g++ 4.8.2 has shadow warnings by default, but eigen and fadbad do not

these warnings, so supress then in these cases.

  1. Move check that arguments come before result into on place,

CPPAD_ASSERT_ARG_BEFORE_RESULT (only one argument case so far).

main.cpp: fix shadowing of index variable.
CMakeLists.txt: adapt to change in teuchos library name.
sparse_jacobian.cpp: fix a shadowed variable.
check_svn_id.sh: ignore svn_commit.sh.
gpl_license.sh: ignore svn_commit.sh.

  • Property svn:keywords set to Id
File size: 5.8 KB
Line 
1/* $Id: atan_op.hpp 3320 2014-09-11 23:06:21Z bradbell $ */
2# ifndef CPPAD_ATAN_OP_INCLUDED
3# define CPPAD_ATAN_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 atan_op.hpp
20Forward and reverse mode calculations for z = atan(x).
21*/
22
23
24/*!
25Forward mode Taylor coefficient for result of op = AtanOp.
26
27The C++ source code corresponding to this operation is
28\verbatim
29        z = atan(x)
30\endverbatim
31The auxillary result is
32\verbatim
33        y = 1 + x * x
34\endverbatim
35The value of y, and its derivatives, are computed along with the value
36and derivatives of z.
37
38\copydetails forward_unary2_op
39*/
40template <class Base>
41inline void forward_atan_op(
42        size_t p           ,
43        size_t q           ,
44        size_t i_z         ,
45        size_t i_x         ,
46        size_t cap_order   , 
47        Base*  taylor      )
48{       
49        // check assumptions
50        CPPAD_ASSERT_UNKNOWN( NumArg(AtanOp) == 1 );
51        CPPAD_ASSERT_UNKNOWN( NumRes(AtanOp) == 2 );
52        CPPAD_ASSERT_UNKNOWN( q < cap_order );
53        CPPAD_ASSERT_UNKNOWN( p <= q );
54
55        // Taylor coefficients corresponding to argument and result
56        Base* x = taylor + i_x * cap_order;
57        Base* z = taylor + i_z * cap_order;
58        Base* b = z      -       cap_order;  // called y in documentation
59
60        size_t k;
61        if( p == 0 )
62        {       z[0] = atan( x[0] );
63                b[0] = Base(1) + x[0] * x[0];
64                p++;
65        }
66        for(size_t j = p; j <= q; j++)
67        {
68                b[j] = Base(2) * x[0] * x[j];
69                z[j] = Base(0);
70                for(k = 1; k < j; k++)
71                {       b[j] += x[k] * x[j-k];
72                        z[j] -= Base(k) * z[k] * b[j-k];
73                }
74                z[j] /= Base(j);
75                z[j] += x[j];
76                z[j] /= b[0];
77        }
78}
79
80/*!
81Multiple direction Taylor coefficient for op = AtanOp.
82
83The C++ source code corresponding to this operation is
84\verbatim
85        z = atan(x)
86\endverbatim
87The auxillary result is
88\verbatim
89        y = 1 + x * x
90\endverbatim
91The value of y, and its derivatives, are computed along with the value
92and derivatives of z.
93
94\copydetails forward_unary2_op_dir
95*/
96template <class Base>
97inline void forward_atan_op_dir(
98        size_t q           ,
99        size_t r           ,
100        size_t i_z         ,
101        size_t i_x         ,
102        size_t cap_order   , 
103        Base*  taylor      )
104{       
105        // check assumptions
106        CPPAD_ASSERT_UNKNOWN( NumArg(AtanOp) == 1 );
107        CPPAD_ASSERT_UNKNOWN( NumRes(AtanOp) == 2 );
108        CPPAD_ASSERT_UNKNOWN( 0 < q );
109        CPPAD_ASSERT_UNKNOWN( q < cap_order );
110
111        // Taylor coefficients corresponding to argument and result
112        size_t num_taylor_per_var = (cap_order-1) * r + 1;
113        Base* x = taylor + i_x * num_taylor_per_var;
114        Base* z = taylor + i_z * num_taylor_per_var;
115        Base* b = z      -       num_taylor_per_var; // called y in documentation
116
117        size_t m = (q-1) * r + 1;
118        for(size_t ell = 0; ell < r; ell++)
119        {       b[m+ell] = Base(2) * x[m+ell] * x[0];
120                z[m+ell] = Base(q) * x[m+ell];
121                for(size_t k = 1; k < q; k++)
122                {       b[m+ell] += x[(k-1)*r+1+ell] * x[(q-k-1)*r+1+ell];
123                        z[m+ell] -= Base(k) * z[(k-1)*r+1+ell] * b[(q-k-1)*r+1+ell];
124                }
125                z[m+ell] /= ( Base(q) * b[0] );
126        }
127}
128
129/*!
130Zero order forward mode Taylor coefficient for result of op = AtanOp.
131
132The C++ source code corresponding to this operation is
133\verbatim
134        z = atan(x)
135\endverbatim
136The auxillary result is
137\verbatim
138        y = 1 + x * x
139\endverbatim
140The value of y is computed along with the value of z.
141
142\copydetails forward_unary2_op_0
143*/
144template <class Base>
145inline void forward_atan_op_0(
146        size_t i_z         ,
147        size_t i_x         ,
148        size_t cap_order   , 
149        Base*  taylor      )
150{
151        // check assumptions
152        CPPAD_ASSERT_UNKNOWN( NumArg(AtanOp) == 1 );
153        CPPAD_ASSERT_UNKNOWN( NumRes(AtanOp) == 2 );
154        CPPAD_ASSERT_UNKNOWN( 0 < cap_order );
155
156        // Taylor coefficients corresponding to argument and result
157        Base* x = taylor + i_x * cap_order;
158        Base* z = taylor + i_z * cap_order;
159        Base* b = z      -       cap_order; // called y in documentation
160
161        z[0] = atan( x[0] );
162        b[0] = Base(1) + x[0] * x[0];
163}
164/*!
165Reverse mode partial derivatives for result of op = AtanOp.
166
167The C++ source code corresponding to this operation is
168\verbatim
169        z = atan(x)
170\endverbatim
171The auxillary result is
172\verbatim
173        y = 1 + x * x
174\endverbatim
175The value of y is computed along with the value of z.
176
177\copydetails reverse_unary2_op
178*/
179
180template <class Base>
181inline void reverse_atan_op(
182        size_t      d            ,
183        size_t      i_z          ,
184        size_t      i_x          ,
185        size_t      cap_order    , 
186        const Base* taylor       ,
187        size_t      nc_partial   ,
188        Base*       partial      )
189{
190        // check assumptions
191        CPPAD_ASSERT_UNKNOWN( NumArg(AtanOp) == 1 );
192        CPPAD_ASSERT_UNKNOWN( NumRes(AtanOp) == 2 );
193        CPPAD_ASSERT_UNKNOWN( d < cap_order );
194        CPPAD_ASSERT_UNKNOWN( d < nc_partial );
195
196        // Taylor coefficients and partials corresponding to argument
197        const Base* x  = taylor  + i_x * cap_order;
198        Base* px       = partial + i_x * nc_partial;
199
200        // Taylor coefficients and partials corresponding to first result
201        const Base* z  = taylor  + i_z * cap_order;
202        Base* pz       = partial + i_z * nc_partial;
203
204        // Taylor coefficients and partials corresponding to auxillary result
205        const Base* b  = z  - cap_order; // called y in documentation
206        Base* pb       = pz - nc_partial;
207
208        // number of indices to access
209        size_t j = d;
210        size_t k;
211        while(j)
212        {       // scale partials w.r.t z[j] and b[j]
213                pz[j] /= b[0];
214                pb[j] *= Base(2);
215
216                pb[0] -= pz[j] * z[j]; 
217                px[j] += pz[j] + pb[j] * x[0];
218                px[0] += pb[j] * x[j];
219
220                // more scaling of partials w.r.t z[j]
221                pz[j] /= Base(j);
222
223                for(k = 1; k < j; k++)
224                {       pb[j-k] -= pz[j] * Base(k) * z[k];
225                        pz[k]   -= pz[j] * Base(k) * b[j-k];
226                        px[k]   += pb[j] * x[j-k];
227                }
228                --j;
229        }
230        px[0] += pz[0] / b[0] + pb[0] * Base(2) * x[0];
231}
232
233} // END_CPPAD_NAMESPACE
234# endif
Note: See TracBrowser for help on using the repository browser.