source: trunk/cppad/local/div_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: 15.4 KB
Line 
1/* $Id: div_op.hpp 3301 2014-05-24 05:20:21Z bradbell $ */
2# ifndef CPPAD_DIV_OP_INCLUDED
3# define CPPAD_DIV_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
16namespace CppAD { // BEGIN_CPPAD_NAMESPACE
17/*!
18\file div_op.hpp
19Forward and reverse mode calculations for z = x / y.
20*/
21
22// --------------------------- Divvv -----------------------------------------
23/*!
24Compute forward mode Taylor coefficients for result of op = DivvvOp.
25
26The C++ source code corresponding to this operation is
27\verbatim
28        z = x / y
29\endverbatim
30In the documentation below,
31this operations is for the case where both x and y are variables
32and the argument \a parameter is not used.
33
34\copydetails forward_binary_op
35*/
36
37template <class Base>
38inline void forward_divvv_op(
39        size_t        p           , 
40        size_t        q           , 
41        size_t        i_z         ,
42        const addr_t* arg         ,
43        const Base*   parameter   ,
44        size_t        cap_order   ,
45        Base*         taylor      )
46{
47        // check assumptions
48        CPPAD_ASSERT_UNKNOWN( NumArg(DivvvOp) == 2 );
49        CPPAD_ASSERT_UNKNOWN( NumRes(DivvvOp) == 1 );
50        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < i_z );
51        CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < i_z );
52        CPPAD_ASSERT_UNKNOWN( q < cap_order );
53        CPPAD_ASSERT_UNKNOWN( p <= q );
54
55        // Taylor coefficients corresponding to arguments and result
56        Base* x = taylor + arg[0] * cap_order;
57        Base* y = taylor + arg[1] * cap_order;
58        Base* z = taylor + i_z    * cap_order;
59
60
61        // Using CondExp, it can make sense to divide by zero,
62        // so do not make it an error.
63        size_t k;
64        for(size_t d = p; d <= q; d++)
65        {       z[d] = x[d];
66                for(k = 1; k <= d; k++)
67                        z[d] -= z[d-k] * y[k];
68                z[d] /= y[0];
69        }
70}
71/*!
72Multiple directions forward mode Taylor coefficients for op = DivvvOp.
73
74The C++ source code corresponding to this operation is
75\verbatim
76        z = x / y
77\endverbatim
78In the documentation below,
79this operations is for the case where both x and y are variables
80and the argument \a parameter is not used.
81
82\copydetails forward_binary_op_dir
83*/
84
85template <class Base>
86inline void forward_divvv_op_dir(
87        size_t        q           , 
88        size_t        r           , 
89        size_t        i_z         ,
90        const addr_t* arg         ,
91        const Base*   parameter   ,
92        size_t        cap_order   ,
93        Base*         taylor      )
94{
95        // check assumptions
96        CPPAD_ASSERT_UNKNOWN( NumArg(DivvvOp) == 2 );
97        CPPAD_ASSERT_UNKNOWN( NumRes(DivvvOp) == 1 );
98        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < i_z );
99        CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < i_z );
100        CPPAD_ASSERT_UNKNOWN( 0 < q );
101        CPPAD_ASSERT_UNKNOWN( q < cap_order );
102
103        // Taylor coefficients corresponding to arguments and result
104        size_t num_taylor_per_var = (cap_order-1) * r + 1;
105        Base* x = taylor + arg[0] * num_taylor_per_var;
106        Base* y = taylor + arg[1] * num_taylor_per_var;
107        Base* z = taylor + i_z    * num_taylor_per_var;
108
109
110        // Using CondExp, it can make sense to divide by zero,
111        // so do not make it an error.
112        size_t m = (q-1) * r + 1;
113        for(size_t ell = 0; ell < r; ell++)
114        {       z[m+ell] = x[m+ell] - z[0] * y[m+ell];
115                for(size_t k = 1; k < q; k++)           
116                        z[m+ell] -= z[(q-k-1)*r+1+ell] * y[(k-1)*r+1+ell];
117                z[m+ell] /= y[0];
118        }
119}
120
121
122/*!
123Compute zero order forward mode Taylor coefficients for result of op = DivvvOp.
124
125The C++ source code corresponding to this operation is
126\verbatim
127        z = x / y
128\endverbatim
129In the documentation below,
130this operations is for the case where both x and y are variables
131and the argument \a parameter is not used.
132
133\copydetails forward_binary_op_0
134*/
135
136template <class Base>
137inline void forward_divvv_op_0(
138        size_t        i_z         ,
139        const addr_t* arg         ,
140        const Base*   parameter   ,
141        size_t        cap_order   ,
142        Base*         taylor      )
143{
144        // check assumptions
145        CPPAD_ASSERT_UNKNOWN( NumArg(DivvvOp) == 2 );
146        CPPAD_ASSERT_UNKNOWN( NumRes(DivvvOp) == 1 );
147        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < i_z );
148        CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < i_z );
149
150        // Taylor coefficients corresponding to arguments and result
151        Base* x = taylor + arg[0] * cap_order;
152        Base* y = taylor + arg[1] * cap_order;
153        Base* z = taylor + i_z    * cap_order;
154
155        z[0] = x[0] / y[0];
156}
157
158/*!
159Compute reverse mode partial derivatives for result of op = DivvvOp.
160
161The C++ source code corresponding to this operation is
162\verbatim
163        z = x / y
164\endverbatim
165In the documentation below,
166this operations is for the case where both x and y are variables
167and the argument \a parameter is not used.
168
169\copydetails reverse_binary_op
170*/
171
172template <class Base>
173inline void reverse_divvv_op(
174        size_t        d           , 
175        size_t        i_z         ,
176        const addr_t* arg         ,
177        const Base*   parameter   ,
178        size_t        cap_order   ,
179        const Base*   taylor      ,
180        size_t        nc_partial  ,
181        Base*         partial     )
182{
183        // check assumptions
184        CPPAD_ASSERT_UNKNOWN( NumArg(DivvvOp) == 2 );
185        CPPAD_ASSERT_UNKNOWN( NumRes(DivvvOp) == 1 );
186        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < i_z );
187        CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < i_z );
188        CPPAD_ASSERT_UNKNOWN( d < cap_order );
189        CPPAD_ASSERT_UNKNOWN( d < nc_partial );
190
191        // Arguments
192        const Base* y  = taylor + arg[1] * cap_order;
193        const Base* z  = taylor + i_z    * cap_order;
194
195        // Partial derivatives corresponding to arguments and result
196        Base* px = partial + arg[0] * nc_partial;
197        Base* py = partial + arg[1] * nc_partial;
198        Base* pz = partial + i_z    * nc_partial;
199
200        // Using CondExp, it can make sense to divide by zero
201        // so do not make it an error.
202
203        size_t k;
204        // number of indices to access
205        size_t j = d + 1;
206        while(j)
207        {       --j;
208                // scale partial w.r.t. z[j]
209                pz[j] /= y[0];
210
211                px[j] += pz[j];
212                for(k = 1; k <= j; k++)
213                {       pz[j-k] -= pz[j] * y[k];
214                        py[k]   -= pz[j] * z[j-k];
215                }       
216                py[0] -= pz[j] * z[j];
217        }
218}
219
220// --------------------------- Divpv -----------------------------------------
221/*!
222Compute forward mode Taylor coefficients for result of op = DivpvOp.
223
224The C++ source code corresponding to this operation is
225\verbatim
226        z = x / y
227\endverbatim
228In the documentation below,
229this operations is for the case where x is a parameter and y is a variable.
230
231\copydetails forward_binary_op
232*/
233
234template <class Base>
235inline void forward_divpv_op(
236        size_t        p           , 
237        size_t        q           , 
238        size_t        i_z         ,
239        const addr_t* arg         ,
240        const Base*   parameter   ,
241        size_t        cap_order   ,
242        Base*         taylor      )
243{
244        // check assumptions
245        CPPAD_ASSERT_UNKNOWN( NumArg(DivpvOp) == 2 );
246        CPPAD_ASSERT_UNKNOWN( NumRes(DivpvOp) == 1 );
247        CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < i_z );
248        CPPAD_ASSERT_UNKNOWN( q < cap_order );
249        CPPAD_ASSERT_UNKNOWN( p <= q );
250
251        // Taylor coefficients corresponding to arguments and result
252        Base* y = taylor + arg[1] * cap_order;
253        Base* z = taylor + i_z    * cap_order;
254
255        // Paraemter value
256        Base x = parameter[ arg[0] ];
257
258        // Using CondExp, it can make sense to divide by zero,
259        // so do not make it an error.
260        size_t k;
261        if( p == 0 )
262        {       z[0] = x / y[0];
263                p++;
264        }
265        for(size_t d = p; d <= q; d++)
266        {       z[d] = Base(0);
267                for(k = 1; k <= d; k++)
268                        z[d] -= z[d-k] * y[k];
269                z[d] /= y[0];
270        }
271}
272/*!
273Multiple directions forward mode Taylor coefficients for op = DivpvOp.
274
275The C++ source code corresponding to this operation is
276\verbatim
277        z = x / y
278\endverbatim
279In the documentation below,
280this operations is for the case where x is a parameter and y is a variable.
281
282\copydetails forward_binary_op_dir
283*/
284
285template <class Base>
286inline void forward_divpv_op_dir(
287        size_t        q           , 
288        size_t        r           , 
289        size_t        i_z         ,
290        const addr_t* arg         ,
291        const Base*   parameter   ,
292        size_t        cap_order   ,
293        Base*         taylor      )
294{
295        // check assumptions
296        CPPAD_ASSERT_UNKNOWN( NumArg(DivpvOp) == 2 );
297        CPPAD_ASSERT_UNKNOWN( NumRes(DivpvOp) == 1 );
298        CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < i_z );
299        CPPAD_ASSERT_UNKNOWN( 0 < q );
300        CPPAD_ASSERT_UNKNOWN( q < cap_order );
301
302        // Taylor coefficients corresponding to arguments and result
303        size_t num_taylor_per_var = (cap_order-1) * r + 1;
304        Base* y = taylor + arg[1] * num_taylor_per_var;
305        Base* z = taylor + i_z    * num_taylor_per_var;
306
307        // Using CondExp, it can make sense to divide by zero,
308        // so do not make it an error.
309        size_t m = (q-1) * r + 1;
310        for(size_t ell = 0; ell < r; ell++)
311        {       z[m+ell] = - z[0] * y[m+ell];
312                for(size_t k = 1; k < q; k++)
313                        z[m+ell] -= z[(q-k-1)*r+1+ell] * y[(k-1)*r+1+ell];
314                z[m+ell] /= y[0];
315        }
316}
317
318/*!
319Compute zero order forward mode Taylor coefficient for result of op = DivpvOp.
320
321The C++ source code corresponding to this operation is
322\verbatim
323        z = x / y
324\endverbatim
325In the documentation below,
326this operations is for the case where x is a parameter and y is a variable.
327
328\copydetails forward_binary_op_0
329*/
330
331template <class Base>
332inline void forward_divpv_op_0(
333        size_t        i_z         ,
334        const addr_t* arg         ,
335        const Base*   parameter   ,
336        size_t        cap_order   ,
337        Base*         taylor      )
338{
339        // check assumptions
340        CPPAD_ASSERT_UNKNOWN( NumArg(DivpvOp) == 2 );
341        CPPAD_ASSERT_UNKNOWN( NumRes(DivpvOp) == 1 );
342        CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < i_z );
343
344        // Paraemter value
345        Base x = parameter[ arg[0] ];
346
347        // Taylor coefficients corresponding to arguments and result
348        Base* y = taylor + arg[1] * cap_order;
349        Base* z = taylor + i_z    * cap_order;
350
351        z[0] = x / y[0];
352}
353
354/*!
355Compute reverse mode partial derivative for result of op = DivpvOp.
356
357The C++ source code corresponding to this operation is
358\verbatim
359        z = x / y
360\endverbatim
361In the documentation below,
362this operations is for the case where x is a parameter and y is a variable.
363
364\copydetails reverse_binary_op
365*/
366
367template <class Base>
368inline void reverse_divpv_op(
369        size_t        d           , 
370        size_t        i_z         ,
371        const addr_t* arg         ,
372        const Base*   parameter   ,
373        size_t        cap_order   ,
374        const Base*   taylor      ,
375        size_t        nc_partial  ,
376        Base*         partial     )
377{
378        // check assumptions
379        CPPAD_ASSERT_UNKNOWN( NumArg(DivvvOp) == 2 );
380        CPPAD_ASSERT_UNKNOWN( NumRes(DivvvOp) == 1 );
381        CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < i_z );
382        CPPAD_ASSERT_UNKNOWN( d < cap_order );
383        CPPAD_ASSERT_UNKNOWN( d < nc_partial );
384
385        // Arguments
386        const Base* y = taylor + arg[1] * cap_order;
387        const Base* z = taylor + i_z    * cap_order;
388
389        // Partial derivatives corresponding to arguments and result
390        Base* py = partial + arg[1] * nc_partial;
391        Base* pz = partial + i_z    * nc_partial;
392
393        // Using CondExp, it can make sense to divide by zero so do not
394        // make it an error.
395
396        size_t k;
397        // number of indices to access
398        size_t j = d + 1;
399        while(j)
400        {       --j;
401                // scale partial w.r.t z[j]
402                pz[j] /= y[0];
403
404                for(k = 1; k <= j; k++)
405                {       pz[j-k] -= pz[j] * y[k];
406                        py[k]   -= pz[j] * z[j-k];
407                }       
408                py[0] -= pz[j] * z[j];
409        }
410}
411
412
413// --------------------------- Divvp -----------------------------------------
414/*!
415Compute forward mode Taylor coefficients for result of op = DivvvOp.
416
417The C++ source code corresponding to this operation is
418\verbatim
419        z = x / y
420\endverbatim
421In the documentation below,
422this operations is for the case where x is a variable and y is a parameter.
423
424\copydetails forward_binary_op
425*/
426
427template <class Base>
428inline void forward_divvp_op(
429        size_t        p           , 
430        size_t        q           , 
431        size_t        i_z         ,
432        const addr_t* arg         ,
433        const Base*   parameter   ,
434        size_t        cap_order   ,
435        Base*         taylor      )
436{
437        // check assumptions
438        CPPAD_ASSERT_UNKNOWN( NumArg(DivvpOp) == 2 );
439        CPPAD_ASSERT_UNKNOWN( NumRes(DivvpOp) == 1 );
440        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < i_z );
441        CPPAD_ASSERT_UNKNOWN( q < cap_order );
442        CPPAD_ASSERT_UNKNOWN( p <= q );
443
444        // Taylor coefficients corresponding to arguments and result
445        Base* x = taylor + arg[0] * cap_order;
446        Base* z = taylor + i_z    * cap_order;
447
448        // Parameter value
449        Base y = parameter[ arg[1] ];
450
451        // Using CondExp and multiple levels of AD, it can make sense
452        // to divide by zero so do not make it an error.
453        for(size_t d = p; d <= q; d++)
454                z[d] = x[d] / y;
455}
456/*!
457Multiple direction forward mode Taylor coefficients for op = DivvvOp.
458
459The C++ source code corresponding to this operation is
460\verbatim
461        z = x / y
462\endverbatim
463In the documentation below,
464this operations is for the case where x is a variable and y is a parameter.
465
466\copydetails forward_binary_op_dir
467*/
468
469template <class Base>
470inline void forward_divvp_op_dir(
471        size_t        q           , 
472        size_t        r           , 
473        size_t        i_z         ,
474        const addr_t* arg         ,
475        const Base*   parameter   ,
476        size_t        cap_order   ,
477        Base*         taylor      )
478{
479        // check assumptions
480        CPPAD_ASSERT_UNKNOWN( NumArg(DivvpOp) == 2 );
481        CPPAD_ASSERT_UNKNOWN( NumRes(DivvpOp) == 1 );
482        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < i_z );
483        CPPAD_ASSERT_UNKNOWN( q < cap_order );
484        CPPAD_ASSERT_UNKNOWN( 0 < q  );
485
486        // Taylor coefficients corresponding to arguments and result
487        size_t num_taylor_per_var = (cap_order-1) * r + 1;
488        Base* x = taylor + arg[0] * num_taylor_per_var;
489        Base* z = taylor +    i_z * num_taylor_per_var;
490
491        // Parameter value
492        Base y = parameter[ arg[1] ];
493
494        // Using CondExp and multiple levels of AD, it can make sense
495        // to divide by zero so do not make it an error.
496        size_t m = (q-1)*r + 1;
497        for(size_t ell = 0; ell < r; ell++)
498                z[m + ell] = x[m + ell] / y;
499}
500
501
502/*!
503Compute zero order forward mode Taylor coefficients for result of op = DivvvOp.
504
505The C++ source code corresponding to this operation is
506\verbatim
507        z = x / y
508\endverbatim
509In the documentation below,
510this operations is for the case where x is a variable and y is a parameter.
511
512\copydetails forward_binary_op_0
513*/
514
515template <class Base>
516inline void forward_divvp_op_0(
517        size_t        i_z         ,
518        const addr_t* arg         ,
519        const Base*   parameter   ,
520        size_t        cap_order   ,
521        Base*         taylor      )
522{
523        // check assumptions
524        CPPAD_ASSERT_UNKNOWN( NumArg(DivvpOp) == 2 );
525        CPPAD_ASSERT_UNKNOWN( NumRes(DivvpOp) == 1 );
526        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < i_z );
527
528        // Parameter value
529        Base y = parameter[ arg[1] ];
530
531        // Taylor coefficients corresponding to arguments and result
532        Base* x = taylor + arg[0] * cap_order;
533        Base* z = taylor + i_z    * cap_order;
534
535        z[0] = x[0] / y;
536}
537
538/*!
539Compute reverse mode partial derivative for result of op = DivvpOp.
540
541The C++ source code corresponding to this operation is
542\verbatim
543        z = x / y
544\endverbatim
545In the documentation below,
546this operations is for the case where x is a variable and y is a parameter.
547
548\copydetails reverse_binary_op
549*/
550
551template <class Base>
552inline void reverse_divvp_op(
553        size_t        d           , 
554        size_t        i_z         ,
555        const addr_t* arg         ,
556        const Base*   parameter   ,
557        size_t        cap_order   ,
558        const Base*   taylor      ,
559        size_t        nc_partial  ,
560        Base*         partial     )
561{
562        // check assumptions
563        CPPAD_ASSERT_UNKNOWN( NumArg(DivvpOp) == 2 );
564        CPPAD_ASSERT_UNKNOWN( NumRes(DivvpOp) == 1 );
565        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < i_z );
566        CPPAD_ASSERT_UNKNOWN( d < cap_order );
567        CPPAD_ASSERT_UNKNOWN( d < nc_partial );
568
569        // Argument values
570        Base  y = parameter[ arg[1] ];
571
572        // Partial derivatives corresponding to arguments and result
573        Base* px = partial + arg[0] * nc_partial;
574        Base* pz = partial + i_z    * nc_partial;
575
576        // Using CondExp, it can make sense to divide by zero
577        // so do not make it an error.
578
579        // number of indices to access
580        size_t j = d + 1;
581        while(j)
582        {       --j;
583                px[j] += pz[j] / y;
584        }
585}
586
587} // END_CPPAD_NAMESPACE
588# endif
Note: See TracBrowser for help on using the repository browser.