source: trunk/cppad/local/atomic_base.hpp @ 3805

Last change on this file since 3805 was 3805, checked in by bradbell, 4 years ago

merge to branch: trunk
from repository: https://github.com/coin-or/CppAD
start hash code: 624601e5e709d2872b5873f628fae76d5f575a9b
end hash code: 256147abd1ffb46688686cc6cd53908c9a64d957

commit 256147abd1ffb46688686cc6cd53908c9a64d957
Author: Brad Bell <bradbell@…>
Date: Tue Mar 22 12:56:16 2016 -0700

Change atomic_ode.cpp -> checkpoint_ode.cpp,
and atomic_extended_ode.cpp -> checkpoint_extended_ode.cpp.

commit 42eb37c2c38a377bf165ea2dcdd9721d251da800
Author: Brad Bell <bradbell@…>
Date: Tue Mar 22 09:05:32 2016 -0700

Add the eigen_mat_mul.cpp example (still under construction).


ode.cpp: minor corrections to typos in documentation.

commit 6126e82c2a9b3f10d15ee792bf9d32796eb9561c
Author: Brad Bell <bradbell@…>
Date: Tue Mar 22 04:57:55 2016 -0700

Reduce level of indent on public and private commands in atomic examples.

commit 50b8a65b81a0d4e3b0d9d03b188966eb5b366bec
Author: Brad Bell <bradbell@…>
Date: Tue Mar 22 04:46:22 2016 -0700

  1. Advance version to 20160322.
  2. Move atomic_mat_mul_xam.cpp -> atomic_mat_mul.cpp.

commit bb5a5300f882891dafff643ab6899a131fa477b0
Author: Brad Bell <bradbell@…>
Date: Mon Mar 21 05:40:28 2016 -0700

test_one.sh.in: use path to file instead of dir file.
mat_mul.hpp: add missing virtual (did not matter).

commit 955919bc43857a7f9b953af96be30f0cf524ae1c
Merge: ceb9087 3eece08
Author: Brad Bell <bradbell@…>
Date: Sun Mar 20 12:52:01 2016 -0700

Do not know how master got out of sync with local copy, but this merge
does not see to have any changes.

commit ceb9087c686d99279687134dc0af4da0150a0bc5
Author: Brad Bell <bradbell@…>
Date: Sun Mar 20 11:23:46 2016 -0700

Use utility/set_union to simplify atomic examples and tests.

commit 60e9db895ee11bddfb02b1e2f8a96b9dd8a01e49
Author: Brad Bell <bradbell@…>
Date: Sun Mar 20 11:02:21 2016 -0700

  1. Sort some lists.
  2. Add some references to set_union (missing in previous commit).

commit 3eece08b8e9af5f9cf2633bc300c4e5477acd2b1
Author: Brad Bell <bradbell@…>
Date: Sun Mar 20 10:00:56 2016 -0700

  1. Sort some lists.
  2. Add some references to set_union (missing in previous commit).

commit 309e461bf32d95b9330f65bbffaf4fa24d5f11d6
Author: Brad Bell <bradbell@…>
Date: Sun Mar 20 09:34:14 2016 -0700

Add utility/set_union.hpp.

  • Property svn:keywords set to Id
File size: 43.5 KB
Line 
1// $Id: atomic_base.hpp 3805 2016-03-22 20:01:22Z bradbell $
2# ifndef CPPAD_LOCAL_ATOMIC_BASE_HPP
3# define CPPAD_LOCAL_ATOMIC_BASE_HPP
4
5/* --------------------------------------------------------------------------
6CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-16 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# include <set>
17# include <cppad/local/cppad_assert.hpp>
18// needed before one can use CPPAD_ASSERT_FIRST_CALL_NOT_PARALLEL
19# include <cppad/utility/thread_alloc.hpp>
20
21namespace CppAD { // BEGIN_CPPAD_NAMESPACE
22/*!
23\file atomic_base.hpp
24Base class for atomic user operations.
25*/
26
27template <class Base>
28class atomic_base {
29// ===================================================================
30public:
31        enum option_enum {
32                pack_sparsity_enum   ,
33                bool_sparsity_enum   ,
34                set_sparsity_enum
35        };
36private:
37        // ------------------------------------------------------
38        // constants
39        //
40        /// index of this object in class_object
41        const size_t index_;
42
43        // -----------------------------------------------------
44        // variables
45        //
46        /// sparsity pattern this object is currently using
47        /// (set by constructor and option member functions)
48        option_enum sparsity_;
49
50        /// temporary work space used afun, declared here to avoid memory
51        /// allocation/deallocation for each call to afun
52        vector<bool>  afun_vx_[CPPAD_MAX_NUM_THREADS];
53        vector<bool>  afun_vy_[CPPAD_MAX_NUM_THREADS];
54        vector<Base>  afun_tx_[CPPAD_MAX_NUM_THREADS];
55        vector<Base>  afun_ty_[CPPAD_MAX_NUM_THREADS];
56
57        // -----------------------------------------------------
58        // static member functions
59        //
60        /// List of all the object in this class
61        static std::vector<atomic_base *>& class_object(void)
62        {       CPPAD_ASSERT_FIRST_CALL_NOT_PARALLEL;
63                static std::vector<atomic_base *> list_;
64                return list_;
65        }
66        /// List of names for each object in this class
67        static std::vector<std::string>& class_name(void)
68        {       CPPAD_ASSERT_FIRST_CALL_NOT_PARALLEL;
69                static std::vector<std::string> list_;
70                return list_;
71        }
72        // =====================================================================
73public:
74        // -----------------------------------------------------
75        // member functions not in user API
76        //
77        /// current sparsity setting
78        option_enum sparsity(void) const
79        {       return sparsity_; }
80
81        /// Name corresponding to a base_atomic object
82        const std::string& afun_name(void) const
83        {       return class_name()[index_]; }
84/*
85$begin atomic_ctor$$
86$spell
87        enum
88        sq
89        std
90        afun
91        arg
92        CppAD
93        bool
94        ctor
95        const
96        mat_mul_xam.cpp
97        hpp
98$$
99
100$section Atomic Function Constructor$$
101
102$head Syntax$$
103$icode%atomic_user afun%(%ctor_arg_list%)
104%$$
105$codei%atomic_base<%Base%>(%name%, %sparsity%)
106%$$
107
108$head atomic_user$$
109
110$subhead ctor_arg_list$$
111Is a list of arguments for the $icode atomic_user$$ constructor.
112
113$subhead afun$$
114The object $icode afun$$ must stay in scope for as long
115as the corresponding atomic function is used.
116This includes use by any $cref/ADFun<Base>/ADFun/$$ that
117has this $icode atomic_user$$ operation in its
118$cref/operation sequence/glossary/Operation/Sequence/$$.
119
120$subhead Implementation$$
121The user defined $icode atomic_user$$ class is a publicly derived class of
122$codei%atomic_base<%Base%>%$$.
123It should be declared as follows:
124$codei%
125        class %atomic_user% : public CppAD::atomic_base<%Base%> {
126        public:
127                %atomic_user%(%ctor_arg_list%) : atomic_base<%Base%>(%name%, %sparsity%)
128        %...%
129        };
130%$$
131where $icode ...$$
132denotes the rest of the implementation of the derived class.
133This includes completing the constructor and
134all the virtual functions that have their
135$code atomic_base$$ implementations replaced by
136$icode atomic_user$$ implementations.
137
138$head atomic_base$$
139
140$subhead Restrictions$$
141The $code atomic_base$$ constructor cannot be called in
142$cref/parallel/ta_in_parallel/$$ mode.
143
144$subhead Base$$
145The template parameter determines the
146$icode Base$$ type for this $codei%AD<%Base%>%$$ atomic operation.
147
148$subhead name$$
149This $code atomic_base$$ constructor argument has the following prototype
150$codei%
151        const std::string& %name%
152%$$
153It is the name for this atomic function and is used for error reporting.
154The suggested value for $icode name$$ is $icode afun$$ or $icode atomic_user$$,
155i.e., the name of the corresponding atomic object or class.
156
157$subhead sparsity$$
158This $code atomic_base$$ constructor argument has prototype
159$codei%
160        atomic_base<%Base%>::option_enum %sparsity%
161%$$
162The current $icode sparsity$$ for an $code atomic_base$$ object
163determines which type of sparsity patterns it uses
164and its value is one of the following:
165$table
166$icode sparsity$$   $cnext sparsity patterns $rnext
167$codei%atomic_base<%Base%>::pack_sparsity_enum%$$ $pre  $$ $cnext
168        $cref/vectorBool/CppAD_vector/vectorBool/$$
169$rnext
170$codei%atomic_base<%Base%>::bool_sparsity_enum%$$ $pre  $$ $cnext
171        $cref/vector/CppAD_vector/$$$code <bool>$$
172$rnext
173$codei%atomic_base<%Base%>::set_sparsity_enum%$$ $pre  $$ $cnext
174        $cref/vector/CppAD_vector/$$$code <std::set<std::size_t> >$$
175$tend
176There is a default value for $icode sparsity$$ if it is not
177included in the constructor (which may be either the bool or set option).
178
179$head Examples$$
180
181$subhead Define Constructor$$
182The following are links to user atomic function constructor definitions:
183$cref%get_started.cpp%atomic_get_started.cpp%Constructor%$$,
184$cref%norm_sq.cpp%atomic_norm_sq.cpp%Constructor%$$,
185$cref%reciprocal.cpp%atomic_reciprocal.cpp%Constructor%$$,
186$cref%tangent.cpp%atomic_tangent.cpp%Constructor%$$,
187$cref%mat_mul.hpp%atomic_mat_mul.hpp%Constructor%$$.
188
189$subhead Use Constructor$$
190The following are links to user atomic function constructor uses:
191$cref%get_started.cpp%atomic_get_started.cpp%Use Atomic Function%Constructor%$$,
192$cref%norm_sq.cpp%atomic_norm_sq.cpp%Use Atomic Function%Constructor%$$,
193$cref%reciprocal.cpp%atomic_reciprocal.cpp%Use Atomic Function%Constructor%$$,
194$cref%tangent.cpp%atomic_tangent.cpp%Use Atomic Function%Constructor%$$,
195$cref%mat_mul.cpp%atomic_mat_mul.cpp%Use Atomic Function%Constructor%$$.
196
197
198$end
199*/
200/*!
201Base class for atomic_user functions.
202
203\tparam Base
204This class is used for defining an AD<Base> atomic operation y = f(x).
205*/
206/// make sure user does not invoke the default constructor
207atomic_base(void)
208{       CPPAD_ASSERT_KNOWN(false,
209                "Attempt to use the atomic_base default constructor"
210        );
211}
212/*!
213Constructor
214
215\param name
216name used for error reporting
217
218\param sparsity [in]
219what type of sparsity patterns are computed by this function,
220bool_sparsity_enum or set_sparsity_enum. Default value is
221bool sparsity patterns.
222*/
223atomic_base(
224                const std::string&     name,
225                option_enum            sparsity = bool_sparsity_enum
226) :
227index_( class_object().size() )     ,
228sparsity_( sparsity )
229{       CPPAD_ASSERT_KNOWN(
230                ! thread_alloc::in_parallel() ,
231                "atomic_base: constructor cannot be called in parallel mode."
232        );
233        class_object().push_back(this);
234        class_name().push_back(name);
235        CPPAD_ASSERT_UNKNOWN( class_object().size() == class_name().size() );
236}
237/// destructor informs CppAD that this atomic function with this index
238/// has dropped out of scope by setting its pointer to null
239virtual ~atomic_base(void)
240{       CPPAD_ASSERT_UNKNOWN( class_object().size() > index_ );
241        // change object pointer to null, but leave name for error reporting
242        class_object()[index_] = CPPAD_NULL;
243}
244/// atomic_base function object corresponding to a certain index
245static atomic_base* class_object(size_t index)
246{       CPPAD_ASSERT_UNKNOWN( class_object().size() > index );
247        return class_object()[index];
248}
249/// atomic_base function name corresponding to a certain index
250static const std::string& class_name(size_t index)
251{       CPPAD_ASSERT_UNKNOWN( class_name().size() > index );
252        return class_name()[index];
253}
254/*
255$begin atomic_option$$
256$spell
257        sq
258        enum
259        afun
260        bool
261        CppAD
262        std
263        typedef
264$$
265
266$section Set Atomic Function Options$$
267
268$head Syntax$$
269$icode%afun%.option(%option_value%)%$$
270These settings do not apply to individual $icode afun$$ calls,
271but rather all subsequent uses of the corresponding atomic operation
272in an $cref ADFun$$ object.
273
274$head atomic_sparsity$$
275Note that, if you use $cref optimize$$, these sparsity patterns are used
276to determine the $cref/dependency/dependency.cpp/$$ relationship between
277argument and result variables.
278
279$subhead pack_sparsity_enum$$
280If $icode option_value$$ is $codei%atomic_base<%Base%>::pack_sparsity_enum%$$,
281then the type used by $icode afun$$ for
282$cref/sparsity patterns/glossary/Sparsity Pattern/$$,
283(after the option is set) will be
284$codei%
285        typedef CppAD::vectorBool %atomic_sparsity%
286%$$
287If $icode r$$ is a sparsity pattern
288for a matrix $latex R \in B^{p \times q}$$:
289$icode%r%.size() == %p% * %q%$$.
290
291$subhead bool_sparsity_enum$$
292If $icode option_value$$ is $codei%atomic_base<%Base%>::bool_sparsity_enum%$$,
293then the type used by $icode afun$$ for
294$cref/sparsity patterns/glossary/Sparsity Pattern/$$,
295(after the option is set) will be
296$codei%
297        typedef CppAD::vector<bool> %atomic_sparsity%
298%$$
299If $icode r$$ is a sparsity pattern
300for a matrix $latex R \in B^{p \times q}$$:
301$icode%r%.size() == %p% * %q%$$.
302
303$subhead set_sparsity_enum$$
304If $icode option_value$$ is $icode%atomic_base<%Base%>::set_sparsity_enum%$$,
305then the type used by $icode afun$$ for
306$cref/sparsity patterns/glossary/Sparsity Pattern/$$,
307(after the option is set) will be
308$codei%
309        typedef CppAD::vector< std::set<size_t> > %atomic_sparsity%
310%$$
311If $icode r$$ is a sparsity pattern
312for a matrix $latex R \in B^{p \times q}$$:
313$icode%r%.size() == %p%$$, and for $latex i = 0 , \ldots , p-1$$,
314the elements of $icode%r%[%i%]%$$ are between zero and $latex q-1$$ inclusive.
315
316$end
317*/
318void option(enum option_enum option_value)
319{       switch( option_value )
320        {       case pack_sparsity_enum:
321                case bool_sparsity_enum:
322                case set_sparsity_enum:
323                sparsity_ = option_value;
324                break;
325
326                default:
327                CPPAD_ASSERT_KNOWN(
328                        false,
329                        "atoic_base::option: option_value is not valid"
330                );
331        }
332        return;
333}
334/*
335-----------------------------------------------------------------------------
336$begin atomic_afun$$
337
338$spell
339        sq
340        mul
341        afun
342        const
343        CppAD
344        mat_mul.cpp
345$$
346
347$section Using AD Version of Atomic Function$$
348
349$head Syntax$$
350$icode%afun%(%ax%, %ay%)%$$
351
352$head Purpose$$
353Given $icode ax$$,
354this call computes the corresponding value of $icode ay$$.
355If $codei%AD<%Base%>%$$ operations are being recorded,
356it enters the computation as an atomic operation in the recording;
357see $cref/start recording/Independent/Start Recording/$$.
358
359$head ADVector$$
360The type $icode ADVector$$ must be a
361$cref/simple vector class/SimpleVector/$$ with elements of type
362$codei%AD<%Base%>%$$; see $cref/Base/atomic_ctor/atomic_base/Base/$$.
363
364$head afun$$
365is a $cref/atomic_user/atomic_ctor/atomic_user/$$ object
366and this $icode afun$$ function call is implemented by the
367$cref/atomic_base/atomic_ctor/atomic_base/$$ class.
368
369$head ax$$
370This argument has prototype
371$codei%
372        const %ADVector%& %ax%
373%$$
374and size must be equal to $icode n$$.
375It specifies vector $latex x \in B^n$$
376at which an $codei%AD<%Base%>%$$ version of
377$latex y = f(x)$$ is to be evaluated; see
378$cref/Base/atomic_ctor/atomic_base/Base/$$.
379
380$head ay$$
381This argument has prototype
382$codei%
383        %ADVector%& %ay%
384%$$
385and size must be equal to $icode m$$.
386The input values of its elements
387are not specified (must not matter).
388Upon return, it is an $codei%AD<%Base%>%$$ version of
389$latex y = f(x)$$.
390
391$head Examples$$
392The following files contain example uses of
393the AD version of atomic functions during recording:
394$cref%get_started.cpp%atomic_get_started.cpp%Use Atomic Function%Recording%$$,
395$cref%norm_sq.cpp%atomic_norm_sq.cpp%Use Atomic Function%Recording%$$,
396$cref%reciprocal.cpp%atomic_reciprocal.cpp%Use Atomic Function%Recording%$$,
397$cref%tangent.cpp%atomic_tangent.cpp%Use Atomic Function%Recording%$$,
398$cref%mat_mul.cpp%atomic_mat_mul.cpp%Use Atomic Function%Recording%$$.
399
400$end
401-----------------------------------------------------------------------------
402*/
403/*!
404Implement the user call to <tt>afun(ax, ay)</tt> and old_atomic call to
405<tt>afun(ax, ay, id)</tt>.
406
407\tparam ADVector
408A simple vector class with elements of type <code>AD<Base></code>.
409
410\param id
411optional extra information vector that is just passed through by CppAD,
412and used by old_atomic derived class (not other derived classes).
413This is an extra parameter to the virtual callbacks for old_atomic;
414see the set_id member function.
415
416\param ax
417is the argument vector for this call,
418<tt>ax.size()</tt> determines the number of arguments.
419
420\param ay
421is the result vector for this call,
422<tt>ay.size()</tt> determines the number of results.
423*/
424template <class ADVector>
425void operator()(
426        const ADVector&  ax     ,
427              ADVector&  ay     ,
428        size_t           id = 0 )
429{       size_t i, j;
430        size_t n = ax.size();
431        size_t m = ay.size();
432# ifndef NDEBUG
433        bool ok;
434        std::string msg = "atomic_base: " + afun_name() + ".eval: ";
435        if( (n == 0) | (m == 0) )
436        {       msg += "ax.size() or ay.size() is zero";
437                CPPAD_ASSERT_KNOWN(false, msg.c_str() );
438        }
439# endif
440        size_t thread = thread_alloc::thread_num();
441        vector <Base>& tx  = afun_tx_[thread];
442        vector <Base>& ty  = afun_ty_[thread];
443        vector <bool>& vx  = afun_vx_[thread];
444        vector <bool>& vy  = afun_vy_[thread];
445        //
446        if( vx.size() != n )
447        {       vx.resize(n);
448                tx.resize(n);
449        }
450        if( vy.size() != m )
451        {       vy.resize(m);
452                ty.resize(m);
453        }
454        //
455        // Determine tape corresponding to variables in ax
456        tape_id_t     tape_id  = 0;
457        ADTape<Base>* tape     = CPPAD_NULL;
458        for(j = 0; j < n; j++)
459        {       tx[j]  = ax[j].value_;
460                vx[j]  = Variable( ax[j] );
461                if( vx[j] )
462                {
463                        if( tape_id == 0 )
464                        {       tape    = ax[j].tape_this();
465                                tape_id = ax[j].tape_id_;
466                                CPPAD_ASSERT_UNKNOWN( tape != CPPAD_NULL );
467                        }
468# ifndef NDEBUG
469                        if( tape_id != ax[j].tape_id_ )
470                        {       msg += afun_name() +
471                                ": ax contains variables from different threads.";
472                                CPPAD_ASSERT_KNOWN(false, msg.c_str());
473                        }
474# endif
475                }
476        }
477        // Use zero order forward mode to compute values
478        size_t p = 0, q = 0;
479        set_id(id);
480# ifdef NDEBUG
481        forward(p, q, vx, vy, tx, ty);
482# else
483        ok = forward(p, q, vx, vy, tx, ty);
484        if( ! ok )
485        {       msg += afun_name() + ": ok is false for "
486                        "zero order forward mode calculation.";
487                CPPAD_ASSERT_KNOWN(false, msg.c_str());
488        }
489# endif
490        bool record_operation = false;
491        for(i = 0; i < m; i++)
492        {
493                // pass back values
494                ay[i].value_ = ty[i];
495
496                // initialize entire vector parameters (not in tape)
497                ay[i].tape_id_ = 0;
498                ay[i].taddr_   = 0;
499
500                // we need to record this operation if
501                // any of the eleemnts of ay are variables,
502                record_operation |= vy[i];
503        }
504# ifndef NDEBUG
505        if( record_operation & (tape == CPPAD_NULL) )
506        {       msg +=
507                "all elements of vx are false but vy contains a true element";
508                CPPAD_ASSERT_KNOWN(false, msg.c_str() );
509        }
510# endif
511        // if tape is not null, ay is on the tape
512        if( record_operation )
513        {
514                // Operator that marks beginning of this atomic operation
515                CPPAD_ASSERT_UNKNOWN( NumRes(UserOp) == 0 );
516                CPPAD_ASSERT_UNKNOWN( NumArg(UserOp) == 4 );
517                tape->Rec_.PutArg(index_, id, n, m);
518                tape->Rec_.PutOp(UserOp);
519
520                // Now put n operators, one for each element of argument vector
521                CPPAD_ASSERT_UNKNOWN( NumRes(UsravOp) == 0 );
522                CPPAD_ASSERT_UNKNOWN( NumRes(UsrapOp) == 0 );
523                CPPAD_ASSERT_UNKNOWN( NumArg(UsravOp) == 1 );
524                CPPAD_ASSERT_UNKNOWN( NumArg(UsrapOp) == 1 );
525                for(j = 0; j < n; j++)
526                {       if( vx[j] )
527                        {       // information for an argument that is a variable
528                                tape->Rec_.PutArg(ax[j].taddr_);
529                                tape->Rec_.PutOp(UsravOp);
530                        }
531                        else
532                        {       // information for an argument that is parameter
533                                addr_t par = tape->Rec_.PutPar(ax[j].value_);
534                                tape->Rec_.PutArg(par);
535                                tape->Rec_.PutOp(UsrapOp);
536                        }
537                }
538
539                // Now put m operators, one for each element of result vector
540                CPPAD_ASSERT_UNKNOWN( NumArg(UsrrpOp) == 1 );
541                CPPAD_ASSERT_UNKNOWN( NumRes(UsrrpOp) == 0 );
542                CPPAD_ASSERT_UNKNOWN( NumArg(UsrrvOp) == 0 );
543                CPPAD_ASSERT_UNKNOWN( NumRes(UsrrvOp) == 1 );
544                for(i = 0; i < m; i++)
545                {       if( vy[i] )
546                        {       ay[i].taddr_    = tape->Rec_.PutOp(UsrrvOp);
547                                ay[i].tape_id_  = tape_id;
548                        }
549                        else
550                        {       addr_t par = tape->Rec_.PutPar(ay[i].value_);
551                                tape->Rec_.PutArg(par);
552                                tape->Rec_.PutOp(UsrrpOp);
553                        }
554                }
555
556                // Put a duplicate UserOp at end of UserOp sequence
557                tape->Rec_.PutArg(index_, id, n, m);
558                tape->Rec_.PutOp(UserOp);
559        }
560        return;
561}
562/*
563-----------------------------------------------------------------------------
564$begin atomic_forward$$
565$spell
566        sq
567        mul.hpp
568        hes
569        afun
570        vx
571        vy
572        ty
573        Taylor
574        const
575        CppAD
576        bool
577$$
578
579$section Atomic Forward Mode$$
580$mindex callback virtual$$
581
582
583$head Syntax$$
584$icode%ok% = %afun%.forward(%p%, %q%, %vx%, %vy%, %tx%, %ty%)%$$
585
586$head Purpose$$
587This virtual function is used by $cref atomic_afun$$
588to evaluate function values.
589It is also used buy $cref/forward/Forward/$$
590to compute function vales and derivatives.
591
592$head Implementation$$
593This virtual function must be defined by the
594$cref/atomic_user/atomic_ctor/atomic_user/$$ class.
595It can just return $icode%ok% == false%$$
596(and not compute anything) for values
597of $icode%q% > 0%$$ that are greater than those used by your
598$cref/forward/Forward/$$ mode calculations.
599
600$head p$$
601The argument $icode p$$ has prototype
602$codei%
603        size_t %p%
604%$$
605It specifies the lowest order Taylor coefficient that we are evaluating.
606During calls to $cref atomic_afun$$, $icode%p% == 0%$$.
607
608$head q$$
609The argument $icode q$$ has prototype
610$codei%
611        size_t %q%
612%$$
613It specifies the highest order Taylor coefficient that we are evaluating.
614During calls to $cref atomic_afun$$, $icode%q% == 0%$$.
615
616$head vx$$
617The $code forward$$ argument $icode vx$$ has prototype
618$codei%
619        const CppAD::vector<bool>& %vx%
620%$$
621The case $icode%vx%.size() > 0%$$ only occurs while evaluating a call to
622$cref atomic_afun$$.
623In this case,
624$icode%p% == %q% == 0%$$,
625$icode%vx%.size() == %n%$$, and
626for $latex j = 0 , \ldots , n-1$$,
627$icode%vx%[%j%]%$$ is true if and only if
628$icode%ax%[%j%]%$$ is a $cref/variable/glossary/Variable/$$
629in the corresponding call to
630$codei%
631        %afun%(%ax%, %ay%)
632%$$
633If $icode%vx%.size() == 0%$$,
634then $icode%vy%.size() == 0%$$ and neither of these vectors
635should be used.
636
637$head vy$$
638The $code forward$$ argument $icode vy$$ has prototype
639$codei%
640        CppAD::vector<bool>& %vy%
641%$$
642If $icode%vy%.size() == 0%$$, it should not be used.
643Otherwise,
644$icode%q% == 0%$$ and $icode%vy%.size() == %m%$$.
645The input values of the elements of $icode vy$$
646are not specified (must not matter).
647Upon return, for $latex j = 0 , \ldots , m-1$$,
648$icode%vy%[%i%]%$$ is true if and only if
649$icode%ay%[%i%]%$$ is a variable
650(CppAD uses $icode vy$$ to reduce the necessary computations).
651
652$head tx$$
653The argument $icode tx$$ has prototype
654$codei%
655        const CppAD::vector<%Base%>& %tx%
656%$$
657and $icode%tx%.size() == (%q%+1)*%n%$$.
658For $latex j = 0 , \ldots , n-1$$ and $latex k = 0 , \ldots , q$$,
659we use the Taylor coefficient notation
660$latex \[
661\begin{array}{rcl}
662        x_j^k    & = & tx [ j * ( q + 1 ) + k ]
663        \\
664        X_j (t)  & = & x_j^0 + x_j^1 t^1 + \cdots + x_j^q t^q
665\end{array}
666\] $$
667Note that superscripts represent an index for $latex x_j^k$$
668and an exponent for $latex t^k$$.
669Also note that the Taylor coefficients for $latex X(t)$$ correspond
670to the derivatives of $latex X(t)$$ at $latex t = 0$$ in the following way:
671$latex \[
672        x_j^k = \frac{1}{ k ! } X_j^{(k)} (0)
673\] $$
674
675$head ty$$
676The argument $icode ty$$ has prototype
677$codei%
678        CppAD::vector<%Base%>& %ty%
679%$$
680and $icode%tx%.size() == (%q%+1)*%m%$$.
681Upon return,
682For $latex i = 0 , \ldots , m-1$$ and $latex k = 0 , \ldots , q$$,
683$latex \[
684\begin{array}{rcl}
685        Y_i (t)  & = & f_i [ X(t) ]
686        \\
687        Y_i (t)  & = & y_i^0 + y_i^1 t^1 + \cdots + y_i^q t^q + o ( t^q )
688        \\
689        ty [ i * ( q + 1 ) + k ] & = & y_i^k
690\end{array}
691\] $$
692where $latex o( t^q ) / t^q \rightarrow 0$$ as $latex t \rightarrow 0$$.
693Note that superscripts represent an index for $latex y_j^k$$
694and an exponent for $latex t^k$$.
695Also note that the Taylor coefficients for $latex Y(t)$$ correspond
696to the derivatives of $latex Y(t)$$ at $latex t = 0$$ in the following way:
697$latex \[
698        y_j^k = \frac{1}{ k ! } Y_j^{(k)} (0)
699\] $$
700If $latex p > 0$$,
701for $latex i = 0 , \ldots , m-1$$ and $latex k = 0 , \ldots , p-1$$,
702the input of $icode ty$$ satisfies
703$latex \[
704        ty [ i * ( q + 1 ) + k ] = y_i^k
705\]$$
706and hence the corresponding elements need not be recalculated.
707
708$head ok$$
709If the required results are calculated, $icode ok$$ should be true.
710Otherwise, it should be false.
711
712$head Discussion$$
713For example, suppose that $icode%q% == 2%$$,
714and you know how to compute the function $latex f(x)$$,
715its first derivative $latex f^{(1)} (x)$$,
716and it component wise Hessian $latex f_i^{(2)} (x)$$.
717Then you can compute $icode ty$$ using the following formulas:
718$latex \[
719\begin{array}{rcl}
720y_i^0 & = & Y(0)
721        = f_i ( x^0 )
722\\
723y_i^1 & = & Y^{(1)} ( 0 )
724        = f_i^{(1)} ( x^0 ) X^{(1)} ( 0 )
725        = f_i^{(1)} ( x^0 ) x^1
726\\
727y_i^2
728& = & \frac{1}{2 !} Y^{(2)} (0)
729\\
730& = & \frac{1}{2} X^{(1)} (0)^\R{T} f_i^{(2)} ( x^0 ) X^{(1)} ( 0 )
731  +   \frac{1}{2} f_i^{(1)} ( x^0 ) X^{(2)} ( 0 )
732\\
733& = & \frac{1}{2} (x^1)^\R{T} f_i^{(2)} ( x^0 ) x^1
734  +    f_i^{(1)} ( x^0 ) x^2
735\end{array}
736\] $$
737For $latex i = 0 , \ldots , m-1$$, and $latex k = 0 , 1 , 2$$,
738$latex \[
739        ty [ i * (q + 1) + k ] = y_i^k
740\] $$
741
742$children%
743        example/atomic/forward.cpp
744%$$
745$head Examples$$
746The file $cref atomic_forward.cpp$$ contains an example and test
747that uses this routine.
748It returns true if the test passes and false if it fails.
749
750$end
751-----------------------------------------------------------------------------
752*/
753/*!
754Link from atomic_base to forward mode
755
756\param p [in]
757lowerest order for this forward mode calculation.
758
759\param q [in]
760highest order for this forward mode calculation.
761
762\param vx [in]
763if size not zero, which components of \c x are variables
764
765\param vy [out]
766if size not zero, which components of \c y are variables
767
768\param tx [in]
769Taylor coefficients corresponding to \c x for this calculation.
770
771\param ty [out]
772Taylor coefficient corresponding to \c y for this calculation
773
774See the forward mode in user's documentation for base_atomic
775*/
776virtual bool forward(
777        size_t                    p  ,
778        size_t                    q  ,
779        const vector<bool>&       vx ,
780              vector<bool>&       vy ,
781        const vector<Base>&       tx ,
782              vector<Base>&       ty )
783{       return false; }
784/*
785-----------------------------------------------------------------------------
786$begin atomic_reverse$$
787$spell
788        sq
789        mul.hpp
790        afun
791        ty
792        px
793        py
794        Taylor
795        const
796        CppAD
797$$
798
799$section Atomic Reverse Mode$$
800$spell
801        bool
802$$
803
804$head Syntax$$
805$icode%ok% = %afun%.reverse(%q%, %tx%, %ty%, %px%, %py%)%$$
806
807$head Purpose$$
808This function is used by $cref/reverse/Reverse/$$
809to compute derivatives.
810
811$head Implementation$$
812If you are using
813$cref/reverse/Reverse/$$ mode,
814this virtual function must be defined by the
815$cref/atomic_user/atomic_ctor/atomic_user/$$ class.
816It can just return $icode%ok% == false%$$
817(and not compute anything) for values
818of $icode q$$ that are greater than those used by your
819$cref/reverse/Reverse/$$ mode calculations.
820
821$head q$$
822The argument $icode q$$ has prototype
823$codei%
824        size_t %q%
825%$$
826It specifies the highest order Taylor coefficient that
827computing the derivative of.
828
829$head tx$$
830The argument $icode tx$$ has prototype
831$codei%
832        const CppAD::vector<%Base%>& %tx%
833%$$
834and $icode%tx%.size() == (%q%+1)*%n%$$.
835For $latex j = 0 , \ldots , n-1$$ and $latex k = 0 , \ldots , q$$,
836we use the Taylor coefficient notation
837$latex \[
838\begin{array}{rcl}
839        x_j^k    & = & tx [ j * ( q + 1 ) + k ]
840        \\
841        X_j (t)  & = & x_j^0 + x_j^1 t^1 + \cdots + x_j^q t^q
842\end{array}
843\] $$
844Note that superscripts represent an index for $latex x_j^k$$
845and an exponent for $latex t^k$$.
846Also note that the Taylor coefficients for $latex X(t)$$ correspond
847to the derivatives of $latex X(t)$$ at $latex t = 0$$ in the following way:
848$latex \[
849        x_j^k = \frac{1}{ k ! } X_j^{(k)} (0)
850\] $$
851
852$head ty$$
853The argument $icode ty$$ has prototype
854$codei%
855        const CppAD::vector<%Base%>& %ty%
856%$$
857and $icode%tx%.size() == (%q%+1)*%m%$$.
858For $latex i = 0 , \ldots , m-1$$ and $latex k = 0 , \ldots , q$$,
859we use the Taylor coefficient notation
860$latex \[
861\begin{array}{rcl}
862        Y_i (t)  & = & f_i [ X(t) ]
863        \\
864        Y_i (t)  & = & y_i^0 + y_i^1 t^1 + \cdots + y_i^q t^q + o ( t^q )
865        \\
866        y_i^k    & = & ty [ i * ( q + 1 ) + k ]
867\end{array}
868\] $$
869where $latex o( t^q ) / t^q \rightarrow 0$$ as $latex t \rightarrow 0$$.
870Note that superscripts represent an index for $latex y_j^k$$
871and an exponent for $latex t^k$$.
872Also note that the Taylor coefficients for $latex Y(t)$$ correspond
873to the derivatives of $latex Y(t)$$ at $latex t = 0$$ in the following way:
874$latex \[
875        y_j^k = \frac{1}{ k ! } Y_j^{(k)} (0)
876\] $$
877
878
879$head F$$
880We use the notation $latex \{ x_j^k \} \in B^{n \times (q+1)}$$ for
881$latex \[
882        \{ x_j^k \W{:} j = 0 , \ldots , n-1, k = 0 , \ldots , q \}
883\]$$
884We use the notation $latex \{ y_i^k \} \in B^{m \times (q+1)}$$ for
885$latex \[
886        \{ y_i^k \W{:} i = 0 , \ldots , m-1, k = 0 , \ldots , q \}
887\]$$
888We define the function
889$latex F : B^{n \times (q+1)} \rightarrow B^{m \times (q+1)}$$ by
890$latex \[
891        y_i^k = F_i^k [ \{ x_j^k \} ]
892\] $$
893Note that
894$latex \[
895        F_i^0 ( \{ x_j^k \} ) = f_i ( X(0) )  = f_i ( x^0 )
896\] $$
897We also note that
898$latex F_i^\ell ( \{ x_j^k \} )$$ is a function of
899$latex x^0 , \ldots , x^\ell$$
900and is determined by the derivatives of $latex f_i (x)$$
901up to order $latex \ell$$.
902
903
904$head G, H$$
905We use $latex G : B^{m \times (q+1)} \rightarrow B$$
906to denote an arbitrary scalar valued function of $latex \{ y_i^k \}$$.
907We use $latex H : B^{n \times (q+1)} \rightarrow B$$
908defined by
909$latex \[
910        H ( \{ x_j^k \} ) = G[ F( \{ x_j^k \} ) ]
911\] $$
912
913$head py$$
914The argument $icode py$$ has prototype
915$codei%
916        const CppAD::vector<%Base%>& %py%
917%$$
918and $icode%py%.size() == m * (%q%+1)%$$.
919For $latex i = 0 , \ldots , m-1$$, $latex k = 0 , \ldots , q$$,
920$latex \[
921        py[ i * (q + 1 ) + k ] = \partial G / \partial y_i^k
922\] $$
923
924$subhead px$$
925The $icode px$$ has prototype
926$codei%
927        CppAD::vector<%Base%>& %px%
928%$$
929and $icode%px%.size() == n * (%q%+1)%$$.
930The input values of the elements of $icode px$$
931are not specified (must not matter).
932Upon return,
933for $latex j = 0 , \ldots , n-1$$ and $latex \ell = 0 , \ldots , q$$,
934$latex \[
935\begin{array}{rcl}
936px [ j * (q + 1) + \ell ] & = & \partial H / \partial x_j^\ell
937\\
938& = &
939( \partial G / \partial \{ y_i^k \} ) \cdot
940        ( \partial \{ y_i^k \} / \partial x_j^\ell )
941\\
942& = &
943\sum_{k=0}^q
944\sum_{i=0}^{m-1}
945( \partial G / \partial y_i^k ) ( \partial y_i^k / \partial x_j^\ell )
946\\
947& = &
948\sum_{k=\ell}^q
949\sum_{i=0}^{m-1}
950py[ i * (q + 1 ) + k ] ( \partial F_i^k / \partial x_j^\ell )
951\end{array}
952\] $$
953Note that we have used the fact that for $latex k < \ell$$,
954$latex \partial F_i^k / \partial x_j^\ell = 0$$.
955
956$head ok$$
957The return value $icode ok$$ has prototype
958$codei%
959        bool %ok%
960%$$
961If it is $code true$$, the corresponding evaluation succeeded,
962otherwise it failed.
963
964$children%
965        example/atomic/reverse.cpp
966%$$
967$head Examples$$
968The file $cref atomic_forward.cpp$$ contains an example and test
969that uses this routine.
970It returns true if the test passes and false if it fails.
971
972$end
973-----------------------------------------------------------------------------
974*/
975/*!
976Link from reverse mode sweep to users routine.
977
978\param q [in]
979highest order for this reverse mode calculation.
980
981\param tx [in]
982Taylor coefficients corresponding to \c x for this calculation.
983
984\param ty [in]
985Taylor coefficient corresponding to \c y for this calculation
986
987\param px [out]
988Partials w.r.t. the \c x Taylor coefficients.
989
990\param py [in]
991Partials w.r.t. the \c y Taylor coefficients.
992
993See atomic_reverse mode use documentation
994*/
995virtual bool reverse(
996        size_t                    q  ,
997        const vector<Base>&       tx ,
998        const vector<Base>&       ty ,
999              vector<Base>&       px ,
1000        const vector<Base>&       py )
1001{       return false; }
1002/*
1003-------------------------------------- ---------------------------------------
1004$begin atomic_for_sparse_jac$$
1005$spell
1006        sq
1007        mul.hpp
1008        afun
1009        Jacobian
1010        jac
1011        const
1012        CppAD
1013        std
1014        bool
1015        std
1016$$
1017
1018$section Atomic Forward Jacobian Sparsity Patterns$$
1019
1020$head Syntax$$
1021$icode%ok% = %afun%.for_sparse_jac(%q%, %r%, %s%)%$$
1022
1023$head Purpose$$
1024This function is used by $cref ForSparseJac$$ to compute
1025Jacobian sparsity patterns.
1026For a fixed matrix $latex R \in B^{n \times q}$$,
1027the Jacobian of $latex f( x + R * u)$$ with respect to $latex u \in B^q$$ is
1028$latex \[
1029        S(x) = f^{(1)} (x) * R
1030\] $$
1031Given a $cref/sparsity pattern/glossary/Sparsity Pattern/$$ for $latex R$$,
1032$code for_sparse_jac$$ computes a sparsity pattern for $latex S(x)$$.
1033
1034$head Implementation$$
1035If you are using
1036$cref ForSparseJac$$,
1037$cref ForSparseHes$$, or
1038$cref RevSparseHes$$,
1039this virtual function must be defined by the
1040$cref/atomic_user/atomic_ctor/atomic_user/$$ class.
1041
1042$subhead q$$
1043The argument $icode q$$ has prototype
1044$codei%
1045     size_t %q%
1046%$$
1047It specifies the number of columns in
1048$latex R \in B^{n \times q}$$ and the Jacobian
1049$latex S(x) \in B^{m \times q}$$.
1050
1051$subhead r$$
1052This argument has prototype
1053$codei%
1054     const %atomic_sparsity%& %r%
1055%$$
1056and is a $cref/atomic_sparsity/atomic_option/atomic_sparsity/$$ pattern for
1057$latex R \in B^{n \times q}$$.
1058
1059$subhead s$$
1060This argument has prototype
1061$codei%
1062        %atomic_sparsity%& %s%
1063%$$
1064The input values of its elements
1065are not specified (must not matter).
1066Upon return, $icode s$$ is a
1067$cref/atomic_sparsity/atomic_option/atomic_sparsity/$$ pattern for
1068$latex S(x) \in B^{m \times q}$$.
1069
1070$head ok$$
1071The return value $icode ok$$ has prototype
1072$codei%
1073        bool %ok%
1074%$$
1075If it is $code true$$, the corresponding evaluation succeeded,
1076otherwise it failed.
1077
1078$children%
1079        example/atomic/for_sparse_jac.cpp
1080%$$
1081$head Examples$$
1082The file $cref atomic_for_sparse_jac.cpp$$ contains an example and test
1083that uses this routine.
1084It returns true if the test passes and false if it fails.
1085
1086$end
1087-----------------------------------------------------------------------------
1088*/
1089/*!
1090Link from forward Jacobian sparsity sweep to atomic_base
1091
1092\param q
1093is the column dimension for the Jacobian sparsity partterns.
1094
1095\param r
1096is the Jacobian sparsity pattern for the argument vector x
1097
1098\param s
1099is the Jacobian sparsity pattern for the result vector y
1100*/
1101virtual bool for_sparse_jac(
1102        size_t                                  q  ,
1103        const vector< std::set<size_t> >&       r  ,
1104              vector< std::set<size_t> >&       s  )
1105{       return false; }
1106virtual bool for_sparse_jac(
1107        size_t                                  q  ,
1108        const vector<bool>&                     r  ,
1109              vector<bool>&                     s  )
1110{       return false; }
1111virtual bool for_sparse_jac(
1112        size_t                                  q  ,
1113        const vectorBool&                       r  ,
1114              vectorBool&                       s  )
1115{       return false; }
1116/*
1117-------------------------------------- ---------------------------------------
1118$begin atomic_rev_sparse_jac$$
1119$spell
1120        sq
1121        mul.hpp
1122        rt
1123        afun
1124        Jacobian
1125        jac
1126        CppAD
1127        std
1128        bool
1129        const
1130        hes
1131$$
1132
1133$section Atomic Reverse Jacobian Sparsity Patterns$$
1134
1135$head Syntax$$
1136$icode%ok% = %afun%.rev_sparse_jac(%q%, %rt%, %st%)%$$
1137
1138$head Purpose$$
1139This function is used by $cref RevSparseJac$$ to compute
1140Jacobian sparsity patterns.
1141For a fixed matrix $latex R \in B^{q \times m}$$,
1142the Jacobian of $latex R * f( x )$$ with respect to $latex x \in B^q$$ is
1143$latex \[
1144        S(x) = R * f^{(1)} (x)
1145\] $$
1146Given a $cref/sparsity pattern/glossary/Sparsity Pattern/$$ for $latex R$$,
1147$code rev_sparse_jac$$ computes a sparsity pattern for $latex S(x)$$.
1148
1149$head Implementation$$
1150If you are using
1151$cref RevSparseJac$$ or $cref ForSparseHes$$,
1152this virtual function must be defined by the
1153$cref/atomic_user/atomic_ctor/atomic_user/$$ class.
1154
1155$subhead q$$
1156The argument $icode q$$ has prototype
1157$codei%
1158     size_t %q%
1159%$$
1160It specifies the number of rows in
1161$latex R \in B^{q \times m}$$ and the Jacobian
1162$latex S(x) \in B^{q \times n}$$.
1163
1164$subhead rt$$
1165This argument has prototype
1166$codei%
1167     const %atomic_sparsity%& %rt%
1168%$$
1169and is a
1170$cref/atomic_sparsity/atomic_option/atomic_sparsity/$$ pattern for
1171$latex R^\R{T} \in B^{m \times q}$$.
1172
1173$subhead st$$
1174This argument has prototype
1175$codei%
1176        %atomic_sparsity%& %st%
1177%$$
1178The input value of its elements
1179are not specified (must not matter).
1180Upon return, $icode s$$ is a
1181$cref/atomic_sparsity/atomic_option/atomic_sparsity/$$ pattern for
1182$latex S(x)^\R{T} \in B^{n \times q}$$.
1183
1184$head ok$$
1185The return value $icode ok$$ has prototype
1186$codei%
1187        bool %ok%
1188%$$
1189If it is $code true$$, the corresponding evaluation succeeded,
1190otherwise it failed.
1191
1192$children%
1193        example/atomic/rev_sparse_jac.cpp
1194%$$
1195$head Examples$$
1196The file $cref atomic_rev_sparse_jac.cpp$$ contains an example and test
1197that uses this routine.
1198It returns true if the test passes and false if it fails.
1199
1200$end
1201-----------------------------------------------------------------------------
1202*/
1203/*!
1204Link from reverse Jacobian sparsity sweep to atomic_base
1205
1206\param q [in]
1207is the row dimension for the Jacobian sparsity partterns
1208
1209\param rt [out]
1210is the tansposed Jacobian sparsity pattern w.r.t to range variables y
1211
1212\param st [in]
1213is the tansposed Jacobian sparsity pattern for the argument variables x
1214*/
1215virtual bool rev_sparse_jac(
1216        size_t                                  q  ,
1217        const vector< std::set<size_t> >&       rt ,
1218              vector< std::set<size_t> >&       st )
1219{       return false; }
1220virtual bool rev_sparse_jac(
1221        size_t                                  q  ,
1222        const vector<bool>&                     rt ,
1223              vector<bool>&                     st )
1224{       return false; }
1225virtual bool rev_sparse_jac(
1226        size_t                                  q  ,
1227        const vectorBool&                       rt ,
1228              vectorBool&                       st )
1229{       return false; }
1230/*
1231-------------------------------------- ---------------------------------------
1232$begin atomic_for_sparse_hes$$
1233$spell
1234        sq
1235        mul.hpp
1236        vx
1237        afun
1238        Jacobian
1239        jac
1240        CppAD
1241        std
1242        bool
1243        hes
1244        const
1245$$
1246
1247$section Atomic Forward Hessian Sparsity Patterns$$
1248
1249$head Syntax$$
1250$icode%ok% = %afun%.for_sparse_hes(%vx%, %r%, %s%, %h%)%$$
1251
1252$head Purpose$$
1253This function is used by $cref ForSparseHes$$ to compute
1254Hessian sparsity patterns.
1255Given a $cref/sparsity pattern/glossary/Sparsity Pattern/$$ for
1256a diagonal matrix $latex R \in B^{n \times n}$$, and
1257a row vector $latex S \in B^{1 \times m}$$,
1258this routine computes the sparsity pattern for
1259$latex \[
1260        H(x) = R^\R{T} \cdot (S \cdot f)^{(2)}( x ) \cdot R
1261\] $$
1262
1263$head Implementation$$
1264If you are using and $cref ForSparseHes$$,
1265this virtual function must be defined by the
1266$cref/atomic_user/atomic_ctor/atomic_user/$$ class.
1267
1268$subhead vx$$
1269The argument $icode vx$$ has prototype
1270$codei%
1271     const CppAD:vector<bool>& %vx%
1272%$$
1273$icode%vx%.size() == %n%$$, and
1274for $latex j = 0 , \ldots , n-1$$,
1275$icode%vx%[%j%]%$$ is true if and only if
1276$icode%ax%[%j%]%$$ is a $cref/variable/glossary/Variable/$$
1277in the corresponding call to
1278$codei%
1279        %afun%(%ax%, %ay%)
1280%$$
1281
1282$subhead r$$
1283This argument has prototype
1284$codei%
1285     const CppAD:vector<bool>& %r%
1286%$$
1287and is a $cref/atomic_sparsity/atomic_option/atomic_sparsity/$$ pattern for
1288the diagonal of $latex R \in B^{n \times n}$$.
1289
1290$subhead s$$
1291The argument $icode s$$ has prototype
1292$codei%
1293     const CppAD:vector<bool>& %s%
1294%$$
1295and its size is $icode m$$.
1296It is a sparsity pattern for $latex S \in B^{1 \times m}$$.
1297
1298$subhead h$$
1299This argument has prototype
1300$codei%
1301     %atomic_sparsity%& %h%
1302%$$
1303The input value of its elements
1304are not specified (must not matter).
1305Upon return, $icode v$$ is a
1306$cref/atomic_sparsity/atomic_option/atomic_sparsity/$$ pattern for
1307$latex H(x) \in B^{n \times n}$$ which is defined above.
1308
1309$children%
1310        example/atomic/for_sparse_hes.cpp
1311%$$
1312$head Examples$$
1313The file $cref atomic_for_sparse_hes.cpp$$ contains an example and test
1314that uses this routine.
1315It returns true if the test passes and false if it fails.
1316
1317$end
1318-----------------------------------------------------------------------------
1319*/
1320/*!
1321Link from forward Hessian sparsity sweep to base_atomic
1322
1323\param vx [in]
1324which componens of x are variables.
1325
1326\param r [in]
1327is the forward Jacobian sparsity pattern w.r.t the argument vector x.
1328
1329\param s [in]
1330is the reverse Jacobian sparsity pattern w.r.t the result vector y.
1331
1332\param h [out]
1333is the Hessian sparsity pattern w.r.t the argument vector x.
1334*/
1335virtual bool for_sparse_hes(
1336        const vector<bool>&             vx ,
1337        const vector<bool>&             r  ,
1338        const vector<bool>&             s  ,
1339        vector< std::set<size_t> >&     h  )
1340{       return false; }
1341virtual bool for_sparse_hes(
1342        const vector<bool>&             vx ,
1343        const vector<bool>&             r  ,
1344        const vector<bool>&             s  ,
1345        vector<bool>&                   h  )
1346{       return false; }
1347virtual bool for_sparse_hes(
1348        const vector<bool>&             vx ,
1349        const vector<bool>&             r  ,
1350        const vector<bool>&             s  ,
1351        vectorBool&                     h  )
1352{       return false; }
1353/*
1354-------------------------------------- ---------------------------------------
1355$begin atomic_rev_sparse_hes$$
1356$spell
1357        sq
1358        mul.hpp
1359        vx
1360        afun
1361        Jacobian
1362        jac
1363        CppAD
1364        std
1365        bool
1366        hes
1367        const
1368$$
1369
1370$section Atomic Reverse Hessian Sparsity Patterns$$
1371
1372$head Syntax$$
1373$icode%ok% = %afun%.rev_sparse_hes(%vx%, %s%, %t%, %q%, %r%, %u%, %v%)%$$
1374
1375$head Purpose$$
1376This function is used by $cref RevSparseHes$$ to compute
1377Hessian sparsity patterns.
1378There is an unspecified scalar valued function
1379$latex g : B^m \rightarrow B$$.
1380Given a $cref/sparsity pattern/glossary/Sparsity Pattern/$$ for
1381$latex R \in B^{n \times q}$$,
1382and information about the function $latex z = g(y)$$,
1383this routine computes the sparsity pattern for
1384$latex \[
1385        V(x) = (g \circ f)^{(2)}( x ) R
1386\] $$
1387
1388$head Implementation$$
1389If you are using and $cref RevSparseHes$$,
1390this virtual function must be defined by the
1391$cref/atomic_user/atomic_ctor/atomic_user/$$ class.
1392
1393$subhead vx$$
1394The argument $icode vx$$ has prototype
1395$codei%
1396     const CppAD:vector<bool>& %vx%
1397%$$
1398$icode%vx%.size() == %n%$$, and
1399for $latex j = 0 , \ldots , n-1$$,
1400$icode%vx%[%j%]%$$ is true if and only if
1401$icode%ax%[%j%]%$$ is a $cref/variable/glossary/Variable/$$
1402in the corresponding call to
1403$codei%
1404        %afun%(%ax%, %ay%)
1405%$$
1406
1407$subhead s$$
1408The argument $icode s$$ has prototype
1409$codei%
1410     const CppAD:vector<bool>& %s%
1411%$$
1412and its size is $icode m$$.
1413It is a sparsity pattern for
1414$latex S(x) = g^{(1)} (y) \in B^{1 \times m}$$.
1415
1416$subhead t$$
1417This argument has prototype
1418$codei%
1419     CppAD:vector<bool>& %t%
1420%$$
1421and its size is $icode m$$.
1422The input values of its elements
1423are not specified (must not matter).
1424Upon return, $icode t$$ is a
1425sparsity pattern for
1426$latex T(x) \in B^{1 \times n}$$ where
1427$latex \[
1428        T(x) = (g \circ f)^{(1)} (x) = S(x) * f^{(1)} (x)
1429\]$$
1430
1431$subhead q$$
1432The argument $icode q$$ has prototype
1433$codei%
1434     size_t %q%
1435%$$
1436It specifies the number of columns in
1437$latex R \in B^{n \times q}$$,
1438$latex U(x) \in B^{m \times q}$$, and
1439$latex V(x) \in B^{n \times q}$$.
1440
1441$subhead r$$
1442This argument has prototype
1443$codei%
1444     const %atomic_sparsity%& %r%
1445%$$
1446and is a $cref/atomic_sparsity/atomic_option/atomic_sparsity/$$ pattern for
1447$latex R \in B^{n \times q}$$.
1448
1449$head u$$
1450This argument has prototype
1451$codei%
1452     const %atomic_sparsity%& %u%
1453%$$
1454and is a $cref/atomic_sparsity/atomic_option/atomic_sparsity/$$ pattern for
1455$latex U(x) \in B^{m \times q}$$ which is defined by
1456$latex \[
1457\begin{array}{rcl}
1458U(x)
1459& = &
1460\partial_u \{ \partial_y g[ y + f^{(1)} (x) R u ] \}_{u=0}
1461\\
1462& = &
1463\partial_u \{ g^{(1)} [ y + f^{(1)} (x) R u ] \}_{u=0}
1464\\
1465& = &
1466g^{(2)} (y) f^{(1)} (x) R
1467\end{array}
1468\] $$
1469
1470$subhead v$$
1471This argument has prototype
1472$codei%
1473     %atomic_sparsity%& %v%
1474%$$
1475The input value of its elements
1476are not specified (must not matter).
1477Upon return, $icode v$$ is a
1478$cref/atomic_sparsity/atomic_option/atomic_sparsity/$$ pattern for
1479$latex V(x) \in B^{n \times q}$$ which is defined by
1480$latex \[
1481\begin{array}{rcl}
1482V(x)
1483& = &
1484\partial_u [ \partial_x (g \circ f) ( x + R u )  ]_{u=0}
1485\\
1486& = &
1487\partial_u [ (g \circ f)^{(1)}( x + R u )  ]_{u=0}
1488\\
1489& = &
1490(g \circ f)^{(2)}( x ) R
1491\\
1492& = &
1493f^{(1)} (x)^\R{T} g^{(2)} ( y ) f^{(1)} (x)  R
1494+
1495\sum_{i=1}^m g_i^{(1)} (y) \; f_i^{(2)} (x) R
1496\\
1497& = &
1498f^{(1)} (x)^\R{T} U(x)
1499+
1500\sum_{i=1}^m S_i (x) \; f_i^{(2)} (x) R
1501\end{array}
1502\] $$
1503
1504$children%
1505        example/atomic/rev_sparse_hes.cpp
1506%$$
1507$head Examples$$
1508The file $cref atomic_rev_sparse_hes.cpp$$ contains an example and test
1509that uses this routine.
1510It returns true if the test passes and false if it fails.
1511
1512$end
1513-----------------------------------------------------------------------------
1514*/
1515/*!
1516Link from reverse Hessian sparsity sweep to base_atomic
1517
1518\param vx [in]
1519which componens of x are variables.
1520
1521\param s [in]
1522is the reverse Jacobian sparsity pattern w.r.t the result vector y.
1523
1524\param t [out]
1525is the reverse Jacobian sparsity pattern w.r.t the argument vector x.
1526
1527\param q [in]
1528is the column dimension for the sparsity partterns.
1529
1530\param r [in]
1531is the forward Jacobian sparsity pattern w.r.t the argument vector x
1532
1533\param u [in]
1534is the Hessian sparsity pattern w.r.t the result vector y.
1535
1536\param v [out]
1537is the Hessian sparsity pattern w.r.t the argument vector x.
1538*/
1539virtual bool rev_sparse_hes(
1540        const vector<bool>&                     vx ,
1541        const vector<bool>&                     s  ,
1542              vector<bool>&                     t  ,
1543        size_t                                  q  ,
1544        const vector< std::set<size_t> >&       r  ,
1545        const vector< std::set<size_t> >&       u  ,
1546              vector< std::set<size_t> >&       v  )
1547{       return false; }
1548virtual bool rev_sparse_hes(
1549        const vector<bool>&                     vx ,
1550        const vector<bool>&                     s  ,
1551              vector<bool>&                     t  ,
1552        size_t                                  q  ,
1553        const vector<bool>&                     r  ,
1554        const vector<bool>&                     u  ,
1555              vector<bool>&                     v  )
1556{       return false; }
1557virtual bool rev_sparse_hes(
1558        const vector<bool>&                     vx ,
1559        const vector<bool>&                     s  ,
1560              vector<bool>&                     t  ,
1561        size_t                                  q  ,
1562        const vectorBool&                       r  ,
1563        const vectorBool&                       u  ,
1564              vectorBool&                       v  )
1565{       return false; }
1566/*
1567------------------------------------------------------------------------------
1568$begin atomic_base_clear$$
1569$spell
1570        sq
1571        alloc
1572$$
1573
1574$section Free Static Variables$$
1575$mindex clear$$
1576
1577$head Syntax$$
1578$codei%atomic_base<%Base%>::clear()%$$
1579
1580$head Purpose$$
1581Each $code atomic_base$$ objects holds onto work space in order to
1582avoid repeated memory allocation calls and thereby increase speed
1583(until it is deleted).
1584If an the $code atomic_base$$ object is global or static because,
1585the it does not get deleted.
1586This is a problem when using
1587$code thread_alloc$$ $cref/free_all/ta_free_all/$$
1588to check that all allocated memory has been freed.
1589Calling this $code clear$$ function will free all the
1590memory currently being held onto by the
1591$codei%atomic_base<%Base%>%$$ class.
1592
1593$head Future Use$$
1594If there is future use of an $code atomic_base$$ object,
1595after a call to $code clear$$,
1596the work space will be reallocated and held onto.
1597
1598$head Restriction$$
1599This routine cannot be called
1600while in $cref/parallel/ta_in_parallel/$$ execution mode.
1601
1602$end
1603------------------------------------------------------------------------------
1604*/
1605/*!
1606Free all thread_alloc static memory held by atomic_base (avoids reallocations).
1607(This does not include class_object() which is an std::vector.)
1608*/
1609/// Free vector memory used by this class (work space)
1610static void clear(void)
1611{       CPPAD_ASSERT_KNOWN(
1612                ! thread_alloc::in_parallel() ,
1613                "cannot use atomic_base clear during parallel execution"
1614        );
1615        size_t i = class_object().size();
1616        while(i--)
1617        {       size_t thread = CPPAD_MAX_NUM_THREADS;
1618                while(thread--)
1619                {
1620                        atomic_base* op = class_object()[i];
1621                        if( op != CPPAD_NULL )
1622                        {       op->afun_vx_[thread].clear();
1623                                op->afun_vy_[thread].clear();
1624                                op->afun_tx_[thread].clear();
1625                                op->afun_ty_[thread].clear();
1626                        }
1627                }
1628        }
1629        return;
1630}
1631// -------------------------------------------------------------------------
1632/*!
1633Set value of id (used by deprecated old_atomic class)
1634
1635This function is called just before calling any of the virtual funcitons
1636and has the corresponding id of the corresponding virtual call.
1637*/
1638virtual void set_id(size_t id)
1639{ }
1640// ---------------------------------------------------------------------------
1641};
1642} // END_CPPAD_NAMESPACE
1643# endif
Note: See TracBrowser for help on using the repository browser.