source: trunk/cppad/local/store_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: 16.3 KB
Line 
1/* $Id: store_op.hpp 3301 2014-05-24 05:20:21Z bradbell $ */
2# ifndef CPPAD_STORE_OP_INCLUDED
3# define CPPAD_STORE_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 store_op.hpp
19Changing the current value of a VecAD element.
20*/
21/*!
22Shared documentation for zero order forward implementation of
23op = StppOp, StpvOp, StvpOp, or StvvOp (not called).
24
25The C++ source code corresponding to this operation is
26\verbatim
27        v[x] = y
28\endverbatim
29where v is a VecAD<Base> vector, x is an AD<Base> object,
30and y is AD<Base> or Base objects.
31We define the index corresponding to v[x] by
32\verbatim
33        i_v_x = index_by_ind[ arg[0] + i_vec ]
34\endverbatim
35where i_vec is defined under the heading arg[1] below:
36
37\tparam Base
38base type for the operator; i.e., this operation was recorded
39using AD<Base> and computations by this routine are done using type Base.
40
41\param i_z
42is the index corresponding to the previous variable on the tape
43(only used for error checking).
44
45\param arg
46\n
47arg[0]
48\n
49is the offset of this VecAD vector relative to the beginning
50of the isvar_by_ind and index_by_ind arrays.
51\n
52\n
53arg[1]
54\n
55If this is a StppOp or StpvOp operation (if x is a parameter),
56i_vec is defined by
57\verbatim
58        i_vec = arg[1]
59\endverbatim
60If this is a StvpOp or StvvOp operation (if x is a variable),
61i_vec is defined by
62\verbatim
63        i_vec = floor( taylor[ arg[1] * cap_order + 0 ] )
64\endverbatim
65where floor(c) is the greatest integer less that or equal c.
66\n
67\n
68arg[2]
69\n
70index corresponding to the third operand for this operator;
71i.e. the index corresponding to y.
72
73\param num_par
74is the total number of parameters on the tape
75(only used for error checking).
76
77\param cap_order
78number of columns in the matrix containing the Taylor coefficients.
79
80\param taylor
81In StvpOp and StvvOp cases, <code><taylor[ arg[1] * cap_order + 0 ]</code>
82is used to compute the index in the definition of i_vec above.
83
84\param isvar_by_ind
85If y is a varable (StpvOp and StvvOp cases),
86<code>isvar_by_ind[ arg[0] + i_vec ] </code> is set to true.
87Otherwise y is a paraemter (StppOp and StvpOp cases) and
88<code>isvar_by_ind[ arg[0] + i_vec ] </code> is set to false.
89
90\param index_by_ind
91<code>index_by_ind[ arg[0] - 1 ]</code>
92is the number of elements in the user vector containing this element.
93The value <code>index_by_ind[ arg[0] + i_vec]</code>
94is set equal to arg[2].
95
96\par Check User Errors
97\li Check that the index is with in range; i.e.
98<code>i_vec < index_by_ind[ arg[0] - 1 ]</code>
99Note that, if x is a parameter,
100the corresponding vector index and it does not change.
101In this case, the error above should be detected during tape recording.
102
103\par Checked Assertions
104\li NumArg(op) == 3
105\li NumRes(op) == 0
106\li 0 <  arg[0]
107\li if y is a parameter, arg[2] < num_par
108\li if x is a variable, arg[1] <= i_z
109\li if y is a variable, arg[2] <= i_z
110*/
111template <class Base>
112inline void forward_store_op_0(
113        size_t         i_z         ,
114        const addr_t*  arg         , 
115        size_t         num_par     ,
116        size_t         cap_order   ,
117        Base*          taylor      ,
118        bool*          isvar_by_ind   ,
119        size_t*        index_by_ind   )
120{
121        // This routine is only for documentaiton, it should not be used
122        CPPAD_ASSERT_UNKNOWN( false );
123}
124/*!
125Shared documnetation for sparsity operations corresponding to
126op = StpvOp or StvvOp (not called).
127
128<!-- define sparse_store_op -->
129The C++ source code corresponding to this operation is
130\verbatim
131        v[x] = y
132\endverbatim
133where v is a VecAD<Base> vector, x is an AD<Base> object,
134and y is AD<Base> or Base objects.
135We define the index corresponding to v[x] by
136\verbatim
137        i_v_x = combined[ arg[0] + i_vec ]
138\endverbatim
139where i_vec is defined under the heading \a arg[1] below:
140
141\tparam Vector_set
142is the type used for vectors of sets. It can be either
143\c sparse_pack, \c sparse_set, or \c sparse_list.
144
145\param op
146is the code corresponding to this operator; i.e., StpvOp or StvvOp
147(only used for error checking).
148
149\param arg
150\n
151\a arg[0]
152is the offset corresponding to this VecAD vector in the combined array.
153\n
154\n
155\a arg[2]
156\n
157The set with index \a arg[2] in \a var_sparsity
158is the sparsity pattern corresponding to y.
159(Note that \a arg[2] > 0 because y is a variable.)
160
161\param num_combined
162is the total number of elements in the VecAD address array.
163
164\param combined
165\a combined [ arg[0] - 1 ]
166is the index of the set in \a vecad_sparsity corresponding
167to the sparsity pattern for the vector v.
168We use the notation i_v below which is defined by
169\verbatim
170        i_v = combined[ \a arg[0] - 1 ]
171\endverbatim
172
173\param var_sparsity
174The set  with index \a arg[2] in \a var_sparsity
175is the sparsity pattern for y.
176This is an input for forward mode operations.
177For reverse mode operations:
178The sparsity pattern for v is added to the spartisy pattern for y.
179
180\param vecad_sparsity
181The set with index \a i_v in \a vecad_sparsity
182is the sparsity pattern for v.
183This is an input for reverse mode operations.
184For forward mode operations, the sparsity pattern for y is added
185to the sparsity pattern for the vector v.
186
187\par Checked Assertions
188\li NumArg(op) == 3
189\li NumRes(op) == 0
190\li 0 <  \a arg[0]
191\li \a arg[0] < \a num_combined
192\li \a arg[2] < \a var_sparsity.n_set()
193\li i_v       < \a vecad_sparsity.n_set()
194<!-- end sparse_store_op -->
195*/
196template <class Vector_set>
197inline void sparse_store_op(
198        OpCode         op             ,
199        const addr_t*  arg            , 
200        size_t         num_combined   ,
201        const size_t*  combined       ,
202        Vector_set&    var_sparsity   ,
203        Vector_set&    vecad_sparsity )
204{
205        // This routine is only for documentaiton, it should not be used
206        CPPAD_ASSERT_UNKNOWN( false );
207}
208
209
210/*!
211Zero order forward mode implementation of op = StppOp.
212
213\copydetails forward_store_op_0
214*/
215template <class Base>
216inline void forward_store_pp_op_0(
217        size_t         i_z         ,
218        const addr_t*  arg         , 
219        size_t         num_par     ,
220        size_t         cap_order   ,
221        Base*          taylor      ,
222        bool*          isvar_by_ind   ,
223        size_t*        index_by_ind   )
224{       size_t i_vec = arg[1];
225
226        // Because the index is a parameter, this indexing error should be
227        // caught and reported to the user when the tape is recording.
228        CPPAD_ASSERT_UNKNOWN( i_vec < index_by_ind[ arg[0] - 1 ] );
229
230        CPPAD_ASSERT_UNKNOWN( NumArg(StppOp) == 3 );
231        CPPAD_ASSERT_UNKNOWN( NumRes(StppOp) == 0 );
232        CPPAD_ASSERT_UNKNOWN( 0 < arg[0] );
233        CPPAD_ASSERT_UNKNOWN( size_t(arg[2]) < num_par );
234
235        isvar_by_ind[ arg[0] + i_vec ]  = false;
236        index_by_ind[ arg[0] + i_vec ]  = arg[2];
237}
238
239/*!
240Zero order forward mode implementation of op = StpvOp.
241
242\copydetails forward_store_op_0
243*/
244template <class Base>
245inline void forward_store_pv_op_0(
246        size_t         i_z         ,
247        const addr_t*  arg         , 
248        size_t         num_par     ,
249        size_t         cap_order   ,
250        Base*          taylor      ,
251        bool*          isvar_by_ind   ,
252        size_t*        index_by_ind   )
253{       size_t i_vec = arg[1];
254
255        // Because the index is a parameter, this indexing error should be
256        // caught and reported to the user when the tape is recording.
257        CPPAD_ASSERT_UNKNOWN( i_vec < index_by_ind[ arg[0] - 1 ] );
258
259        CPPAD_ASSERT_UNKNOWN( NumArg(StpvOp) == 3 );
260        CPPAD_ASSERT_UNKNOWN( NumRes(StpvOp) == 0 );
261        CPPAD_ASSERT_UNKNOWN( 0 < arg[0] );
262        CPPAD_ASSERT_UNKNOWN( size_t(arg[2]) <= i_z );
263
264        isvar_by_ind[ arg[0] + i_vec ]  = true;
265        index_by_ind[ arg[0] + i_vec ]  = arg[2];
266}
267
268/*!
269Zero order forward mode implementation of op = StvpOp.
270
271\copydetails forward_store_op_0
272*/
273template <class Base>
274inline void forward_store_vp_op_0(
275        size_t         i_z         ,
276        const addr_t*  arg         , 
277        size_t         num_par     ,
278        size_t         cap_order   ,
279        Base*          taylor      ,
280        bool*          isvar_by_ind   ,
281        size_t*        index_by_ind   )
282{       
283        CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) <= i_z );
284        size_t i_vec = Integer( taylor[ arg[1] * cap_order + 0 ] );
285        CPPAD_ASSERT_KNOWN( 
286                i_vec < index_by_ind[ arg[0] - 1 ] ,
287                "VecAD: index during zero order forward sweep is out of range"
288        );
289
290        CPPAD_ASSERT_UNKNOWN( NumArg(StvpOp) == 3 );
291        CPPAD_ASSERT_UNKNOWN( NumRes(StvpOp) == 0 );
292        CPPAD_ASSERT_UNKNOWN( 0 < arg[0] );
293        CPPAD_ASSERT_UNKNOWN( size_t(arg[2]) < num_par );
294
295        isvar_by_ind[ arg[0] + i_vec ]  = false;
296        index_by_ind[ arg[0] + i_vec ]  = arg[2];
297}
298
299/*!
300Zero order forward mode implementation of op = StvvOp.
301
302\copydetails forward_store_op_0
303*/
304template <class Base>
305inline void forward_store_vv_op_0(
306        size_t         i_z         ,
307        const addr_t*  arg         , 
308        size_t         num_par     ,
309        size_t         cap_order   ,
310        Base*          taylor      ,
311        bool*          isvar_by_ind   ,
312        size_t*        index_by_ind   )
313{       
314        CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) <= i_z );
315        size_t i_vec = Integer( taylor[ arg[1] * cap_order + 0 ] );
316        CPPAD_ASSERT_KNOWN( 
317                i_vec < index_by_ind[ arg[0] - 1 ] ,
318                "VecAD: index during zero order forward sweep is out of range"
319        );
320
321        CPPAD_ASSERT_UNKNOWN( NumArg(StvpOp) == 3 );
322        CPPAD_ASSERT_UNKNOWN( NumRes(StvpOp) == 0 );
323        CPPAD_ASSERT_UNKNOWN( 0 < arg[0] );
324        CPPAD_ASSERT_UNKNOWN( size_t(arg[2]) <= i_z );
325
326        isvar_by_ind[ arg[0] + i_vec ]  = true;
327        index_by_ind[ arg[0] + i_vec ]  = arg[2];
328}
329
330/*!
331Forward mode sparsity operations for StpvOp and StvvOp
332
333\copydetails sparse_store_op
334*/
335template <class Vector_set>
336inline void forward_sparse_store_op(
337        OpCode              op             ,
338        const addr_t*       arg            , 
339        size_t              num_combined   ,
340        const size_t*       combined       ,
341        Vector_set&         var_sparsity   ,
342        Vector_set&         vecad_sparsity )
343{
344        CPPAD_ASSERT_UNKNOWN( NumArg(op) == 3 );
345        CPPAD_ASSERT_UNKNOWN( NumRes(op) == 0 );
346        CPPAD_ASSERT_UNKNOWN( 0 < arg[0] );
347        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_combined );
348        size_t i_v = combined[ arg[0] - 1 ];
349        CPPAD_ASSERT_UNKNOWN( i_v < vecad_sparsity.n_set() );
350        CPPAD_ASSERT_UNKNOWN( size_t(arg[2]) < var_sparsity.n_set() );
351
352        vecad_sparsity.binary_union(i_v, i_v, arg[2], var_sparsity);
353
354        return;
355}
356
357/*!
358Reverse mode sparsity operations for StpvOp and StvvOp
359
360This routine is given the sparsity patterns for
361G(v[x], y , w , u ... ) and it uses them to compute the
362sparsity patterns for 
363\verbatim
364        H(y , w , u , ... ) = G[ v[x], y , w , u , ... ]
365\endverbatim
366
367<!-- replace sparse_store_op -->
368The C++ source code corresponding to this operation is
369\verbatim
370        v[x] = y
371\endverbatim
372where v is a VecAD<Base> vector, x is an AD<Base> object,
373and y is AD<Base> or Base objects.
374We define the index corresponding to v[x] by
375\verbatim
376        i_v_x = combined[ arg[0] + i_vec ]
377\endverbatim
378where i_vec is defined under the heading \a arg[1] below:
379
380\tparam Vector_set
381is the type used for vectors of sets. It can be either
382\c sparse_pack, \c sparse_set, or \c sparse_list.
383
384\param op
385is the code corresponding to this operator; i.e., StpvOp or StvvOp
386(only used for error checking).
387
388\param arg
389\n
390\a arg[0]
391is the offset corresponding to this VecAD vector in the combined array.
392\n
393\n
394\a arg[2]
395\n
396The set with index \a arg[2] in \a var_sparsity
397is the sparsity pattern corresponding to y.
398(Note that \a arg[2] > 0 because y is a variable.)
399
400\param num_combined
401is the total number of elements in the VecAD address array.
402
403\param combined
404\a combined [ arg[0] - 1 ]
405is the index of the set in \a vecad_sparsity corresponding
406to the sparsity pattern for the vector v.
407We use the notation i_v below which is defined by
408\verbatim
409        i_v = combined[ \a arg[0] - 1 ]
410\endverbatim
411
412\param var_sparsity
413The set  with index \a arg[2] in \a var_sparsity
414is the sparsity pattern for y.
415This is an input for forward mode operations.
416For reverse mode operations:
417The sparsity pattern for v is added to the spartisy pattern for y.
418
419\param vecad_sparsity
420The set with index \a i_v in \a vecad_sparsity
421is the sparsity pattern for v.
422This is an input for reverse mode operations.
423For forward mode operations, the sparsity pattern for y is added
424to the sparsity pattern for the vector v.
425
426\par Checked Assertions
427\li NumArg(op) == 3
428\li NumRes(op) == 0
429\li 0 <  \a arg[0]
430\li \a arg[0] < \a num_combined
431\li \a arg[2] < \a var_sparsity.n_set()
432\li i_v       < \a vecad_sparsity.n_set()
433<!-- end sparse_store_op -->
434*/
435template <class Vector_set>
436inline void reverse_sparse_jacobian_store_op(
437        OpCode             op              ,
438        const addr_t*      arg             , 
439        size_t             num_combined    ,
440        const size_t*      combined        ,
441        Vector_set&        var_sparsity    ,
442        Vector_set&        vecad_sparsity  )
443{
444        CPPAD_ASSERT_UNKNOWN( NumArg(op) == 3 );
445        CPPAD_ASSERT_UNKNOWN( NumRes(op) == 0 );
446        CPPAD_ASSERT_UNKNOWN( 0 < arg[0] );
447        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_combined );
448        size_t i_v = combined[ arg[0] - 1 ];
449        CPPAD_ASSERT_UNKNOWN( i_v < vecad_sparsity.n_set() );
450        CPPAD_ASSERT_UNKNOWN( size_t(arg[2]) < var_sparsity.n_set() );
451
452        var_sparsity.binary_union(arg[2], arg[2], i_v, vecad_sparsity);
453
454        return;
455}
456
457/*!
458Reverse mode sparsity operations for StpvOp and StvvOp
459
460This routine is given the sparsity patterns for
461G(v[x], y , w , u ... )
462and it uses them to compute the sparsity patterns for
463\verbatim
464        H(y , w , u , ... ) = G[ v[x], y , w , u , ... ]
465\endverbatim
466
467<!-- replace sparse_store_op -->
468The C++ source code corresponding to this operation is
469\verbatim
470        v[x] = y
471\endverbatim
472where v is a VecAD<Base> vector, x is an AD<Base> object,
473and y is AD<Base> or Base objects.
474We define the index corresponding to v[x] by
475\verbatim
476        i_v_x = combined[ arg[0] + i_vec ]
477\endverbatim
478where i_vec is defined under the heading \a arg[1] below:
479
480\tparam Vector_set
481is the type used for vectors of sets. It can be either
482\c sparse_pack, \c sparse_set, or \c sparse_list.
483
484\param op
485is the code corresponding to this operator; i.e., StpvOp or StvvOp
486(only used for error checking).
487
488\param arg
489\n
490\a arg[0]
491is the offset corresponding to this VecAD vector in the combined array.
492\n
493\n
494\a arg[2]
495\n
496The set with index \a arg[2] in \a var_sparsity
497is the sparsity pattern corresponding to y.
498(Note that \a arg[2] > 0 because y is a variable.)
499
500\param num_combined
501is the total number of elements in the VecAD address array.
502
503\param combined
504\a combined [ arg[0] - 1 ]
505is the index of the set in \a vecad_sparsity corresponding
506to the sparsity pattern for the vector v.
507We use the notation i_v below which is defined by
508\verbatim
509        i_v = combined[ \a arg[0] - 1 ]
510\endverbatim
511
512\param var_sparsity
513The set  with index \a arg[2] in \a var_sparsity
514is the sparsity pattern for y.
515This is an input for forward mode operations.
516For reverse mode operations:
517The sparsity pattern for v is added to the spartisy pattern for y.
518
519\param vecad_sparsity
520The set with index \a i_v in \a vecad_sparsity
521is the sparsity pattern for v.
522This is an input for reverse mode operations.
523For forward mode operations, the sparsity pattern for y is added
524to the sparsity pattern for the vector v.
525
526\par Checked Assertions
527\li NumArg(op) == 3
528\li NumRes(op) == 0
529\li 0 <  \a arg[0]
530\li \a arg[0] < \a num_combined
531\li \a arg[2] < \a var_sparsity.n_set()
532\li i_v       < \a vecad_sparsity.n_set()
533<!-- end sparse_store_op -->
534
535\param var_jacobian
536\a var_jacobian[ \a arg[2] ]
537is false (true) if the Jacobian of G with respect to y is always zero
538(may be non-zero).
539
540\param vecad_jacobian
541\a vecad_jacobian[i_v]
542is false (true) if the Jacobian with respect to x is always zero
543(may be non-zero).
544On input, it corresponds to the function G,
545and on output it corresponds to the function H.
546*/
547template <class Vector_set>
548inline void reverse_sparse_hessian_store_op(
549        OpCode             op           ,
550        const addr_t*      arg          , 
551        size_t             num_combined ,
552        const size_t*      combined     ,
553        Vector_set&        var_sparsity ,
554        Vector_set&        vecad_sparsity ,
555        bool*              var_jacobian   ,
556        bool*              vecad_jacobian )
557{
558        CPPAD_ASSERT_UNKNOWN( NumArg(op) == 3 );
559        CPPAD_ASSERT_UNKNOWN( NumRes(op) == 0 );
560        CPPAD_ASSERT_UNKNOWN( 0 < arg[0] );
561        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_combined );
562        size_t i_v = combined[ arg[0] - 1 ];
563        CPPAD_ASSERT_UNKNOWN( i_v < vecad_sparsity.n_set() );
564        CPPAD_ASSERT_UNKNOWN( size_t(arg[2]) < var_sparsity.n_set() );
565
566        var_sparsity.binary_union(arg[2], arg[2], i_v, vecad_sparsity);
567
568        var_jacobian[ arg[2] ] |= vecad_jacobian[i_v];
569
570        return;
571}
572
573} // END_CPPAD_NAMESPACE
574# endif
Note: See TracBrowser for help on using the repository browser.