source: trunk/cppad/local/acos_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: 6.5 KB
Line 
1/* $Id: acos_op.hpp 3320 2014-09-11 23:06:21Z bradbell $ */
2# ifndef CPPAD_ACOS_OP_INCLUDED
3# define CPPAD_ACOS_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 acos_op.hpp
20Forward and reverse mode calculations for z = acos(x).
21*/
22
23
24/*!
25Compute forward mode Taylor coefficient for result of op = AcosOp.
26
27The C++ source code corresponding to this operation is
28\verbatim
29        z = acos(x)
30\endverbatim
31The auxillary result is
32\verbatim
33        y = sqrt(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_acos_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(AcosOp) == 1 );
51        CPPAD_ASSERT_UNKNOWN( NumRes(AcosOp) == 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        Base uj;
62        if( p == 0 )
63        {       z[0] = acos( x[0] );
64                uj   = Base(1) - x[0] * x[0];
65                b[0] = sqrt( uj );
66                p++;
67        }
68        for(size_t j = p; j <= q; j++)
69        {       uj = Base(0);
70                for(k = 0; k <= j; k++)
71                        uj -= x[k] * x[j-k];
72                b[j] = Base(0);
73                z[j] = Base(0);
74                for(k = 1; k < j; k++)
75                {       b[j] -= Base(k) * b[k] * b[j-k];
76                        z[j] -= Base(k) * z[k] * b[j-k];
77                }
78                b[j] /= Base(j);
79                z[j] /= Base(j);
80                //
81                b[j] += uj / Base(2);
82                z[j] -= x[j];
83                //
84                b[j] /= b[0];
85                z[j] /= b[0];
86        }
87}
88/*!
89Multiple directions forward mode Taylor coefficient for op = AcosOp.
90
91The C++ source code corresponding to this operation is
92\verbatim
93        z = acos(x)
94\endverbatim
95The auxillary result is
96\verbatim
97        y = sqrt(1 - x * x)
98\endverbatim
99The value of y, and its derivatives, are computed along with the value
100and derivatives of z.
101
102\copydetails forward_unary2_op_dir
103*/
104template <class Base>
105inline void forward_acos_op_dir(
106        size_t q           ,
107        size_t r           ,
108        size_t i_z         ,
109        size_t i_x         ,
110        size_t cap_order   , 
111        Base*  taylor      )
112{       
113        // check assumptions
114        CPPAD_ASSERT_UNKNOWN( NumArg(AcosOp) == 1 );
115        CPPAD_ASSERT_UNKNOWN( NumRes(AcosOp) == 2 );
116        CPPAD_ASSERT_UNKNOWN( 0 < q );
117        CPPAD_ASSERT_UNKNOWN( q < cap_order );
118
119        // Taylor coefficients corresponding to argument and result
120        size_t num_taylor_per_var = (cap_order-1) * r + 1;
121        Base* x = taylor + i_x * num_taylor_per_var;
122        Base* z = taylor + i_z * num_taylor_per_var;
123        Base* b = z - num_taylor_per_var;  // called y in documentation
124
125        size_t k, ell;
126        size_t m = (q-1) * r + 1;
127        for(ell = 0; ell < r; ell ++)
128        {       Base uq = - 2.0 * x[m + ell] * x[0];
129                for(k = 1; k < q; k++)
130                        uq -= x[(k-1)*r+1+ell] * x[(q-k-1)*r+1+ell]; 
131                b[m+ell] = Base(0);
132                z[m+ell] = Base(0);
133                for(k = 1; k < q; k++)
134                {       b[m+ell] += Base(k) * b[(k-1)*r+1+ell] * b[(q-k-1)*r+1+ell]; 
135                        z[m+ell] += Base(k) * z[(k-1)*r+1+ell] * b[(q-k-1)*r+1+ell]; 
136                }
137                b[m+ell] =  ( uq / Base(2) - b[m+ell] / Base(q) ) / b[0];
138                z[m+ell] = -( x[m+ell]     + z[m+ell] / Base(q) ) / b[0];
139        }
140}
141
142/*!
143Compute zero order forward mode Taylor coefficient for result of op = AcosOp.
144
145The C++ source code corresponding to this operation is
146\verbatim
147        z = acos(x)
148\endverbatim
149The auxillary result is
150\verbatim
151        y = sqrt( 1 - x * x )
152\endverbatim
153The value of y is computed along with the value of z.
154
155\copydetails forward_unary2_op_0
156*/
157template <class Base>
158inline void forward_acos_op_0(
159        size_t i_z         ,
160        size_t i_x         ,
161        size_t cap_order   , 
162        Base*  taylor      )
163{
164        // check assumptions
165        CPPAD_ASSERT_UNKNOWN( NumArg(AcosOp) == 1 );
166        CPPAD_ASSERT_UNKNOWN( NumRes(AcosOp) == 2 );
167        CPPAD_ASSERT_UNKNOWN( 0 < cap_order );
168
169        // Taylor coefficients corresponding to argument and result
170        Base* x = taylor + i_x * cap_order;
171        Base* z = taylor + i_z * cap_order;
172        Base* b = z      -       cap_order; // called y in documentation
173
174        z[0] = acos( x[0] );
175        b[0] = sqrt( Base(1) - x[0] * x[0] );
176}
177/*!
178Compute reverse mode partial derivatives for result of op = AcosOp.
179
180The C++ source code corresponding to this operation is
181\verbatim
182        z = acos(x)
183\endverbatim
184The auxillary result is
185\verbatim
186        y = sqrt( 1 - x * x )
187\endverbatim
188The value of y is computed along with the value of z.
189
190\copydetails reverse_unary2_op
191*/
192
193template <class Base>
194inline void reverse_acos_op(
195        size_t      d            ,
196        size_t      i_z          ,
197        size_t      i_x          ,
198        size_t      cap_order    , 
199        const Base* taylor       ,
200        size_t      nc_partial   ,
201        Base*       partial      )
202{
203        // check assumptions
204        CPPAD_ASSERT_UNKNOWN( NumArg(AcosOp) == 1 );
205        CPPAD_ASSERT_UNKNOWN( NumRes(AcosOp) == 2 );
206        CPPAD_ASSERT_UNKNOWN( d < cap_order );
207        CPPAD_ASSERT_UNKNOWN( d < nc_partial );
208
209        // Taylor coefficients and partials corresponding to argument
210        const Base* x  = taylor  + i_x * cap_order;
211        Base* px       = partial + i_x * nc_partial;
212
213        // Taylor coefficients and partials corresponding to first result
214        const Base* z  = taylor  + i_z * cap_order;
215        Base* pz       = partial + i_z * nc_partial;
216
217        // Taylor coefficients and partials corresponding to auxillary result
218        const Base* b  = z  - cap_order; // called y in documentation
219        Base* pb       = pz - nc_partial;
220
221        // number of indices to access
222        size_t j = d;
223        size_t k;
224        while(j)
225        {
226                // scale partials w.r.t b[j] by 1 / b[0]
227                pb[j] /= b[0];
228
229                // scale partials w.r.t z[j] by 1 / b[0]
230                pz[j] /= b[0];
231
232                // update partials w.r.t b^0
233                pb[0] -= pz[j] * z[j] + pb[j] * b[j]; 
234
235                // update partial w.r.t. x^0
236                px[0] -= pb[j] * x[j];
237
238                // update partial w.r.t. x^j
239                px[j] -= pz[j] + pb[j] * x[0];
240
241                // further scale partial w.r.t. z[j] by 1 / j
242                pz[j] /= Base(j);
243
244                for(k = 1; k < j; k++)
245                {       // update partials w.r.t b^(j-k)
246                        pb[j-k] -= Base(k) * pz[j] * z[k] + pb[j] * b[k];
247
248                        // update partials w.r.t. x^k
249                        px[k]   -= pb[j] * x[j-k];
250
251                        // update partials w.r.t. z^k
252                        pz[k]   -= pz[j] * Base(k) * b[j-k];
253                }
254                --j;
255        }
256
257        // j == 0 case
258        px[0] -= ( pz[0] + pb[0] * x[0]) / b[0];
259}
260
261} // END_CPPAD_NAMESPACE
262# endif
Note: See TracBrowser for help on using the repository browser.