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