Changeset 4024


Ignore:
Timestamp:
Sep 23, 2018 1:19:16 PM (9 months ago)
Author:
bradbell
Message:

merge to branch: trunk
from repository: https://github.com/coin-or/CppAD
start hash code: 46be2e2dbcaba47d8a38e8d9110efc48c7245754
end hash code: 50c63912234ddd32949460e8843399ac7f6df9b5

commit 50c63912234ddd32949460e8843399ac7f6df9b5
Author: Brad Bell <bradbell@…>
Date: Sun Sep 23 09:41:09 2018 -0700

master: add base2ad capability to checkpoint forward and reverse.

commit 423b330668e0e3f6d5eb97259ddcc5e2aac7189e
Author: Brad Bell <bradbell@…>
Date: Sun Sep 23 06:52:10 2018 -0700

master: Advance to cppad-20180923

commit 259744c0819c30bb648b0eb8ac362cdcb86f6a98
Author: Brad Bell <bradbell@…>
Date: Sun Sep 23 06:51:59 2018 -0700

master: atomic: improve see also, mention memory in atomic_base.

commit e7fe98386156854f28cc041ffc370aaa3c360ad3
Author: Brad Bell <bradbell@…>
Date: Sat Sep 22 12:24:38 2018 -0700

master: checkpoint: improve memory allocation documentation.

commit 1540c7cdac18ee3ff2b10f7ec3989c991be8dc6c
Author: Brad Bell <bradbell@…>
Date: Sat Sep 22 10:27:12 2018 -0700

master: checkpoint.hpp: remove some doxygen formattig commands.

commit e271a776a1725caf354346abc077a2e04345707f
Author: Brad Bell <bradbell@…>
Date: Sat Sep 22 09:27:39 2018 -0700

master: checkpoint: remove special case for tmb (no longer necessary).

commit 4e319c0b93d77723b36f8c77dc868dbbfb5161c5
Author: Brad Bell <bradbell@…>
Date: Sat Sep 22 09:11:18 2018 -0700

master: multi_thread/checkpoint.cpp test now passes.

commit b72d4c47b980fd92b2180aa0ead3fbb358835ef9
Author: Brad Bell <bradbell@…>
Date: Sat Sep 22 08:39:15 2018 -0700

master: working on geting multi_thread/checkpoint.cpp to pass its tests.


thread_test.cpp: print result for each number of threads.

commit 93adf79454964553c06e3a6d70e5a529d63e41c2
Author: Brad Bell <bradbell@…>
Date: Sat Sep 22 07:16:37 2018 -0700

master: Advance to cppad-20180922

commit 4bcd023ed4af9d4f52ebed766a8022af59f5574b
Author: Brad Bell <bradbell@…>
Date: Sat Sep 22 07:16:28 2018 -0700

master: checkpoint: attempt to get example/multi_thread/checkpoint.cpp to work.

commit b911efa5a513a3a414f001d05570b05f7e75732d
Author: Brad Bell <bradbell@…>
Date: Fri Sep 21 09:20:27 2018 -0700

master: multi_thread: add checkpoint example (Under Construction).
This example has been added because there is a suspected bug using checkpoint
with multiple threads.

commit 0dfb5649037c3a5c0ccef376c0b0449bde9d737f
Author: Brad Bell <bradbell@…>
Date: Fri Sep 21 05:35:13 2018 -0700

master: Advance to cppad-20180921

commit 536b06f59630a8d16f3ec85b1c01675d49855bd1
Author: Brad Bell <bradbell@…>
Date: Fri Sep 21 05:35:03 2018 -0700

master: make checkpoint/ctor.hpp a separate file.

commit d4c2647a46b926da8cf55126ffa09e036bb6ebef
Author: Brad Bell <bradbell@…>
Date: Thu Sep 20 16:23:15 2018 -0700

master: make checkpoint/reverse.hpp a separate file.

commit 459b8769a165c779e2b553b88136a15166738919
Author: Brad Bell <bradbell@…>
Date: Thu Sep 20 16:16:44 2018 -0700

master: make checkpoint/forward.hpp a separate file.

commit 55c6f34a9ef82a837e11a74c9da49ead54807031
Author: Brad Bell <bradbell@…>
Date: Thu Sep 20 13:20:24 2018 -0700

master: make checkpoint/rev_sparse_hes.hpp a separate file.

commit d15f401dec253340bb2df404c93b0a0ce6c0fd42
Author: Brad Bell <bradbell@…>
Date: Thu Sep 20 13:16:13 2018 -0700

master: checkpoint.hpp: make checkpoint/rev_sparse_jac.hpp a separate file.

commit 1b7561be8e9378784523931591ad0c69b47ef1c2
Author: Brad Bell <bradbell@…>
Date: Thu Sep 20 13:10:12 2018 -0700

master: checkpoint.hpp: make checkpoint/for_sparse_jac.hpp a separate file.

commit 8a44207afae9c2b4b527a988f6abc5354fc9bd03
Author: Brad Bell <bradbell@…>
Date: Thu Sep 20 11:49:34 2018 -0700

master: atomic_base.hpp: add atomic_base2ad.cpp example (deleted by mistake)

commit 60fdd10df33ffd30ca1b729fb3f0aceed47b2b7b
Author: Brad Bell <bradbell@…>
Date: Thu Sep 20 11:02:03 2018 -0700

master: checkpoint.hpp: make set_hes_sparse_bool.hpp a separate file.

commit 63d57df665178dd1eb93386f614228e9fc8adca3
Author: Brad Bell <bradbell@…>
Date: Thu Sep 20 10:59:56 2018 -0700

master: checkpoint.hpp: make set_hes_sparse_set.hpp a separate file.

commit 50f19625a3f6b7828a317fe819411873de88a6ae
Author: Brad Bell <bradbell@…>
Date: Thu Sep 20 10:56:54 2018 -0700

master: checkpoint.hpp: make set_jac_sparse_bool.hpp a separate file.

commit b5b029dbc88bd078b9999e92d3a58639fa77536f
Author: Brad Bell <bradbell@…>
Date: Thu Sep 20 10:23:53 2018 -0700

master: Advance to cppad-20180920

commit 2498480606fa639de5faa4c8796b4012a6828d14
Author: Brad Bell <bradbell@…>
Date: Thu Sep 20 10:19:55 2018 -0700

master: checkpoint.hpp: make set_jac_sparse_set.hpp a separate file.

commit d26d949132df354a62356988668b16123e0a0c49
Author: Brad Bell <bradbell@…>
Date: Wed Sep 19 22:05:00 2018 -0700

master: fix atomics in sparsity pattern calculations by ADFun< AD<Base>, Base > objects

commit 5c4eb81045564b2677773221f245678820d14436
Author: Brad Bell <bradbell@…>
Date: Wed Sep 19 14:57:51 2018 -0700

master: base2ad.cpp: edit comments in example.

commit ee6139b8b698e7d0751c699b2f283c7cb1b15d31
Merge: a19cdd86 48d85185
Author: Brad Bell <bradbell@…>
Date: Wed Sep 19 14:55:18 2018 -0700

Merge remote-tracking branch 'origin/master'

commit a19cdd86794a7dc3632a7d1869cc6089b724e472
Author: Brad Bell <bradbell@…>
Date: Wed Sep 19 13:21:29 2018 -0700

master: Advance to cppad-20180919

commit 3c09e3c2199cfe8999333960658dd965d6f07075
Author: Brad Bell <bradbell@…>
Date: Wed Sep 19 13:21:18 2018 -0700

master:
batch_edit.sh: add sorting to plan.
subgraph_reverse.hpp: add missing not_used_rec_base.
whats_new_18.omh: user's view of merge of base2ad.

commit a99fe1dfae5c231c5bdb7334c3d34f4e10c2d6be
Author: Brad Bell <bradbell@…>
Date: Wed Sep 19 11:18:48 2018 -0700

base2ad: extend base2ad atomics to work with reverse mode.

commit 335003af3576234710b793f389111c7e9a09a271
Author: Brad Bell <bradbell@…>
Date: Wed Sep 19 09:08:41 2018 -0700

base2ad: first commit where atomic function forward works with base2ad.

commit 48d85185a34398d35cd4152eaec2423783fb8af9
Author: Brad Bell <bradbell@…>
Date: Wed Sep 19 07:09:17 2018 -0700

master: batch_edit.sh: add to plan.

commit 4d54341db7fdf36e204d065f900582516019e1d5
Author: Brad Bell <bradbell@…>
Date: Wed Sep 19 07:07:39 2018 -0700

base2ad: cahnge result for base2ad to be ADFun< AD<Base>, Base>.

commit 22cb2a4a8566d30802d68a21f19f174675519cd5
Author: Brad Bell <bradbell@…>
Date: Tue Sep 18 10:34:06 2018 -0700

base2ad: Add RecBase? template parameter to ADFun type.


example_list.omh: add missing reference to base2ad.cpp.

commit f96708115f3c09730ce2e30082dfb0255627b476
Author: Brad Bell <bradbell@…>
Date: Mon Sep 17 13:45:26 2018 -0700

master: group user API and non user API functions.

commit 0008c201290460f47361ec8627efcb550f4f9a8b
Merge: d98dbd14 575a5c3f
Author: Brad Bell <bradbell@…>
Date: Mon Sep 17 09:50:57 2018 -0700

Merge branch 'master' into base2ad

commit 575a5c3fd26642da987bb99a374b7967dd634791
Author: Brad Bell <bradbell@…>
Date: Mon Sep 17 08:32:58 2018 -0700

master: user's view of spliting atomic functions into separate files.

commit 1559a985cb77dffbdf59d956264a8d22dd7e1940
Author: Brad Bell <bradbell@…>
Date: Mon Sep 17 07:56:07 2018 -0700

master: Advance to cppad-20180917

commit 74a019e6ff5bf7c60bb66729c2b57439ee611057
Author: Brad Bell <bradbell@…>
Date: Mon Sep 17 07:55:50 2018 -0700

master: split user API part of atomic_base.hpp into separate files.

commit d98dbd14da4993708418f8e001a554bab2e5efd3
Merge: e906f02e 47cc2b4c
Author: Brad Bell <bradbell@…>
Date: Sun Sep 16 17:49:39 2018 -0700

Merge branch 'master' into base2ad

commit 47cc2b4cd8777573b2b42b92a5bdf7962fa05b7c
Author: Brad Bell <bradbell@…>
Date: Sun Sep 16 17:48:49 2018 -0700

master: remove con_ad_type (corresponds to tape_id does not match).

commit e906f02ea898ffda314749ecfb2bcf4152a1da01
Author: Brad Bell <bradbell@…>
Date: Sun Sep 16 07:56:12 2018 -0700

base2ad: change ode_taylor -> taylor_ode.

commit 706cda6e746c150c196352526945e3a46f050c48
Author: Brad Bell <bradbell@…>
Date: Sun Sep 16 07:53:43 2018 -0700

base2ad: make mul_level_adolc_ode.cpp more like base2ad.cpp.

commit c438d447499d2d76b6676155f3c93b86a965dc06
Author: Brad Bell <bradbell@…>
Date: Sun Sep 16 07:33:00 2018 -0700

base2ad: mul_level_ode.cpp: use typedefs and make more like base2ad.cpp.

commit 3f0b37e04057b9a790654c123ec887d805456a39
Merge: 7d7f695e cc02600d
Author: Brad Bell <bradbell@…>
Date: Sun Sep 16 06:46:17 2018 -0700

Merge branch 'master' into base2ad

commit cc02600d325bfd72348a194c45a87ee982ed1320
Author: Brad Bell <bradbell@…>
Date: Sun Sep 16 06:40:43 2018 -0700

master: improve ode_taylor.cpp and move it to taylor_ode.cpp.

commit 97bfaac023564cf34a265c80ddc6b2e4f9dac2c5
Author: Brad Bell <bradbell@…>
Date: Sun Sep 16 06:31:10 2018 -0700

master: change ode_taylor.cpp to use same notation as taylor_ode.omh.

commit 7d7f695e3cb41efb7958632fe5cfd6699cdb2acb
Author: Brad Bell <bradbell@…>
Date: Sun Sep 16 06:00:35 2018 -0700

base2ad: put base2ad in separate file and document it.

commit ac3fdbc24d4836537aa2f845a9ec26b188388991
Author: Brad Bell <bradbell@…>
Date: Sun Sep 16 05:58:45 2018 -0700

base2ad: first working version of base2ad.

commit 78014e8ef632e9fc804ebea74422acf4c64f0c69
Author: Brad Bell <bradbell@…>
Date: Sun Sep 16 04:34:57 2018 -0700

master: Advance to cppad-20180916

commit 0cd2a999f9566655fda723482d2fd060b6867e8e
Author: Brad Bell <bradbell@…>
Date: Sun Sep 16 04:34:47 2018 -0700

master: example general.cpp: include base_require during testing.

commit fdc230eba87cf19312cddd92d5dd66fe2434aa60
Author: Brad Bell <bradbell@…>
Date: Sat Sep 15 15:34:33 2018 -0700

master: move AD theory for Taylor's ODE method to separate section.

commit c02959c71c4ba48ceb54d3caa644c4d4d3e3e320
Author: Brad Bell <bradbell@…>
Date: Sat Sep 15 14:28:45 2018 -0700

master: batch_edit.sh: add to plan.

commit 3e72aad050219ad4993fccc3e2cc94f0abd7c496
Author: Brad Bell <bradbell@…>
Date: Sat Sep 15 08:51:41 2018 -0700

master: fun_construct.hpp: use sparse_pack and sparse_list assignment operators.

commit eb5906dfd75e5fa557ac307dadd0b860c33733ed
Author: Brad Bell <bradbell@…>
Date: Sat Sep 15 07:31:57 2018 -0700

master: improve move semantic version of ADFun assignment operator.

commit b183932db18501467b5ea2e238ccfa737d76020a
Author: Brad Bell <bradbell@…>
Date: Sat Sep 15 06:54:02 2018 -0700

master: add move semantic version of function assignment.

commit 38f8a5a553d2dbbbe7ed8158077548cc80dfe643
Author: Brad Bell <bradbell@…>
Date: Sat Sep 15 06:05:47 2018 -0700

master: ADFun assignment operator: remove resizes that are no longer necessary.

commit 2d4a7245e06a6dc5e3e83ad3e09726340201821e
Author: Brad Bell <bradbell@…>
Date: Sat Sep 15 05:33:18 2018 -0700

master: remove Erase from recorder and player (no longer used).

commit 75736474de93a6064329a9f578ae4f8383b0fa78
Author: Brad Bell <bradbell@…>
Date: Sat Sep 15 05:21:18 2018 -0700

master: ad_fun.hpp: convert dep_parameter_ from vector -> pod_vector.

commit 5ed186c4ecef2958913325694100e5384f03d8ec
Author: Brad Bell <bradbell@…>
Date: Sat Sep 15 05:12:25 2018 -0700

master: Advance to cppad-20180915

commit b99ba06f766d64233b1c0ae66e590d3a189467a6
Author: Brad Bell <bradbell@…>
Date: Sat Sep 15 05:12:14 2018 -0700

master: ad_fun.hpp: change dep_taddr_ and ind_taddr_ from vector -> pod_vector.

commit 7872c7f69a6618176ae76dd2bf64a0fa71a65d14
Author: Brad Bell <bradbell@…>
Date: Fri Sep 14 08:32:22 2018 -0700

master: ad_fun.hpp: add more comments about location of doxygen documentation.

commit 73c0d52cb060ca7bbdc0d2fe6a956f0936fdeed2
Author: Brad Bell <bradbell@…>
Date: Fri Sep 14 07:40:24 2018 -0700

master: Advance to cppad-20180914

commit d24872e359b53a1863bd43af4f3a61ec77c149bf
Author: Brad Bell <bradbell@…>
Date: Fri Sep 14 07:40:14 2018 -0700

master: ad_fun.hpp: use full path when referencing files where doxygen is.

commit bf73b3fc3fa88abab2fc1913515b5525c62d5290
Author: Brad Bell <bradbell@…>
Date: Thu Sep 13 07:54:50 2018 -0700

master: Advance to cppad-20180913

commit 5fb67b72e237e00e021e4e40f08dc4ba6d53c224
Author: Brad Bell <bradbell@…>
Date: Thu Sep 13 07:54:29 2018 -0700

master: wish_list.omh: add AD<Base> Functions from Base Functions item.

commit 690ed2af5feb2d4a36267acf3669d0d02aef6158
Author: Brad Bell <bradbell@…>
Date: Mon Aug 27 12:19:14 2018 -0700

master: move eigen_vector from global namspace to CppAD namespace.

commit 3ac03943e02e60832b1f8a3311a96cf45b810c21
Author: Brad Bell <bradbell@…>
Date: Mon Aug 27 06:29:29 2018 -0700

master: Advance to cppad-20180827

commit a74a90ca37c3c290bbe4c71fb6fffccc5f8b44e3
Author: Brad Bell <bradbell@…>
Date: Mon Aug 27 06:29:17 2018 -0700

master: dynamic parameters: improve discussion and example.

Location:
trunk
Files:
30 added
1 deleted
96 edited

Legend:

Unmodified
Added
Removed
  • trunk/CMakeLists.txt

    r4023 r4024  
    2525#
    2626# cppad_version is used by version.sh to get the version number.
    27 SET(cppad_version "20180824")
     27SET(cppad_version "20180923")
    2828SET(cppad_url          "http://www.coin-or.org/CppAD" )
    2929SET(cppad_description  "Differentiation of C++ Algorithms" )
  • trunk/bin/batch_edit.sh

    r4023 r4024  
    2727# 14. Remove all inlines for functions that depend on template parameters.
    2828# 15. Change 'ind[0-9] -> arg[0-9] when used as operator arguments.
     29# 16. Remove all lines that have 'SHORT COPYRIGHT' on them.
     30# 17. Remove 'It returns true if it succeeds and false otherwise.'
     31# 18. Change template <class *> -> template <typename *>.
     32# 19. Create check_sort.sh and use it to sort all alphabetical lists.
    2933# -----------------------------------------------------------------------------
    3034spell_list='
     
    3438move_list='
    3539'
    36 move_sed='s|par_var_dyn|con_dyn_var|'
     40move_sed='s|ode_taylor.cpp|taylor_ode.cpp|'
    3741#
    3842cat << EOF > junk.sed
    39 s/if( ! IdenticalZero( *right.value_ *) )/if( dyn_right | (! IdenticalZero(right.value_) ) )/
    40 #
    41 s|if( IdenticalZero( *right.value_ *) )|if( (! dyn_right) \\& IdenticalZero(right.value_) )|
    42 s|if( IdenticalOne( *right.value_ *) )|if( (! dyn_right) \\& IdenticalOne(right.value_) )|
    43 #
    44 s|if( IdenticalZero( *left *) )|if( (! dyn_left) \\& IdenticalZero(left) )|
    45 s|if( IdenticalZero( *left.value_ *) )|if( (! dyn_left) \\& IdenticalZero(left.value_) )|
    46 s|if( IdenticalOne( *left *) )|if( (! dyn_left) \\& IdenticalOne(left) )|
    47 s|if( IdenticalOne( *left.value_ *) )|if( (! dyn_left) \\& IdenticalOne(left.value_) )|
     43/template *<[a-z]* *Base>/! b end
     44N
     45N
     46/ADFun<Base>::/! b end
     47s|template *<class Base>|template <class Base, class RecBase>|
     48s|template <typename Base>|template <typename Base, typename RecBase>|
     49s|ADFun<Base>::|ADFun<Base,RecBase>::|g
     50: end
    4851EOF
    4952# -----------------------------------------------------------------------------
  • trunk/bin/check_all.sh

    r4023 r4024  
    246246        echo_log_eval ./$program multi_atomic 1 4 100
    247247        #
     248        # test_time=1,max_thread=4,num_solve=100
     249        next_program
     250        echo_log_eval ./$program checkpoint 1 4 100
     251        #
    248252        # test_time=2,max_thread=4,num_zero=20,num_sub=30,num_sum=50,use_ad=true
    249253        next_program
  • trunk/configure

    r4023 r4024  
    11#! /bin/sh
    22# Guess values for system-dependent variables and create Makefiles.
    3 # Generated by GNU Autoconf 2.69 for cppad 20180824.
     3# Generated by GNU Autoconf 2.69 for cppad 20180923.
    44#
    55# Report bugs to <cppad@list.coin-or.org>.
     
    580580PACKAGE_NAME='cppad'
    581581PACKAGE_TARNAME='cppad'
    582 PACKAGE_VERSION='20180824'
    583 PACKAGE_STRING='cppad 20180824'
     582PACKAGE_VERSION='20180923'
     583PACKAGE_STRING='cppad 20180923'
    584584PACKAGE_BUGREPORT='cppad@list.coin-or.org'
    585585PACKAGE_URL=''
     
    13751375  # This message is too long to be a string in the A/UX 3.1 sh.
    13761376  cat <<_ACEOF
    1377 \`configure' configures cppad 20180824 to adapt to many kinds of systems.
     1377\`configure' configures cppad 20180923 to adapt to many kinds of systems.
    13781378
    13791379Usage: $0 [OPTION]... [VAR=VALUE]...
     
    14451445if test -n "$ac_init_help"; then
    14461446  case $ac_init_help in
    1447      short | recursive ) echo "Configuration of cppad 20180824:";;
     1447     short | recursive ) echo "Configuration of cppad 20180923:";;
    14481448   esac
    14491449  cat <<\_ACEOF
     
    15791579if $ac_init_version; then
    15801580  cat <<\_ACEOF
    1581 cppad configure 20180824
     1581cppad configure 20180923
    15821582generated by GNU Autoconf 2.69
    15831583
     
    19521952running configure, to aid debugging if configure makes a mistake.
    19531953
    1954 It was created by cppad $as_me 20180824, which was
     1954It was created by cppad $as_me 20180923, which was
    19551955generated by GNU Autoconf 2.69.  Invocation command line was
    19561956
     
    28422842# Define the identity of the package.
    28432843 PACKAGE='cppad'
    2844  VERSION='20180824'
     2844 VERSION='20180923'
    28452845
    28462846
     
    78697869# values after options handling.
    78707870ac_log="
    7871 This file was extended by cppad $as_me 20180824, which was
     7871This file was extended by cppad $as_me 20180923, which was
    78727872generated by GNU Autoconf 2.69.  Invocation command line was
    78737873
     
    79267926ac_cs_config="`$as_echo "$ac_configure_args" | sed 's/^ //; s/[\\""\`\$]/\\\\&/g'`"
    79277927ac_cs_version="\\
    7928 cppad config.status 20180824
     7928cppad config.status 20180923
    79297929configured by $0, generated by GNU Autoconf 2.69,
    79307930  with options \\"\$ac_cs_config\\"
  • trunk/configure.ac

    r4023 r4024  
    1212dnl Process this file with autoconf to produce a configure script.
    1313dnl   package   version              bug-report
    14 AC_INIT([cppad], [20180824], [cppad@list.coin-or.org])
     14AC_INIT([cppad], [20180923], [cppad@list.coin-or.org])
    1515AM_SILENT_RULES([yes])
    1616
  • trunk/cppad/core/abs_normal_fun.hpp

    r4023 r4024  
    293293# endif
    294294
    295 template <class Base>
    296 void ADFun<Base>::abs_normal_fun(ADFun<Base>& g, ADFun<Base>& a) const
     295template <class Base, class RecBase>
     296void ADFun<Base,RecBase>::abs_normal_fun(ADFun& g, ADFun& a) const
    297297{       using namespace local;
    298298
  • trunk/cppad/core/ad_fun.hpp

    r4022 r4024  
    6666
    6767
    68 template <class Base>
     68template <typename Base, typename RecBase>
    6969class ADFun {
    70 // ------------------------------------------------------------
    71 // Private member variables
     70        // ADFun<Base> must be a friend of ADFun< AD<Base> > for base2ad to work.
     71        template <typename Base2, typename RecBase2> friend class ADFun;
    7272private:
     73        // ------------------------------------------------------------
     74        // Private member variables
     75        // ------------------------------------------------------------
    7376        /// Has this ADFun object been optmized
    7477        bool has_been_optimized_;
     
    104107
    105108        /// tape address for the independent variables
    106         CppAD::vector<size_t> ind_taddr_;
     109        local::pod_vector<size_t> ind_taddr_;
    107110
    108111        /// tape address and parameter flag for the dependent variables
    109         CppAD::vector<size_t> dep_taddr_;
     112        local::pod_vector<size_t> dep_taddr_;
    110113
    111114        /// which dependent variables are actually parameters
    112         CppAD::vector<bool>  dep_parameter_;
     115        local::pod_vector<bool> dep_parameter_;
    113116
    114117        /// results of the forward mode calculations
     
    143146        local::pod_vector_maybe<Base> subgraph_partial_;
    144147
    145 // ------------------------------------------------------------
    146 // Private member functions
     148        // ------------------------------------------------------------
     149        // Private member functions
     150        // ------------------------------------------------------------
    147151
    148152        /// change the operation sequence corresponding to this object
     
    150154        void Dependent(local::ADTape<Base> *tape, const ADvector &y);
    151155
    152         // ------------------------------------------------------------
    153156        // vector of bool version of ForSparseJac
    154         // (see doxygen in for_sparse_jac.hpp)
     157        // (doxygen in cppad/core/for_sparse_jac.hpp)
    155158        template <class VectorSet>
    156159        void ForSparseJacCase(
     
    162165                VectorSet&         s
    163166        );
     167
    164168        // vector of std::set<size_t> version of ForSparseJac
    165         // (see doxygen in for_sparse_jac.hpp)
     169        // (doxygen in cppad/core/for_sparse_jac.hpp)
    166170        template <class VectorSet>
    167171        void ForSparseJacCase(
     
    173177                VectorSet&               s
    174178        );
    175         // ------------------------------------------------------------
     179
    176180        // vector of bool version of RevSparseJac
    177         // (see doxygen in rev_sparse_jac.hpp)
     181        // (doxygen in cppad/core/rev_sparse_jac.hpp)
    178182        template <class VectorSet>
    179183        void RevSparseJacCase(
     
    185189                VectorSet&         r
    186190        );
     191
    187192        // vector of std::set<size_t> version of RevSparseJac
    188         // (see doxygen in rev_sparse_jac.hpp)
     193        // (doxygen in cppad/core/rev_sparse_jac.hpp)
    189194        template <class VectorSet>
    190195        void RevSparseJacCase(
     
    196201                VectorSet&               r
    197202        );
    198         // ------------------------------------------------------------
     203
    199204        // vector of bool version of ForSparseHes
    200         // (see doxygen in for_sparse_hes.hpp)
     205        // (doxygen in cppad/core/for_sparse_hes.hpp)
    201206        template <class VectorSet>
    202207        void ForSparseHesCase(
     
    206211                VectorSet&         h
    207212        );
     213
    208214        // vector of std::set<size_t> version of ForSparseHes
    209         // (see doxygen in for_sparse_hes.hpp)
     215        // (doxygen in cppad/core/for_sparse_hes.hpp)
    210216        template <class VectorSet>
    211217        void ForSparseHesCase(
     
    215221                VectorSet&               h
    216222        );
    217         // ------------------------------------------------------------
     223
    218224        // vector of bool version of RevSparseHes
    219         // (see doxygen in rev_sparse_hes.hpp)
     225        // (doxygen in cppad/core/rev_sparse_hes.hpp)
    220226        template <class VectorSet>
    221227        void RevSparseHesCase(
     
    226232                VectorSet&         h
    227233        );
     234
    228235        // vector of std::set<size_t> version of RevSparseHes
    229         // (see doxygen in rev_sparse_hes.hpp)
     236        // (doxygen in cppad/core/rev_sparse_hes.hpp)
    230237        template <class VectorSet>
    231238        void RevSparseHesCase(
     
    236243                VectorSet&               h
    237244        );
    238         // ------------------------------------------------------------
     245
    239246        // Forward mode version of SparseJacobian
    240         // (see doxygen in sparse_jacobian.hpp)
     247        // (doxygen in cppad/core/sparse_jacobian.hpp)
    241248        template <class VectorBase, class VectorSet, class VectorSize>
    242249        size_t SparseJacobianFor(
     
    248255                      sparse_jacobian_work& work
    249256        );
     257
    250258        // Reverse mode version of SparseJacobian
    251         // (see doxygen in sparse_jacobian.hpp)
     259        // (doxygen in cppad/core/sparse_jacobian.hpp)
    252260        template <class VectorBase, class VectorSet, class VectorSize>
    253261        size_t SparseJacobianRev(
     
    259267                      sparse_jacobian_work& work
    260268        );
    261         // ------------------------------------------------------------
    262         // combined sparse_list and sparse_pack version of
    263         // SparseHessian (see doxygen in sparse_hessian.hpp)
     269
     270        // combined sparse_list and sparse_pack version of SparseHessian
     271        // (doxygen in cppad/core/sparse_hessian.hpp)
    264272        template <class VectorBase, class VectorSet, class VectorSize>
    265273        size_t SparseHessianCompute(
     
    272280                      sparse_hessian_work&     work
    273281        );
    274 // ------------------------------------------------------------
     282
    275283public:
    276284        /// copy constructor
     
    292300
    293301        // assignment operator
    294         // (see doxygen in fun_construct.hpp)
     302        // (doxygen in cppad/core/fun_construct.hpp)
    295303        void operator=(const ADFun& f);
     304# if CPPAD_USE_CPLUSPLUS_2011
     305        // assignment operator with move semantics
     306        void operator=(ADFun&& f);
     307# endif
     308
     309        // create ADFun< AD<Base> > from this ADFun<Base>
     310        // (doxygen in cppad/core/base2ad.hpp)
     311        ADFun< AD<Base>, RecBase > base2ad(void) const;
    296312
    297313        /// sequence constructor
     
    331347        VectorBase Reverse(size_t p, const VectorBase &v);
    332348
    333         // ---------------------------------------------------------------------
    334         // Jacobian sparsity
     349        // forward Jacobian sparsity pattern
     350        // (doxygen in cppad/core/for_sparse_jac.hpp)
    335351        template <typename VectorSet>
    336352        VectorSet ForSparseJac(
     
    338354                bool dependency = false
    339355        );
     356
     357        // reverse Jacobian sparsity pattern
     358        // (doxygen in cppad/core/rev_sparse_jac.hpp)
    340359        template <typename VectorSet>
    341360        VectorSet RevSparseJac(
     
    343362                bool dependency = false
    344363        );
    345         // ---------------------------------------------------------------------
     364
     365        // subgraph_reverse: select domain
     366        // (doxygen in cppad/core/subgraph_reverse.hpp)
    346367        template <typename VectorBool>
    347368        void subgraph_reverse(
    348369                const VectorBool&                   select_domain
    349370        );
     371
     372        // subgraph_reverse: compute derivative
     373        // (doxygen in cppad/core/subgraph_reverse.hpp)
    350374        template <typename Addr, typename VectorBase, typename SizeVector>
    351375        void subgraph_reverse_helper(
     
    355379                VectorBase&                          dw
    356380        );
     381
     382        // subgraph_reverse: compute derivative
     383        // (doxygen in cppad/core/subgraph_reverse.hpp)
    357384        template <typename VectorBase, typename SizeVector>
    358385        void subgraph_reverse(
     
    362389                VectorBase&                          dw
    363390        );
     391
     392        // subgraph_jac_rev: compute Jacobian
     393        // (doxygen in cppad/core/subgraph_jac_rev.hpp)
    364394        template <typename SizeVector, typename BaseVector>
    365395        void subgraph_jac_rev(
     
    367397                sparse_rcv<SizeVector, BaseVector>&  subset
    368398        );
     399
     400        // subgraph_jac_rev: compute Jacobian
     401        // (doxygen missing in cppad/core/subgraph_jac_rev.hpp)
    369402        template <typename BoolVector, typename SizeVector, typename BaseVector>
    370403        void subgraph_jac_rev(
     
    374407                sparse_rcv<SizeVector, BaseVector>&  matrix_out
    375408        );
     409
     410
     411        // compute sparse Jacobian using forward mode
     412        // (doxygen in cppad/core/sparse_jac.hpp)
    376413        template <typename SizeVector, typename BaseVector>
    377414        size_t sparse_jac_for(
     
    383420                sparse_jac_work&                     work
    384421        );
     422
     423        // compute sparse Jacobian using reverse mode
     424        // (doxygen in cppad/core/sparse_jac.hpp)
    385425        template <typename SizeVector, typename BaseVector>
    386426        size_t sparse_jac_rev(
     
    391431                sparse_jac_work&                     work
    392432        );
     433
     434        // compute sparse Hessian
     435        // (doxygen in cppad/core/sparse_hes.hpp)
    393436        template <typename SizeVector, typename BaseVector>
    394437        size_t sparse_hes(
     
    400443                sparse_hes_work&                     work
    401444        );
    402         // ---------------------------------------------------------------------
     445
     446        // compute sparsity pattern using subgraphs
     447        // (doxygen in cppad/core/subgraph_sparsity.hpp)
    403448        template <typename BoolVector, typename SizeVector>
    404449        void subgraph_sparsity(
     
    408453                sparse_rc<SizeVector>&       pattern_out
    409454        );
     455
     456
     457        // forward mode Jacobian sparsity pattern
     458        // (doxygen in cppad/core/for_jac_sparsity.hpp)
    410459        template <typename SizeVector>
    411460        void for_jac_sparsity(
     
    416465                sparse_rc<SizeVector>&       pattern_out
    417466        );
     467
     468        // reverse mode Jacobian sparsity pattern
     469        // (doxygen in cppad/core/for_jac_sparsity.hpp)
    418470        template <typename SizeVector>
    419471        void rev_jac_sparsity(
     
    424476                sparse_rc<SizeVector>&       pattern_out
    425477        );
     478
     479        // reverse mode Hessian sparsity pattern
     480        // (doxygen in cppad/core/rev_hes_sparsity.hpp)
    426481        template <typename BoolVector, typename SizeVector>
    427482        void rev_hes_sparsity(
     
    431486                sparse_rc<SizeVector>&       pattern_out
    432487        );
     488
     489        // forward mode Hessian sparsity pattern
     490        // (doxygen in cppad/core/for_hes_sparsity.hpp)
    433491        template <typename BoolVector, typename SizeVector>
    434492        void for_hes_sparsity(
     
    438496                sparse_rc<SizeVector>&       pattern_out
    439497        );
    440         // ---------------------------------------------------------------------
    441         // forward mode Hessian sparsity
    442         // (see doxygen documentation in for_sparse_hes.hpp)
     498
     499        // forward mode Hessian sparsity pattern
     500        // (see doxygen in cppad/core/for_sparse_hes.hpp)
    443501        template <typename VectorSet>
    444502        VectorSet ForSparseHes(
    445503                const VectorSet &r, const VectorSet &s
    446504        );
     505
    447506        // internal set sparsity version of ForSparseHes
    448507        // (used by checkpoint functions only)
     
    452511                local::sparse_list&                  h
    453512        );
    454         // reverse mode Hessian sparsity
    455         // (see doxygen documentation in rev_sparse_hes.hpp)
     513
     514        // reverse mode Hessian sparsity pattern
     515        // (see doxygen in cppad/core/rev_sparse_hes.hpp)
    456516        template <typename VectorSet>
    457517        VectorSet RevSparseHes(
    458518                size_t q, const VectorSet &s, bool transpose = false
    459519        );
     520
    460521        // internal set sparsity version of RevSparseHes
     522        // (doxygen in cppad/core/rev_sparse_hes.hpp)
    461523        // (used by checkpoint functions only)
    462524        void RevSparseHesCheckpoint(
     
    466528                local::sparse_list&                  h
    467529        );
     530
    468531        // internal set sparsity version of RevSparseJac
     532        // (doxygen in cppad/core/rev_sparse_jac.hpp)
    469533        // (used by checkpoint functions only)
    470534        void RevSparseJacCheckpoint(
     
    475539                local::sparse_list&                  s
    476540        );
    477     // internal set sparsity version of RevSparseJac
    478     // (used by checkpoint functions only)
     541
     542        // internal set sparsity version of RevSparseJac
     543        // (doxygen in cppad/core/for_sparse_jac.hpp)
     544        // (used by checkpoint functions only)
    479545        void ForSparseJacCheckpoint(
    480546        size_t                        q          ,
     
    776842# include <cppad/core/dependent.hpp>
    777843# include <cppad/core/fun_construct.hpp>
     844# include <cppad/core/base2ad.hpp>
    778845# include <cppad/core/abort_recording.hpp>
    779846# include <cppad/core/fun_eval.hpp>
  • trunk/cppad/core/ad_valued.hpp

    r3845 r4024  
    1 // $Id$
    21# ifndef CPPAD_CORE_AD_VALUED_HPP
    32# define CPPAD_CORE_AD_VALUED_HPP
    43
    54/* --------------------------------------------------------------------------
    6 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-16 Bradley M. Bell
     5CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-18 Bradley M. Bell
    76
    87CppAD is distributed under multiple licenses. This distribution is under
     
    2221$section AD Valued Operations and Functions$$
    2322
    24 $comment atomic.omh includes atomic_base.omh which atomic_base.hpp$$
     23$comment atomic.omh includes atomic_base.hpp$$
    2524$childtable%
    2625        cppad/core/arithmetic.hpp%
  • trunk/cppad/core/atomic_base.hpp

    r4023 r4024  
    1212Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
    1313-------------------------------------------------------------------------- */
     14/*
     15$begin atomic_base$$
     16$spell
     17        ctor
     18        cppad
     19        hpp
     20        afun
     21        arg
     22        vx
     23        vy
     24        tx
     25        ty
     26        px
     27        py
     28        jac
     29        hes
     30        CppAD
     31        checkpointing
     32        algo
     33$$
     34
     35$section User Defined Atomic AD Functions$$
     36
     37
     38$head Syntax$$
     39
     40$codei%
     41%atomic_user% %afun%(%ctor_arg_list%)
     42%afun%(%ax%, %ay%)
     43%ok% = %afun%.forward(%p%, %q%, %vx%, %vy%, %tx%, %ty%)
     44%ok% = %afun%.reverse(%q%, %tx%, %ty%, %px%, %py%)
     45%ok% = %afun%.for_sparse_jac(%q%, %r%, %s%, %x%)
     46%ok% = %afun%.rev_sparse_jac(%q%, %r%, %s%, %x%)
     47%ok% = %afun%.for_sparse_hes(%vx%, %r%, %s%, %h%, %x%)
     48%ok% = %afun%.rev_sparse_hes(%vx%, %s%, %t%, %q%, %r%, %u%, %v%, %x%)
     49atomic_base<%Base%>::clear()%$$
     50
     51$head See Also$$
     52$cref checkpoint$$
     53
     54$head Purpose$$
     55
     56$subhead Speed$$
     57In some cases, the user knows how to compute derivatives of a function
     58$latex \[
     59        y = f(x) \; {\rm where} \; f : B^n \rightarrow B^m
     60\] $$
     61more efficiently than by coding it using $codei%AD<%Base%>%$$
     62$cref/atomic/glossary/Operation/Atomic/$$ operations
     63and letting CppAD do the rest.
     64In this case $codei%atomic_base%<%Base%>%$$ can use
     65the user code for $latex f(x)$$, and its derivatives,
     66as $codei%AD<%Base%>%$$ atomic operations.
     67
     68$subhead Reduce Memory$$
     69If the function $latex f(x)$$ is used often,
     70using an atomic version of $latex f(x)$$ remove the need for repeated
     71copies of the corresponding $codei%AD<%Base%>%$$ operations.
     72
     73$head Virtual Functions$$
     74User defined derivatives are implemented by defining the
     75following virtual functions in the $icode base_atomic$$ class:
     76$cref/forward/atomic_forward/$$,
     77$cref/reverse/atomic_reverse/$$,
     78$cref/for_sparse_jac/atomic_for_sparse_jac/$$,
     79$cref/rev_sparse_jac/atomic_rev_sparse_jac/$$, and
     80$cref/rev_sparse_hes/atomic_rev_sparse_hes/$$.
     81These virtual functions have a default implementation
     82that returns $icode%ok% == false%$$.
     83The $code forward$$ function,
     84for the case $icode%q% == 0%$$, must be implemented.
     85Otherwise, only those functions
     86required by the your calculations need to be implemented.
     87For example,
     88$icode forward$$ for the case $icode%q% == 2%$$ can just return
     89$icode%ok% == false%$$ unless you require
     90forward mode calculation of second derivatives.
     91
     92$head Examples$$
     93See $cref atomic_example$$.
     94
     95$childtable%
     96        cppad/core/atomic_base/ctor.hpp%
     97        cppad/core/atomic_base/option.hpp%
     98        cppad/core/atomic_base/afun.hpp%
     99        cppad/core/atomic_base/forward.hpp%
     100        cppad/core/atomic_base/reverse.hpp%
     101        cppad/core/atomic_base/for_sparse_jac.hpp%
     102        cppad/core/atomic_base/rev_sparse_jac.hpp%
     103        cppad/core/atomic_base/for_sparse_hes.hpp%
     104        cppad/core/atomic_base/rev_sparse_hes.hpp%
     105        cppad/core/atomic_base/clear.hpp
     106%$$
     107
     108$end
     109-------------------------------------------------------------------------------
     110$begin atomic_example$$
     111
     112$section Example User Defined Atomic Functions$$
     113
     114$head Getting Started$$
     115The file $cref atomic_get_started.cpp$$ contains an example and test
     116that shows the minimal amount of information required to create
     117a user defined atomic operation.
     118
     119$head Scalar Function$$
     120The file $cref atomic_reciprocal.cpp$$ contains an example and test
     121where the user provides the code for computing derivatives.
     122This example is simple because the domain and range are scalars.
     123
     124$head Vector Range$$
     125The file $cref atomic_tangent.cpp$$ contains another example
     126where the user provides the code for computing derivatives.
     127This example is more complex because the range has two components.
     128
     129$head Hessian Sparsity Patterns$$
     130The file $cref atomic_rev_sparse_hes.cpp$$ contains an minimal example
     131where the user provides the code for computing Hessian sparsity patterns.
     132
     133$head General Case$$
     134The file $cref atomic_mat_mul.cpp$$ contains a more general example
     135where the user provides the code for computing derivatives.
     136This example is more complex because both the domain and range
     137dimensions are arbitrary.
     138
     139$childtable%
     140        example/atomic/get_started.cpp%
     141        example/atomic/norm_sq.cpp%
     142        example/atomic/reciprocal.cpp%
     143        example/atomic/set_sparsity.cpp%
     144        example/atomic/tangent.cpp%
     145        example/atomic/eigen_mat_mul.cpp%
     146        example/atomic/eigen_mat_inv.cpp%
     147        example/atomic/eigen_cholesky.cpp%
     148        example/atomic/mat_mul.cpp%
     149        example/atomic/base2ad.cpp
     150%$$
     151
     152$end
     153-------------------------------------------------------------------------------
     154*/
    14155
    15156# include <set>
    16157# include <cppad/core/cppad_assert.hpp>
    17158# include <cppad/local/sparse_internal.hpp>
     159
    18160// needed before one can use CPPAD_ASSERT_FIRST_CALL_NOT_PARALLEL
    19161# include <cppad/utility/thread_alloc.hpp>
     
    53195                vector<bool>               vx;
    54196                vector<bool>               vy;
     197                //
    55198                vector<Base>               tx;
    56199                vector<Base>               ty;
     200                //
     201                vector< AD<Base> >         atx;
     202                vector< AD<Base> >         aty;
    57203                //
    58204                vector<bool>               bool_t;
     
    92238                return list_;
    93239        }
     240public:
    94241        // =====================================================================
    95 public:
    96         // -----------------------------------------------------
    97         // member functions not in user API
     242        // In User API
     243        // =====================================================================
    98244        //
     245        // ---------------------------------------------------------------------
     246        // ctor: doxygen in atomic_base/ctor.hpp
     247        atomic_base(void);
     248        atomic_base(
     249                const std::string&     name,
     250                option_enum            sparsity = bool_sparsity_enum
     251        );
     252
     253        // option: see doxygen in atomic_base/option.hpp
     254        void option(enum option_enum option_value);
     255
     256        // operator(): see doxygen in atomic_base/afun.hpp
     257        template <class ADVector>
     258        void operator()(
     259                const ADVector&  ax     ,
     260                      ADVector&  ay     ,
     261                size_t           id = 0
     262        );
     263
     264        // ------------------------------------------------------------------------
     265        // forward: see docygen in atomic_base/forward.hpp
     266        virtual bool forward(
     267                size_t                    p  ,
     268                size_t                    q  ,
     269                const vector<bool>&       vx ,
     270                      vector<bool>&       vy ,
     271                const vector<Base>&       tx ,
     272                      vector<Base>&       ty
     273        );
     274        virtual bool forward(
     275                size_t                    p  ,
     276                size_t                    q  ,
     277                const vector<bool>&       vx ,
     278                      vector<bool>&       vy ,
     279                const vector< AD<Base> >& atx ,
     280                      vector< AD<Base> >& aty
     281        );
     282
     283        // ------------------------------------------------------------------------
     284        // reverse: see docygen in atomic_base/reverse.hpp
     285        virtual bool reverse(
     286                size_t                    q  ,
     287                const vector<Base>&       tx ,
     288                const vector<Base>&       ty ,
     289                      vector<Base>&       px ,
     290                const vector<Base>&       py
     291        );
     292        virtual bool reverse(
     293                size_t                    q   ,
     294                const vector< AD<Base> >& atx ,
     295                const vector< AD<Base> >& aty ,
     296                      vector< AD<Base> >& apx ,
     297                const vector< AD<Base> >& apy
     298        );
     299
     300        // ------------------------------------------------------------
     301        // for_sparse_jac: see doxygen in atomic_base/for_sparse_jac.hpp
     302        virtual bool for_sparse_jac(
     303                size_t                                  q  ,
     304                const vector< std::set<size_t> >&       r  ,
     305                      vector< std::set<size_t> >&       s  ,
     306                const vector<Base>&                     x
     307        );
     308        virtual bool for_sparse_jac(
     309                size_t                                  q  ,
     310                const vector<bool>&                     r  ,
     311                      vector<bool>&                     s  ,
     312                const vector<Base>&                     x
     313        );
     314        virtual bool for_sparse_jac(
     315                size_t                                  q  ,
     316                const vectorBool&                       r  ,
     317                      vectorBool&                       s  ,
     318                const vector<Base>&                     x
     319        );
     320        template <class InternalSparsity>
     321        void for_sparse_jac(
     322                const vector<Base>&              x            ,
     323                const local::pod_vector<size_t>& x_index      ,
     324                const local::pod_vector<size_t>& y_index      ,
     325                InternalSparsity&                var_sparsity
     326        );
     327        // deprecated versions
     328        virtual bool for_sparse_jac(
     329                size_t                                  q  ,
     330                const vector< std::set<size_t> >&       r  ,
     331                      vector< std::set<size_t> >&       s
     332        );
     333        virtual bool for_sparse_jac(
     334                size_t                                  q  ,
     335                const vector<bool>&                     r  ,
     336                      vector<bool>&                     s
     337        );
     338        virtual bool for_sparse_jac(
     339                size_t                                  q  ,
     340                const vectorBool&                       r  ,
     341                      vectorBool&                       s
     342        );
     343        // ------------------------------------------------------------
     344        // rev_sparse_jac: see doxygen in atomic_base/rev_sparse_jac.hpp
     345        virtual bool rev_sparse_jac(
     346                size_t                                  q  ,
     347                const vector< std::set<size_t> >&       rt ,
     348                      vector< std::set<size_t> >&       st ,
     349                const vector<Base>&                     x
     350        );
     351        virtual bool rev_sparse_jac(
     352                size_t                                  q  ,
     353                const vector<bool>&                     rt ,
     354                      vector<bool>&                     st ,
     355                const vector<Base>&                     x
     356        );
     357        virtual bool rev_sparse_jac(
     358                size_t                                  q  ,
     359                const vectorBool&                       rt ,
     360                      vectorBool&                       st ,
     361                const vector<Base>&                     x
     362        );
     363        template <class InternalSparsity>
     364        void rev_sparse_jac(
     365                const vector<Base>&        x            ,
     366                const local::pod_vector<size_t>& x_index ,
     367                const local::pod_vector<size_t>& y_index ,
     368                InternalSparsity&          var_sparsity
     369        );
     370        // deprecated versions
     371        virtual bool rev_sparse_jac(
     372                size_t                                  q  ,
     373                const vector< std::set<size_t> >&       rt ,
     374                      vector< std::set<size_t> >&       st
     375        );
     376        virtual bool rev_sparse_jac(
     377                size_t                                  q  ,
     378                const vector<bool>&                     rt ,
     379                      vector<bool>&                     st
     380        );
     381        virtual bool rev_sparse_jac(
     382                size_t                                  q  ,
     383                const vectorBool&                       rt ,
     384                      vectorBool&                       st
     385        );
     386        // ------------------------------------------------------------
     387        // for_sparse_hes: see doxygen in atomic_base/for_sparse_hes.hpp
     388        virtual bool for_sparse_hes(
     389                const vector<bool>&             vx ,
     390                const vector<bool>&             r  ,
     391                const vector<bool>&             s  ,
     392                vector< std::set<size_t> >&     h  ,
     393                const vector<Base>&             x
     394        );
     395        virtual bool for_sparse_hes(
     396                const vector<bool>&             vx ,
     397                const vector<bool>&             r  ,
     398                const vector<bool>&             s  ,
     399                vector<bool>&                   h  ,
     400                const vector<Base>&             x
     401        );
     402        virtual bool for_sparse_hes(
     403                const vector<bool>&             vx ,
     404                const vector<bool>&             r  ,
     405                const vector<bool>&             s  ,
     406                vectorBool&                     h  ,
     407                const vector<Base>&             x
     408        );
     409        template <class InternalSparsity>
     410        void for_sparse_hes(
     411                const vector<Base>&              x                ,
     412                const local::pod_vector<size_t>& x_index          ,
     413                const local::pod_vector<size_t>& y_index          ,
     414                const InternalSparsity&          for_jac_sparsity ,
     415                const InternalSparsity&          rev_jac_sparsity ,
     416                InternalSparsity&                for_hes_sparsity
     417        );
     418        // deprecated versions
     419        virtual bool for_sparse_hes(
     420                const vector<bool>&             vx ,
     421                const vector<bool>&             r  ,
     422                const vector<bool>&             s  ,
     423                vector< std::set<size_t> >&     h
     424        );
     425        virtual bool for_sparse_hes(
     426                const vector<bool>&             vx ,
     427                const vector<bool>&             r  ,
     428                const vector<bool>&             s  ,
     429                vector<bool>&                   h
     430        );
     431        virtual bool for_sparse_hes(
     432                const vector<bool>&             vx ,
     433                const vector<bool>&             r  ,
     434                const vector<bool>&             s  ,
     435                vectorBool&                     h
     436        );
     437        // ------------------------------------------------------------
     438        // rev_sparse_hes: see doxygen in atomic_base/rev_sparse_hes.hpp
     439        virtual bool rev_sparse_hes(
     440                const vector<bool>&                     vx ,
     441                const vector<bool>&                     s  ,
     442                      vector<bool>&                     t  ,
     443                size_t                                  q  ,
     444                const vector< std::set<size_t> >&       r  ,
     445                const vector< std::set<size_t> >&       u  ,
     446                      vector< std::set<size_t> >&       v  ,
     447                const vector<Base>&                     x
     448        );
     449        virtual bool rev_sparse_hes(
     450                const vector<bool>&                     vx ,
     451                const vector<bool>&                     s  ,
     452                      vector<bool>&                     t  ,
     453                size_t                                  q  ,
     454                const vector<bool>&                     r  ,
     455                const vector<bool>&                     u  ,
     456                      vector<bool>&                     v  ,
     457                const vector<Base>&                     x
     458        );
     459        virtual bool rev_sparse_hes(
     460                const vector<bool>&                     vx ,
     461                const vector<bool>&                     s  ,
     462                      vector<bool>&                     t  ,
     463                size_t                                  q  ,
     464                const vectorBool&                       r  ,
     465                const vectorBool&                       u  ,
     466                      vectorBool&                       v  ,
     467                const vector<Base>&                     x
     468        );
     469        template <class InternalSparsity>
     470        void rev_sparse_hes(
     471                const vector<Base>&              x                ,
     472                const local::pod_vector<size_t>& x_index          ,
     473                const local::pod_vector<size_t>& y_index          ,
     474                const InternalSparsity&          for_jac_sparsity ,
     475                bool*                            rev_jac_flag     ,
     476                InternalSparsity&                rev_hes_sparsity
     477        );
     478        // deprecated
     479        virtual bool rev_sparse_hes(
     480                const vector<bool>&                     vx ,
     481                const vector<bool>&                     s  ,
     482                      vector<bool>&                     t  ,
     483                size_t                                  q  ,
     484                const vector< std::set<size_t> >&       r  ,
     485                const vector< std::set<size_t> >&       u  ,
     486                      vector< std::set<size_t> >&       v
     487        );
     488        virtual bool rev_sparse_hes(
     489                const vector<bool>&                     vx ,
     490                const vector<bool>&                     s  ,
     491                      vector<bool>&                     t  ,
     492                size_t                                  q  ,
     493                const vector<bool>&                     r  ,
     494                const vector<bool>&                     u  ,
     495                      vector<bool>&                     v
     496        );
     497        virtual bool rev_sparse_hes(
     498                const vector<bool>&                     vx ,
     499                const vector<bool>&                     s  ,
     500                      vector<bool>&                     t  ,
     501                size_t                                  q  ,
     502                const vectorBool&                       r  ,
     503                const vectorBool&                       u  ,
     504                      vectorBool&                       v
     505        );
     506        // ------------------------------------------------------------
     507        // clear: see doxygen in atomic_base/clear.hpp
     508        static void clear(void);
     509
     510        // =====================================================================
     511        // Not in User API
     512        // =====================================================================
     513
    99514        /// current sparsity setting
    100515        option_enum sparsity(void) const
     
    104519        const std::string& afun_name(void) const
    105520        {       return class_name()[index_]; }
    106 /*
    107 $begin atomic_ctor$$
    108 $spell
    109         enum
    110         sq
    111         std
    112         afun
    113         arg
    114         CppAD
    115         bool
    116         ctor
    117         const
    118         mat_mul_xam.cpp
    119         hpp
    120 $$
    121 
    122 $section Atomic Function Constructor$$
    123 
    124 $head Syntax$$
    125 $icode%atomic_user afun%(%ctor_arg_list%)
    126 %$$
    127 $codei%atomic_base<%Base%>(%name%, %sparsity%)
    128 %$$
    129 
    130 $head atomic_user$$
    131 
    132 $subhead ctor_arg_list$$
    133 Is a list of arguments for the $icode atomic_user$$ constructor.
    134 
    135 $subhead afun$$
    136 The object $icode afun$$ must stay in scope for as long
    137 as the corresponding atomic function is used.
    138 This includes use by any $cref/ADFun<Base>/ADFun/$$ that
    139 has this $icode atomic_user$$ operation in its
    140 $cref/operation sequence/glossary/Operation/Sequence/$$.
    141 
    142 $subhead Implementation$$
    143 The user defined $icode atomic_user$$ class is a publicly derived class of
    144 $codei%atomic_base<%Base%>%$$.
    145 It should be declared as follows:
    146 $codei%
    147         class %atomic_user% : public CppAD::atomic_base<%Base%> {
    148         public:
    149                 %atomic_user%(%ctor_arg_list%) : atomic_base<%Base%>(%name%, %sparsity%)
    150         %...%
    151         };
    152 %$$
    153 where $icode ...$$
    154 denotes the rest of the implementation of the derived class.
    155 This includes completing the constructor and
    156 all the virtual functions that have their
    157 $code atomic_base$$ implementations replaced by
    158 $icode atomic_user$$ implementations.
    159 
    160 $head atomic_base$$
    161 
    162 $subhead Restrictions$$
    163 The $code atomic_base$$ constructor cannot be called in
    164 $cref/parallel/ta_in_parallel/$$ mode.
    165 
    166 $subhead Base$$
    167 The template parameter determines the
    168 $icode Base$$ type for this $codei%AD<%Base%>%$$ atomic operation.
    169 
    170 $subhead name$$
    171 This $code atomic_base$$ constructor argument has the following prototype
    172 $codei%
    173         const std::string& %name%
    174 %$$
    175 It is the name for this atomic function and is used for error reporting.
    176 The suggested value for $icode name$$ is $icode afun$$ or $icode atomic_user$$,
    177 i.e., the name of the corresponding atomic object or class.
    178 
    179 $subhead sparsity$$
    180 This $code atomic_base$$ constructor argument has prototype
    181 $codei%
    182         atomic_base<%Base%>::option_enum %sparsity%
    183 %$$
    184 The current $icode sparsity$$ for an $code atomic_base$$ object
    185 determines which type of sparsity patterns it uses
    186 and its value is one of the following:
    187 $table
    188 $icode sparsity$$   $cnext sparsity patterns $rnext
    189 $codei%atomic_base<%Base%>::pack_sparsity_enum%$$ $pre  $$ $cnext
    190         $cref/vectorBool/CppAD_vector/vectorBool/$$
    191 $rnext
    192 $codei%atomic_base<%Base%>::bool_sparsity_enum%$$ $pre  $$ $cnext
    193         $cref/vector/CppAD_vector/$$$code <bool>$$
    194 $rnext
    195 $codei%atomic_base<%Base%>::set_sparsity_enum%$$ $pre  $$ $cnext
    196         $cref/vector/CppAD_vector/$$$code <std::set<std::size_t> >$$
    197 $tend
    198 There is a default value for $icode sparsity$$ if it is not
    199 included in the constructor (which may be either the bool or set option).
    200 
    201 $head Example$$
    202 
    203 $subhead Define Constructor$$
    204 The following is an example of a user atomic function constructor definitions:
    205 $cref%get_started.cpp%atomic_get_started.cpp%Constructor%$$.
    206 
    207 $subhead Use Constructor$$
    208 The following is an example using a user atomic function constructor:
    209 $cref%get_started.cpp%atomic_get_started.cpp%Use Atomic Function%Constructor%$$.
    210 
    211 $end
    212 */
    213 /*!
    214 Base class for atomic_user functions.
    215 
    216 \tparam Base
    217 This class is used for defining an AD<Base> atomic operation y = f(x).
    218 */
    219 /// make sure user does not invoke the default constructor
    220 atomic_base(void)
    221 {       CPPAD_ASSERT_KNOWN(false,
    222                 "Attempt to use the atomic_base default constructor"
    223         );
    224 }
    225 /*!
    226 Constructor
    227 
    228 \param name
    229 name used for error reporting
    230 
    231 \param sparsity [in]
    232 what type of sparsity patterns are computed by this function,
    233 bool_sparsity_enum or set_sparsity_enum. Default value is
    234 bool sparsity patterns.
    235 */
    236 atomic_base(
    237                 const std::string&     name,
    238                 option_enum            sparsity = bool_sparsity_enum
    239 ) :
    240 index_   ( class_object().size()  )  ,
    241 sparsity_( sparsity               )
    242 {       CPPAD_ASSERT_KNOWN(
    243                 ! thread_alloc::in_parallel() ,
    244                 "atomic_base: constructor cannot be called in parallel mode."
    245         );
    246         class_object().push_back(this);
    247         class_name().push_back(name);
    248         CPPAD_ASSERT_UNKNOWN( class_object().size() == class_name().size() );
    249         //
    250         // initialize work pointers as null;
    251         for(size_t thread = 0; thread < CPPAD_MAX_NUM_THREADS; thread++)
    252                 work_[thread] = CPPAD_NULL;
    253 }
    254 /// destructor informs CppAD that this atomic function with this index
    255 /// has dropped out of scope by setting its pointer to null
    256 virtual ~atomic_base(void)
    257 {       CPPAD_ASSERT_UNKNOWN( class_object().size() > index_ );
    258         // change object pointer to null, but leave name for error reporting
    259         class_object()[index_] = CPPAD_NULL;
    260         //
    261         // free temporary work memory
    262         for(size_t thread = 0; thread < CPPAD_MAX_NUM_THREADS; thread++)
    263                 free_work(thread);
    264 }
    265 /// allocates work_ for a specified thread
    266 void allocate_work(size_t thread)
    267 {       if( work_[thread] == CPPAD_NULL )
    268         {       // allocate the raw memory
    269                 size_t min_bytes = sizeof(work_struct);
    270                 size_t num_bytes;
    271                 void*  v_ptr     = thread_alloc::get_memory(min_bytes, num_bytes);
    272                 // save in work_
    273                 work_[thread]    = reinterpret_cast<work_struct*>( v_ptr );
    274                 // call constructor
    275                 new( work_[thread] ) work_struct;
    276         }
    277         return;
    278 }
    279 /// frees work_ for a specified thread
    280 void free_work(size_t thread)
    281 {       if( work_[thread] != CPPAD_NULL )
    282         {       // call destructor
    283                 work_[thread]->~work_struct();
    284                 // return memory to avialable pool for this thread
    285         thread_alloc::return_memory( reinterpret_cast<void*>(work_[thread]) );
    286                 // mark this thread as not allocated
    287                 work_[thread] = CPPAD_NULL;
    288         }
    289         return;
    290 }
    291 /// atomic_base function object corresponding to a certain index
    292 static atomic_base* class_object(size_t index)
    293 {       CPPAD_ASSERT_UNKNOWN( class_object().size() > index );
    294         return class_object()[index];
    295 }
    296 /// atomic_base function name corresponding to a certain index
    297 static const std::string& class_name(size_t index)
    298 {       CPPAD_ASSERT_UNKNOWN( class_name().size() > index );
    299         return class_name()[index];
    300 }
    301 /*
    302 $begin atomic_option$$
    303 $spell
    304         sq
    305         enum
    306         afun
    307         bool
    308         CppAD
    309         std
    310         typedef
    311 $$
    312 
    313 $section Set Atomic Function Options$$
    314 
    315 $head Syntax$$
    316 $icode%afun%.option(%option_value%)%$$
    317 These settings do not apply to individual $icode afun$$ calls,
    318 but rather all subsequent uses of the corresponding atomic operation
    319 in an $cref ADFun$$ object.
    320 
    321 $head atomic_sparsity$$
    322 Note that, if you use $cref optimize$$, these sparsity patterns are used
    323 to determine the $cref/dependency/dependency.cpp/$$ relationship between
    324 argument and result variables.
    325 
    326 $subhead pack_sparsity_enum$$
    327 If $icode option_value$$ is $codei%atomic_base<%Base%>::pack_sparsity_enum%$$,
    328 then the type used by $icode afun$$ for
    329 $cref/sparsity patterns/glossary/Sparsity Pattern/$$,
    330 (after the option is set) will be
    331 $codei%
    332         typedef CppAD::vectorBool %atomic_sparsity%
    333 %$$
    334 If $icode r$$ is a sparsity pattern
    335 for a matrix $latex R \in B^{p \times q}$$:
    336 $icode%r%.size() == %p% * %q%$$.
    337 
    338 $subhead bool_sparsity_enum$$
    339 If $icode option_value$$ is $codei%atomic_base<%Base%>::bool_sparsity_enum%$$,
    340 then the type used by $icode afun$$ for
    341 $cref/sparsity patterns/glossary/Sparsity Pattern/$$,
    342 (after the option is set) will be
    343 $codei%
    344         typedef CppAD::vector<bool> %atomic_sparsity%
    345 %$$
    346 If $icode r$$ is a sparsity pattern
    347 for a matrix $latex R \in B^{p \times q}$$:
    348 $icode%r%.size() == %p% * %q%$$.
    349 
    350 $subhead set_sparsity_enum$$
    351 If $icode option_value$$ is $icode%atomic_base<%Base%>::set_sparsity_enum%$$,
    352 then the type used by $icode afun$$ for
    353 $cref/sparsity patterns/glossary/Sparsity Pattern/$$,
    354 (after the option is set) will be
    355 $codei%
    356         typedef CppAD::vector< std::set<size_t> > %atomic_sparsity%
    357 %$$
    358 If $icode r$$ is a sparsity pattern
    359 for a matrix $latex R \in B^{p \times q}$$:
    360 $icode%r%.size() == %p%$$, and for $latex i = 0 , \ldots , p-1$$,
    361 the elements of $icode%r%[%i%]%$$ are between zero and $latex q-1$$ inclusive.
    362 
    363 $end
    364 */
    365 void option(enum option_enum option_value)
    366 {       switch( option_value )
    367         {       case pack_sparsity_enum:
    368                 case bool_sparsity_enum:
    369                 case set_sparsity_enum:
    370                 sparsity_ = option_value;
    371                 break;
    372 
    373                 default:
    374                 CPPAD_ASSERT_KNOWN(
    375                         false,
    376                         "atoic_base::option: option_value is not valid"
    377                 );
    378         }
    379         return;
    380 }
    381 /*
    382 -----------------------------------------------------------------------------
    383 $begin atomic_afun$$
    384 
    385 $spell
    386         sq
    387         mul
    388         afun
    389         const
    390         CppAD
    391         mat_mul.cpp
    392 $$
    393 
    394 $section Using AD Version of Atomic Function$$
    395 
    396 $head Syntax$$
    397 $icode%afun%(%ax%, %ay%)%$$
    398 
    399 $head Purpose$$
    400 Given $icode ax$$,
    401 this call computes the corresponding value of $icode ay$$.
    402 If $codei%AD<%Base%>%$$ operations are being recorded,
    403 it enters the computation as an atomic operation in the recording;
    404 see $cref/start recording/Independent/Start Recording/$$.
    405 
    406 $head ADVector$$
    407 The type $icode ADVector$$ must be a
    408 $cref/simple vector class/SimpleVector/$$ with elements of type
    409 $codei%AD<%Base%>%$$; see $cref/Base/atomic_ctor/atomic_base/Base/$$.
    410 
    411 $head afun$$
    412 is a $cref/atomic_user/atomic_ctor/atomic_user/$$ object
    413 and this $icode afun$$ function call is implemented by the
    414 $cref/atomic_base/atomic_ctor/atomic_base/$$ class.
    415 
    416 $head ax$$
    417 This argument has prototype
    418 $codei%
    419         const %ADVector%& %ax%
    420 %$$
    421 and size must be equal to $icode n$$.
    422 It specifies vector $latex x \in B^n$$
    423 at which an $codei%AD<%Base%>%$$ version of
    424 $latex y = f(x)$$ is to be evaluated; see
    425 $cref/Base/atomic_ctor/atomic_base/Base/$$.
    426 
    427 $head ay$$
    428 This argument has prototype
    429 $codei%
    430         %ADVector%& %ay%
    431 %$$
    432 and size must be equal to $icode m$$.
    433 The input values of its elements
    434 are not specified (must not matter).
    435 Upon return, it is an $codei%AD<%Base%>%$$ version of
    436 $latex y = f(x)$$.
    437 
    438 $head Examples$$
    439 The following files contain example uses of
    440 the AD version of atomic functions during recording:
    441 $cref%get_started.cpp%atomic_get_started.cpp%Use Atomic Function%Recording%$$,
    442 $cref%norm_sq.cpp%atomic_norm_sq.cpp%Use Atomic Function%Recording%$$,
    443 $cref%reciprocal.cpp%atomic_reciprocal.cpp%Use Atomic Function%Recording%$$,
    444 $cref%tangent.cpp%atomic_tangent.cpp%Use Atomic Function%Recording%$$,
    445 $cref%mat_mul.cpp%atomic_mat_mul.cpp%Use Atomic Function%Recording%$$.
    446 
    447 $end
    448 -----------------------------------------------------------------------------
    449 */
    450 /*!
    451 Implement the user call to <tt>afun(ax, ay)</tt> and old_atomic call to
    452 <tt>afun(ax, ay, id)</tt>.
    453 
    454 \tparam ADVector
    455 A simple vector class with elements of type <code>AD<Base></code>.
    456 
    457 \param id
    458 optional extra information vector that is just passed through by CppAD,
    459 and used by old_atomic derived class (not other derived classes).
    460 This is an extra parameter to the virtual callbacks for old_atomic;
    461 see the set_old member function.
    462 
    463 \param ax
    464 is the argument vector for this call,
    465 <tt>ax.size()</tt> determines the number of arguments.
    466 
    467 \param ay
    468 is the result vector for this call,
    469 <tt>ay.size()</tt> determines the number of results.
    470 */
    471 template <class ADVector>
    472 void operator()(
    473         const ADVector&  ax     ,
    474               ADVector&  ay     ,
    475         size_t           id = 0 )
    476 {       size_t i, j;
    477         size_t n = ax.size();
    478         size_t m = ay.size();
    479 # ifndef NDEBUG
    480         bool ok;
    481         std::string msg = "atomic_base: " + afun_name() + ".eval: ";
    482         if( (n == 0) | (m == 0) )
    483         {       msg += "ax.size() or ay.size() is zero";
    484                 CPPAD_ASSERT_KNOWN(false, msg.c_str() );
    485         }
    486 # endif
    487         size_t thread = thread_alloc::thread_num();
    488         allocate_work(thread);
    489         vector <Base>& tx  = work_[thread]->tx;
    490         vector <Base>& ty  = work_[thread]->ty;
    491         vector <bool>& vx  = work_[thread]->vx;
    492         vector <bool>& vy  = work_[thread]->vy;
    493         //
    494         if( vx.size() != n )
    495         {       vx.resize(n);
    496                 tx.resize(n);
    497         }
    498         if( vy.size() != m )
    499         {       vy.resize(m);
    500                 ty.resize(m);
    501         }
    502         //
    503         // Determine tape corresponding to variables in ax
    504         tape_id_t            tape_id  = 0;
    505         local::ADTape<Base>* tape     = CPPAD_NULL;
    506         for(j = 0; j < n; j++)
    507         {       tx[j]  = ax[j].value_;
    508                 vx[j]  = ! Constant( ax[j] );
    509                 if( vx[j] )
    510                 {
    511                         if( tape_id == 0 )
    512                         {       tape    = ax[j].tape_this();
    513                                 tape_id = ax[j].tape_id_;
    514                                 CPPAD_ASSERT_UNKNOWN( tape != CPPAD_NULL );
    515                         }
    516 # ifndef NDEBUG
    517                         if( tape_id != ax[j].tape_id_ )
    518                         {       msg += afun_name() +
    519                                 ": ax contains variables from different threads.";
    520                                 CPPAD_ASSERT_KNOWN(false, msg.c_str());
    521                         }
    522 # endif
     521
     522        /// destructor informs CppAD that this atomic function with this index
     523        /// has dropped out of scope by setting its pointer to null
     524        virtual ~atomic_base(void)
     525        {       CPPAD_ASSERT_UNKNOWN( class_object().size() > index_ );
     526                // change object pointer to null, but leave name for error reporting
     527                class_object()[index_] = CPPAD_NULL;
     528                //
     529                // free temporary work memory
     530                for(size_t thread = 0; thread < CPPAD_MAX_NUM_THREADS; thread++)
     531                        free_work(thread);
     532        }
     533        /// allocates work_ for a specified thread
     534        void allocate_work(size_t thread)
     535        {       if( work_[thread] == CPPAD_NULL )
     536                {       // allocate the raw memory
     537                        size_t min_bytes = sizeof(work_struct);
     538                        size_t num_bytes;
     539                        void*  v_ptr     = thread_alloc::get_memory(min_bytes, num_bytes);
     540                        // save in work_
     541                        work_[thread]    = reinterpret_cast<work_struct*>( v_ptr );
     542                        // call constructor
     543                        new( work_[thread] ) work_struct;
    523544                }
    524         }
    525         // Use zero order forward mode to compute values
    526         size_t p = 0, q = 0;
    527         set_old(id);
    528 # ifdef NDEBUG
    529         forward(p, q, vx, vy, tx, ty);
    530 # else
    531         ok = forward(p, q, vx, vy, tx, ty);
    532         if( ! ok )
    533         {       msg += afun_name() + ": ok is false for "
    534                         "zero order forward mode calculation.";
    535                 CPPAD_ASSERT_KNOWN(false, msg.c_str());
    536         }
    537 # endif
    538         bool record_operation = false;
    539         for(i = 0; i < m; i++)
    540         {
    541                 // pass back values
    542                 ay[i].value_ = ty[i];
    543 
    544                 // initialize entire vector parameters (not in tape)
    545                 ay[i].tape_id_ = 0;
    546                 ay[i].taddr_   = 0;
    547 
    548                 // we need to record this operation if
    549                 // any of the elemnts of ay are variables,
    550                 record_operation |= vy[i];
    551         }
    552 # ifndef NDEBUG
    553         if( record_operation & (tape == CPPAD_NULL) )
    554         {       msg +=
    555                 "all elements of vx are false but vy contains a true element";
    556                 CPPAD_ASSERT_KNOWN(false, msg.c_str() );
    557         }
    558 # endif
    559         // if tape is not null, ay is on the tape
    560         if( record_operation )
    561         {
    562                 // Operator that marks beginning of this atomic operation
    563                 CPPAD_ASSERT_UNKNOWN( local::NumRes(local::UserOp) == 0 );
    564                 CPPAD_ASSERT_UNKNOWN( local::NumArg(local::UserOp) == 4 );
    565                 CPPAD_ASSERT_KNOWN(
    566                         size_t( std::numeric_limits<addr_t>::max() ) >=
    567                         std::max( std::max( std::max(index_, id), n), m ),
    568                         "atomic_base: cppad_tape_addr_type maximum not large enough"
    569                 );
    570                 tape->Rec_.PutArg(addr_t(index_), addr_t(id), addr_t(n), addr_t(m));
    571                 tape->Rec_.PutOp(local::UserOp);
    572 
    573                 // Now put n operators, one for each element of argument vector
    574                 CPPAD_ASSERT_UNKNOWN( local::NumRes(local::UsravOp) == 0 );
    575                 CPPAD_ASSERT_UNKNOWN( local::NumRes(local::UsrapOp) == 0 );
    576                 CPPAD_ASSERT_UNKNOWN( local::NumArg(local::UsravOp) == 1 );
    577                 CPPAD_ASSERT_UNKNOWN( local::NumArg(local::UsrapOp) == 1 );
    578                 for(j = 0; j < n; j++)
    579                 {       if( Variable(ax[j]) )
    580                         {       // information for an argument that is a variable
    581                                 tape->Rec_.PutArg(ax[j].taddr_);
    582                                 tape->Rec_.PutOp(local::UsravOp);
    583                         }
    584                         else
    585                         {       // information for an argument that is parameter
    586                                 addr_t par = ax[j].taddr_;
    587                                 if( ! Dynamic( ax[j] ) )
    588                                         par = tape->Rec_.put_con_par(ax[j].value_);
    589                                 tape->Rec_.PutArg(par);
    590                                 tape->Rec_.PutOp(local::UsrapOp);
    591                         }
     545                return;
     546        }
     547        /// frees work_ for a specified thread
     548        void free_work(size_t thread)
     549        {       if( work_[thread] != CPPAD_NULL )
     550                {       // call destructor
     551                        work_[thread]->~work_struct();
     552                        // return memory to avialable pool for this thread
     553                thread_alloc::return_memory(
     554                                reinterpret_cast<void*>(work_[thread])
     555                        );
     556                        // mark this thread as not allocated
     557                        work_[thread] = CPPAD_NULL;
    592558                }
    593 
    594                 // Now put m operators, one for each element of result vector
    595                 CPPAD_ASSERT_UNKNOWN( local::NumArg(local::UsrrpOp) == 1 );
    596                 CPPAD_ASSERT_UNKNOWN( local::NumRes(local::UsrrpOp) == 0 );
    597                 CPPAD_ASSERT_UNKNOWN( local::NumArg(local::UsrrvOp) == 0 );
    598                 CPPAD_ASSERT_UNKNOWN( local::NumRes(local::UsrrvOp) == 1 );
    599                 for(i = 0; i < m; i++)
    600                 {       if( vy[i] )
    601                         {       ay[i].taddr_    = tape->Rec_.PutOp(local::UsrrvOp);
    602                                 ay[i].tape_id_  = tape_id;
    603                                 ay[i].ad_type_  = local::var_ad_type;
    604                         }
    605                         else
    606                         {       CPPAD_ASSERT_UNKNOWN( ! Dynamic( ay[i] ) );
    607                                 addr_t par = tape->Rec_.put_con_par(ay[i].value_);
    608                                 tape->Rec_.PutArg(par);
    609                                 tape->Rec_.PutOp(local::UsrrpOp);
    610                         }
    611                 }
    612 
    613                 // Put a duplicate UserOp at end of UserOp sequence
    614                 CPPAD_ASSERT_KNOWN(
    615                         size_t( std::numeric_limits<addr_t>::max() ) >=
    616                         std::max( std::max( std::max(index_, id), n), m ),
    617                         "atomic_base: cppad_tape_addr_type maximum not large enough"
    618                 );
    619                 tape->Rec_.PutArg(addr_t(index_), addr_t(id), addr_t(n), addr_t(m));
    620                 tape->Rec_.PutOp(local::UserOp);
    621         }
    622         return;
    623 }
    624 /*
    625 -----------------------------------------------------------------------------
    626 $begin atomic_forward$$
    627 $spell
    628         sq
    629         mul.hpp
    630         hes
    631         afun
    632         vx
    633         vy
    634         ty
    635         Taylor
    636         const
    637         CppAD
    638         bool
    639 $$
    640 
    641 $section Atomic Forward Mode$$
    642 $mindex callback virtual$$
    643 
    644 
    645 $head Syntax$$
    646 $icode%ok% = %afun%.forward(%p%, %q%, %vx%, %vy%, %tx%, %ty%)%$$
    647 
    648 $head Purpose$$
    649 This virtual function is used by $cref atomic_afun$$
    650 to evaluate function values.
    651 It is also used buy $cref/forward/Forward/$$
    652 to compute function vales and derivatives.
    653 
    654 $head Implementation$$
    655 This virtual function must be defined by the
    656 $cref/atomic_user/atomic_ctor/atomic_user/$$ class.
    657 It can just return $icode%ok% == false%$$
    658 (and not compute anything) for values
    659 of $icode%q% > 0%$$ that are greater than those used by your
    660 $cref/forward/Forward/$$ mode calculations.
    661 
    662 $head p$$
    663 The argument $icode p$$ has prototype
    664 $codei%
    665         size_t %p%
    666 %$$
    667 It specifies the lowest order Taylor coefficient that we are evaluating.
    668 During calls to $cref atomic_afun$$, $icode%p% == 0%$$.
    669 
    670 $head q$$
    671 The argument $icode q$$ has prototype
    672 $codei%
    673         size_t %q%
    674 %$$
    675 It specifies the highest order Taylor coefficient that we are evaluating.
    676 During calls to $cref atomic_afun$$, $icode%q% == 0%$$.
    677 
    678 $head vx$$
    679 The $code forward$$ argument $icode vx$$ has prototype
    680 $codei%
    681         const CppAD::vector<bool>& %vx%
    682 %$$
    683 The case $icode%vx%.size() > 0%$$ only occurs while evaluating a call to
    684 $cref atomic_afun$$.
    685 In this case,
    686 $icode%p% == %q% == 0%$$,
    687 $icode%vx%.size() == %n%$$, and
    688 for $latex j = 0 , \ldots , n-1$$,
    689 $icode%vx%[%j%]%$$ is true if and only if
    690 $icode%ax%[%j%]%$$ is a $cref/variable/glossary/Variable/$$
    691 or $cref/dynamic parameter/glossary/Parameter/Dynamic/$$
    692 in the corresponding call to
    693 $codei%
    694         %afun%(%ax%, %ay%)
    695 %$$
    696 If $icode%vx%.size() == 0%$$,
    697 then $icode%vy%.size() == 0%$$ and neither of these vectors
    698 should be used.
    699 
    700 $head vy$$
    701 The $code forward$$ argument $icode vy$$ has prototype
    702 $codei%
    703         CppAD::vector<bool>& %vy%
    704 %$$
    705 If $icode%vy%.size() == 0%$$, it should not be used.
    706 Otherwise,
    707 $icode%q% == 0%$$ and $icode%vy%.size() == %m%$$.
    708 The input values of the elements of $icode vy$$
    709 are not specified (must not matter).
    710 Upon return, for $latex j = 0 , \ldots , m-1$$,
    711 $icode%vy%[%i%]%$$ is true if and only if
    712 $icode%ay%[%i%]%$$ is a variable
    713 or dynamic parameter
    714 (CppAD uses $icode vy$$ to reduce the necessary computations).
    715 
    716 $head tx$$
    717 The argument $icode tx$$ has prototype
    718 $codei%
    719         const CppAD::vector<%Base%>& %tx%
    720 %$$
    721 and $icode%tx%.size() == (%q%+1)*%n%$$.
    722 For $latex j = 0 , \ldots , n-1$$ and $latex k = 0 , \ldots , q$$,
    723 we use the Taylor coefficient notation
    724 $latex \[
    725 \begin{array}{rcl}
    726         x_j^k    & = & tx [ j * ( q + 1 ) + k ]
    727         \\
    728         X_j (t)  & = & x_j^0 + x_j^1 t^1 + \cdots + x_j^q t^q
    729 \end{array}
    730 \] $$
    731 Note that superscripts represent an index for $latex x_j^k$$
    732 and an exponent for $latex t^k$$.
    733 Also note that the Taylor coefficients for $latex X(t)$$ correspond
    734 to the derivatives of $latex X(t)$$ at $latex t = 0$$ in the following way:
    735 $latex \[
    736         x_j^k = \frac{1}{ k ! } X_j^{(k)} (0)
    737 \] $$
    738 
    739 $head ty$$
    740 The argument $icode ty$$ has prototype
    741 $codei%
    742         CppAD::vector<%Base%>& %ty%
    743 %$$
    744 and $icode%tx%.size() == (%q%+1)*%m%$$.
    745 Upon return,
    746 For $latex i = 0 , \ldots , m-1$$ and $latex k = 0 , \ldots , q$$,
    747 $latex \[
    748 \begin{array}{rcl}
    749         Y_i (t)  & = & f_i [ X(t) ]
    750         \\
    751         Y_i (t)  & = & y_i^0 + y_i^1 t^1 + \cdots + y_i^q t^q + o ( t^q )
    752         \\
    753         ty [ i * ( q + 1 ) + k ] & = & y_i^k
    754 \end{array}
    755 \] $$
    756 where $latex o( t^q ) / t^q \rightarrow 0$$ as $latex t \rightarrow 0$$.
    757 Note that superscripts represent an index for $latex y_j^k$$
    758 and an exponent for $latex t^k$$.
    759 Also note that the Taylor coefficients for $latex Y(t)$$ correspond
    760 to the derivatives of $latex Y(t)$$ at $latex t = 0$$ in the following way:
    761 $latex \[
    762         y_j^k = \frac{1}{ k ! } Y_j^{(k)} (0)
    763 \] $$
    764 If $latex p > 0$$,
    765 for $latex i = 0 , \ldots , m-1$$ and $latex k = 0 , \ldots , p-1$$,
    766 the input of $icode ty$$ satisfies
    767 $latex \[
    768         ty [ i * ( q + 1 ) + k ] = y_i^k
    769 \]$$
    770 and hence the corresponding elements need not be recalculated.
    771 
    772 $head ok$$
    773 If the required results are calculated, $icode ok$$ should be true.
    774 Otherwise, it should be false.
    775 
    776 $head Discussion$$
    777 For example, suppose that $icode%q% == 2%$$,
    778 and you know how to compute the function $latex f(x)$$,
    779 its first derivative $latex f^{(1)} (x)$$,
    780 and it component wise Hessian $latex f_i^{(2)} (x)$$.
    781 Then you can compute $icode ty$$ using the following formulas:
    782 $latex \[
    783 \begin{array}{rcl}
    784 y_i^0 & = & Y(0)
    785         = f_i ( x^0 )
    786 \\
    787 y_i^1 & = & Y^{(1)} ( 0 )
    788         = f_i^{(1)} ( x^0 ) X^{(1)} ( 0 )
    789         = f_i^{(1)} ( x^0 ) x^1
    790 \\
    791 y_i^2
    792 & = & \frac{1}{2 !} Y^{(2)} (0)
    793 \\
    794 & = & \frac{1}{2} X^{(1)} (0)^\R{T} f_i^{(2)} ( x^0 ) X^{(1)} ( 0 )
    795   +   \frac{1}{2} f_i^{(1)} ( x^0 ) X^{(2)} ( 0 )
    796 \\
    797 & = & \frac{1}{2} (x^1)^\R{T} f_i^{(2)} ( x^0 ) x^1
    798   +    f_i^{(1)} ( x^0 ) x^2
    799 \end{array}
    800 \] $$
    801 For $latex i = 0 , \ldots , m-1$$, and $latex k = 0 , 1 , 2$$,
    802 $latex \[
    803         ty [ i * (q + 1) + k ] = y_i^k
    804 \] $$
    805 
    806 $children%
    807         example/atomic/forward.cpp
    808 %$$
    809 $head Examples$$
    810 The file $cref atomic_forward.cpp$$ contains an example and test
    811 that uses this routine.
    812 It returns true if the test passes and false if it fails.
    813 
    814 $end
    815 -----------------------------------------------------------------------------
    816 */
    817 /*!
    818 Link from atomic_base to forward mode
    819 
    820 \param p [in]
    821 lowerest order for this forward mode calculation.
    822 
    823 \param q [in]
    824 highest order for this forward mode calculation.
    825 
    826 \param vx [in]
    827 if size not zero, which components of \c x are variables
    828 
    829 \param vy [out]
    830 if size not zero, which components of \c y are variables
    831 
    832 \param tx [in]
    833 Taylor coefficients corresponding to \c x for this calculation.
    834 
    835 \param ty [out]
    836 Taylor coefficient corresponding to \c y for this calculation
    837 
    838 See the forward mode in user's documentation for base_atomic
    839 */
    840 virtual bool forward(
    841         size_t                    p  ,
    842         size_t                    q  ,
    843         const vector<bool>&       vx ,
    844               vector<bool>&       vy ,
    845         const vector<Base>&       tx ,
    846               vector<Base>&       ty )
    847 {       return false; }
    848 /*
    849 -----------------------------------------------------------------------------
    850 $begin atomic_reverse$$
    851 $spell
    852         sq
    853         mul.hpp
    854         afun
    855         ty
    856         px
    857         py
    858         Taylor
    859         const
    860         CppAD
    861 $$
    862 
    863 $section Atomic Reverse Mode$$
    864 $spell
    865         bool
    866 $$
    867 
    868 $head Syntax$$
    869 $icode%ok% = %afun%.reverse(%q%, %tx%, %ty%, %px%, %py%)%$$
    870 
    871 $head Purpose$$
    872 This function is used by $cref/reverse/Reverse/$$
    873 to compute derivatives.
    874 
    875 $head Implementation$$
    876 If you are using
    877 $cref/reverse/Reverse/$$ mode,
    878 this virtual function must be defined by the
    879 $cref/atomic_user/atomic_ctor/atomic_user/$$ class.
    880 It can just return $icode%ok% == false%$$
    881 (and not compute anything) for values
    882 of $icode q$$ that are greater than those used by your
    883 $cref/reverse/Reverse/$$ mode calculations.
    884 
    885 $head q$$
    886 The argument $icode q$$ has prototype
    887 $codei%
    888         size_t %q%
    889 %$$
    890 It specifies the highest order Taylor coefficient that
    891 computing the derivative of.
    892 
    893 $head tx$$
    894 The argument $icode tx$$ has prototype
    895 $codei%
    896         const CppAD::vector<%Base%>& %tx%
    897 %$$
    898 and $icode%tx%.size() == (%q%+1)*%n%$$.
    899 For $latex j = 0 , \ldots , n-1$$ and $latex k = 0 , \ldots , q$$,
    900 we use the Taylor coefficient notation
    901 $latex \[
    902 \begin{array}{rcl}
    903         x_j^k    & = & tx [ j * ( q + 1 ) + k ]
    904         \\
    905         X_j (t)  & = & x_j^0 + x_j^1 t^1 + \cdots + x_j^q t^q
    906 \end{array}
    907 \] $$
    908 Note that superscripts represent an index for $latex x_j^k$$
    909 and an exponent for $latex t^k$$.
    910 Also note that the Taylor coefficients for $latex X(t)$$ correspond
    911 to the derivatives of $latex X(t)$$ at $latex t = 0$$ in the following way:
    912 $latex \[
    913         x_j^k = \frac{1}{ k ! } X_j^{(k)} (0)
    914 \] $$
    915 
    916 $head ty$$
    917 The argument $icode ty$$ has prototype
    918 $codei%
    919         const CppAD::vector<%Base%>& %ty%
    920 %$$
    921 and $icode%tx%.size() == (%q%+1)*%m%$$.
    922 For $latex i = 0 , \ldots , m-1$$ and $latex k = 0 , \ldots , q$$,
    923 we use the Taylor coefficient notation
    924 $latex \[
    925 \begin{array}{rcl}
    926         Y_i (t)  & = & f_i [ X(t) ]
    927         \\
    928         Y_i (t)  & = & y_i^0 + y_i^1 t^1 + \cdots + y_i^q t^q + o ( t^q )
    929         \\
    930         y_i^k    & = & ty [ i * ( q + 1 ) + k ]
    931 \end{array}
    932 \] $$
    933 where $latex o( t^q ) / t^q \rightarrow 0$$ as $latex t \rightarrow 0$$.
    934 Note that superscripts represent an index for $latex y_j^k$$
    935 and an exponent for $latex t^k$$.
    936 Also note that the Taylor coefficients for $latex Y(t)$$ correspond
    937 to the derivatives of $latex Y(t)$$ at $latex t = 0$$ in the following way:
    938 $latex \[
    939         y_j^k = \frac{1}{ k ! } Y_j^{(k)} (0)
    940 \] $$
    941 
    942 
    943 $head F$$
    944 We use the notation $latex \{ x_j^k \} \in B^{n \times (q+1)}$$ for
    945 $latex \[
    946         \{ x_j^k \W{:} j = 0 , \ldots , n-1, k = 0 , \ldots , q \}
    947 \]$$
    948 We use the notation $latex \{ y_i^k \} \in B^{m \times (q+1)}$$ for
    949 $latex \[
    950         \{ y_i^k \W{:} i = 0 , \ldots , m-1, k = 0 , \ldots , q \}
    951 \]$$
    952 We define the function
    953 $latex F : B^{n \times (q+1)} \rightarrow B^{m \times (q+1)}$$ by
    954 $latex \[
    955         y_i^k = F_i^k [ \{ x_j^k \} ]
    956 \] $$
    957 Note that
    958 $latex \[
    959         F_i^0 ( \{ x_j^k \} ) = f_i ( X(0) )  = f_i ( x^0 )
    960 \] $$
    961 We also note that
    962 $latex F_i^\ell ( \{ x_j^k \} )$$ is a function of
    963 $latex x^0 , \ldots , x^\ell$$
    964 and is determined by the derivatives of $latex f_i (x)$$
    965 up to order $latex \ell$$.
    966 
    967 
    968 $head G, H$$
    969 We use $latex G : B^{m \times (q+1)} \rightarrow B$$
    970 to denote an arbitrary scalar valued function of $latex \{ y_i^k \}$$.
    971 We use $latex H : B^{n \times (q+1)} \rightarrow B$$
    972 defined by
    973 $latex \[
    974         H ( \{ x_j^k \} ) = G[ F( \{ x_j^k \} ) ]
    975 \] $$
    976 
    977 $head py$$
    978 The argument $icode py$$ has prototype
    979 $codei%
    980         const CppAD::vector<%Base%>& %py%
    981 %$$
    982 and $icode%py%.size() == m * (%q%+1)%$$.
    983 For $latex i = 0 , \ldots , m-1$$, $latex k = 0 , \ldots , q$$,
    984 $latex \[
    985         py[ i * (q + 1 ) + k ] = \partial G / \partial y_i^k
    986 \] $$
    987 
    988 $subhead px$$
    989 The $icode px$$ has prototype
    990 $codei%
    991         CppAD::vector<%Base%>& %px%
    992 %$$
    993 and $icode%px%.size() == n * (%q%+1)%$$.
    994 The input values of the elements of $icode px$$
    995 are not specified (must not matter).
    996 Upon return,
    997 for $latex j = 0 , \ldots , n-1$$ and $latex \ell = 0 , \ldots , q$$,
    998 $latex \[
    999 \begin{array}{rcl}
    1000 px [ j * (q + 1) + \ell ] & = & \partial H / \partial x_j^\ell
    1001 \\
    1002 & = &
    1003 ( \partial G / \partial \{ y_i^k \} ) \cdot
    1004         ( \partial \{ y_i^k \} / \partial x_j^\ell )
    1005 \\
    1006 & = &
    1007 \sum_{k=0}^q
    1008 \sum_{i=0}^{m-1}
    1009 ( \partial G / \partial y_i^k ) ( \partial y_i^k / \partial x_j^\ell )
    1010 \\
    1011 & = &
    1012 \sum_{k=\ell}^q
    1013 \sum_{i=0}^{m-1}
    1014 py[ i * (q + 1 ) + k ] ( \partial F_i^k / \partial x_j^\ell )
    1015 \end{array}
    1016 \] $$
    1017 Note that we have used the fact that for $latex k < \ell$$,
    1018 $latex \partial F_i^k / \partial x_j^\ell = 0$$.
    1019 
    1020 $head ok$$
    1021 The return value $icode ok$$ has prototype
    1022 $codei%
    1023         bool %ok%
    1024 %$$
    1025 If it is $code true$$, the corresponding evaluation succeeded,
    1026 otherwise it failed.
    1027 
    1028 $children%
    1029         example/atomic/reverse.cpp
    1030 %$$
    1031 $head Examples$$
    1032 The file $cref atomic_reverse.cpp$$ contains an example and test
    1033 that uses this routine.
    1034 It returns true if the test passes and false if it fails.
    1035 
    1036 $end
    1037 -----------------------------------------------------------------------------
    1038 */
    1039 /*!
    1040 Link from reverse mode sweep to users routine.
    1041 
    1042 \param q [in]
    1043 highest order for this reverse mode calculation.
    1044 
    1045 \param tx [in]
    1046 Taylor coefficients corresponding to \c x for this calculation.
    1047 
    1048 \param ty [in]
    1049 Taylor coefficient corresponding to \c y for this calculation
    1050 
    1051 \param px [out]
    1052 Partials w.r.t. the \c x Taylor coefficients.
    1053 
    1054 \param py [in]
    1055 Partials w.r.t. the \c y Taylor coefficients.
    1056 
    1057 See atomic_reverse mode use documentation
    1058 */
    1059 virtual bool reverse(
    1060         size_t                    q  ,
    1061         const vector<Base>&       tx ,
    1062         const vector<Base>&       ty ,
    1063               vector<Base>&       px ,
    1064         const vector<Base>&       py )
    1065 {       return false; }
    1066 /*
    1067 -------------------------------------- ---------------------------------------
    1068 $begin atomic_for_sparse_jac$$
    1069 $spell
    1070         sq
    1071         mul.hpp
    1072         afun
    1073         Jacobian
    1074         jac
    1075         const
    1076         CppAD
    1077         std
    1078         bool
    1079         std
    1080 $$
    1081 
    1082 $section Atomic Forward Jacobian Sparsity Patterns$$
    1083 
    1084 $head Syntax$$
    1085 $icode%ok% = %afun%.for_sparse_jac(%q%, %r%, %s%, %x%)
    1086 %$$
    1087 
    1088 $head Deprecated 2016-06-27$$
    1089 $icode%ok% = %afun%.for_sparse_jac(%q%, %r%, %s%)
    1090 %$$
    1091 
    1092 $head Purpose$$
    1093 This function is used by $cref ForSparseJac$$ to compute
    1094 Jacobian sparsity patterns.
    1095 For a fixed matrix $latex R \in B^{n \times q}$$,
    1096 the Jacobian of $latex f( x + R * u)$$ with respect to $latex u \in B^q$$ is
    1097 $latex \[
    1098         S(x) = f^{(1)} (x) * R
    1099 \] $$
    1100 Given a $cref/sparsity pattern/glossary/Sparsity Pattern/$$ for $latex R$$,
    1101 $code for_sparse_jac$$ computes a sparsity pattern for $latex S(x)$$.
    1102 
    1103 $head Implementation$$
    1104 If you are using
    1105 $cref ForSparseJac$$,
    1106 $cref ForSparseHes$$, or
    1107 $cref RevSparseHes$$,
    1108 one of the versions of this
    1109 virtual function must be defined by the
    1110 $cref/atomic_user/atomic_ctor/atomic_user/$$ class.
    1111 
    1112 $subhead q$$
    1113 The argument $icode q$$ has prototype
    1114 $codei%
    1115      size_t %q%
    1116 %$$
    1117 It specifies the number of columns in
    1118 $latex R \in B^{n \times q}$$ and the Jacobian
    1119 $latex S(x) \in B^{m \times q}$$.
    1120 
    1121 $subhead r$$
    1122 This argument has prototype
    1123 $codei%
    1124      const %atomic_sparsity%& %r%
    1125 %$$
    1126 and is a $cref/atomic_sparsity/atomic_option/atomic_sparsity/$$ pattern for
    1127 $latex R \in B^{n \times q}$$.
    1128 
    1129 $subhead s$$
    1130 This argument has prototype
    1131 $codei%
    1132         %atomic_sparsity%& %s%
    1133 %$$
    1134 The input values of its elements
    1135 are not specified (must not matter).
    1136 Upon return, $icode s$$ is a
    1137 $cref/atomic_sparsity/atomic_option/atomic_sparsity/$$ pattern for
    1138 $latex S(x) \in B^{m \times q}$$.
    1139 
    1140 $subhead x$$
    1141 $index deprecated$$
    1142 The argument has prototype
    1143 $codei%
    1144         const CppAD::vector<%Base%>& %x%
    1145 %$$
    1146 and size is equal to the $icode n$$.
    1147 This is the $cref Value$$ value corresponding to the parameters in the
    1148 vector $cref/ax/atomic_afun/ax/$$ (when the atomic function was called).
    1149 To be specific, if
    1150 $codei%
    1151         if( Parameter(%ax%[%i%]) == true )
    1152                 %x%[%i%] = Value( %ax%[%i%] );
    1153         else
    1154                 %x%[%i%] = CppAD::numeric_limits<%Base%>::quiet_NaN();
    1155 %$$
    1156 The version of this function with out the $icode x$$ argument is deprecated;
    1157 i.e., you should include the argument even if you do not use it.
    1158 
    1159 $head ok$$
    1160 The return value $icode ok$$ has prototype
    1161 $codei%
    1162         bool %ok%
    1163 %$$
    1164 If it is $code true$$, the corresponding evaluation succeeded,
    1165 otherwise it failed.
    1166 
    1167 $children%
    1168         example/atomic/for_sparse_jac.cpp
    1169 %$$
    1170 $head Examples$$
    1171 The file $cref atomic_for_sparse_jac.cpp$$ contains an example and test
    1172 that uses this routine.
    1173 It returns true if the test passes and false if it fails.
    1174 
    1175 $end
    1176 -----------------------------------------------------------------------------
    1177 */
    1178 /*!
    1179 Link, after case split, from for_jac_sweep to atomic_base.
    1180 
    1181 \param q
    1182 is the column dimension for the Jacobian sparsity partterns.
    1183 
    1184 \param r
    1185 is the Jacobian sparsity pattern for the argument vector x
    1186 
    1187 \param s
    1188 is the Jacobian sparsity pattern for the result vector y
    1189 
    1190 \param x
    1191 is the integer value for x arguments that are parameters.
    1192 */
    1193 virtual bool for_sparse_jac(
    1194         size_t                                  q  ,
    1195         const vector< std::set<size_t> >&       r  ,
    1196               vector< std::set<size_t> >&       s  ,
    1197         const vector<Base>&                     x  )
    1198 {       return false; }
    1199 virtual bool for_sparse_jac(
    1200         size_t                                  q  ,
    1201         const vector<bool>&                     r  ,
    1202               vector<bool>&                     s  ,
    1203         const vector<Base>&                     x  )
    1204 {       return false; }
    1205 virtual bool for_sparse_jac(
    1206         size_t                                  q  ,
    1207         const vectorBool&                       r  ,
    1208               vectorBool&                       s  ,
    1209         const vector<Base>&                     x  )
    1210 {       return false; }
    1211 // deprecated versions
    1212 virtual bool for_sparse_jac(
    1213         size_t                                  q  ,
    1214         const vector< std::set<size_t> >&       r  ,
    1215               vector< std::set<size_t> >&       s  )
    1216 {       return false; }
    1217 virtual bool for_sparse_jac(
    1218         size_t                                  q  ,
    1219         const vector<bool>&                     r  ,
    1220               vector<bool>&                     s  )
    1221 {       return false; }
    1222 virtual bool for_sparse_jac(
    1223         size_t                                  q  ,
    1224         const vectorBool&                       r  ,
    1225               vectorBool&                       s  )
    1226 {       return false; }
    1227 
    1228 /*!
    1229 Link, before case split, from for_jac_sweep to atomic_base.
    1230 
    1231 \tparam InternalSparsity
    1232 Is the used internaly for sparsity calculations; i.e.,
    1233 sparse_pack or sparse_list.
    1234 
    1235 \param x
    1236 is parameter arguments to the function, other components are nan.
    1237 
    1238 \param x_index
    1239 is the variable index, on the tape, for the arguments to this function.
    1240 This size of x_index is n, the number of arguments to this function.
    1241 
    1242 \param y_index
    1243 is the variable index, on the tape, for the results for this function.
    1244 This size of y_index is m, the number of results for this function.
    1245 
    1246 \param var_sparsity
    1247 On input, for j = 0, ... , n-1, the sparsity pattern with index x_index[j],
    1248 is the sparsity for the j-th argument to this atomic function.
    1249 On input, for i = 0, ... , m-1, the sparsity pattern with index y_index[i],
    1250 is empty. On output, it is the sparsity
    1251 for the j-th result for this atomic function.
    1252 */
    1253 template <class InternalSparsity>
    1254 void for_sparse_jac(
    1255         const vector<Base>&        x            ,
    1256         const vector<size_t>&      x_index      ,
    1257         const vector<size_t>&      y_index      ,
    1258         InternalSparsity&          var_sparsity )
    1259 {
    1260         // intial results are empty during forward mode
    1261         size_t q           = var_sparsity.end();
    1262         bool   input_empty = true;
    1263         bool   zero_empty  = true;
    1264         bool   transpose   = false;
    1265         size_t m           = y_index.size();
    1266         bool   ok          = false;
    1267         size_t thread      = thread_alloc::thread_num();
    1268         allocate_work(thread);
    1269         //
    1270         std::string msg    = ": atomic_base.for_sparse_jac: returned false";
    1271         if( sparsity_ == pack_sparsity_enum )
    1272         {       vectorBool& pack_r ( work_[thread]->pack_r );
    1273                 vectorBool& pack_s ( work_[thread]->pack_s );
    1274                 local::get_internal_sparsity(
    1275                         transpose, x_index, var_sparsity, pack_r
    1276                 );
    1277                 //
    1278                 pack_s.resize(m * q );
    1279                 ok = for_sparse_jac(q, pack_r, pack_s, x);
    1280                 if( ! ok )
    1281                         ok = for_sparse_jac(q, pack_r, pack_s);
    1282                 if( ! ok )
    1283                 {       msg = afun_name() + msg + " sparsity = pack_sparsity_enum";
    1284                         CPPAD_ASSERT_KNOWN(false, msg.c_str());
    1285                 }
    1286                 local::set_internal_sparsity(zero_empty, input_empty,
    1287                         transpose, y_index, var_sparsity, pack_s
    1288                 );
    1289         }
    1290         else if( sparsity_ == bool_sparsity_enum )
    1291         {       vector<bool>& bool_r ( work_[thread]->bool_r );
    1292                 vector<bool>& bool_s ( work_[thread]->bool_s );
    1293                 local::get_internal_sparsity(
    1294                         transpose, x_index, var_sparsity, bool_r
    1295                 );
    1296                 bool_s.resize(m * q );
    1297                 ok = for_sparse_jac(q, bool_r, bool_s, x);
    1298                 if( ! ok )
    1299                         ok = for_sparse_jac(q, bool_r, bool_s);
    1300                 if( ! ok )
    1301                 {       msg = afun_name() + msg + " sparsity = bool_sparsity_enum";
    1302                         CPPAD_ASSERT_KNOWN(false, msg.c_str());
    1303                 }
    1304                 local::set_internal_sparsity(zero_empty, input_empty,
    1305                         transpose, y_index, var_sparsity, bool_s
    1306                 );
    1307         }
    1308         else
    1309         {       CPPAD_ASSERT_UNKNOWN( sparsity_ == set_sparsity_enum );
    1310                 vector< std::set<size_t> >& set_r ( work_[thread]->set_r );
    1311                 vector< std::set<size_t> >& set_s ( work_[thread]->set_s );
    1312                 local::get_internal_sparsity(
    1313                         transpose, x_index, var_sparsity, set_r
    1314                 );
    1315                 //
    1316                 set_s.resize(m);
    1317                 ok = for_sparse_jac(q, set_r, set_s, x);
    1318                 if( ! ok )
    1319                         ok = for_sparse_jac(q, set_r, set_s);
    1320                 if( ! ok )
    1321                 {       msg = afun_name() + msg + " sparsity = set_sparsity_enum";
    1322                         CPPAD_ASSERT_KNOWN(false, msg.c_str());
    1323                 }
    1324                 local::set_internal_sparsity(zero_empty, input_empty,
    1325                         transpose, y_index, var_sparsity, set_s
    1326                 );
    1327         }
    1328         return;
    1329 }
    1330 /*
    1331 -------------------------------------- ---------------------------------------
    1332 $begin atomic_rev_sparse_jac$$
    1333 $spell
    1334         sq
    1335         mul.hpp
    1336         rt
    1337         afun
    1338         Jacobian
    1339         jac
    1340         CppAD
    1341         std
    1342         bool
    1343         const
    1344         hes
    1345 $$
    1346 
    1347 $section Atomic Reverse Jacobian Sparsity Patterns$$
    1348 
    1349 $head Syntax$$
    1350 $icode%ok% = %afun%.rev_sparse_jac(%q%, %rt%, %st%, %x%)
    1351 %$$
    1352 
    1353 $head Deprecated 2016-06-27$$
    1354 $icode%ok% = %afun%.rev_sparse_jac(%q%, %rt%, %st%)
    1355 %$$
    1356 
    1357 $head Purpose$$
    1358 This function is used by
    1359 $cref RevSparseJac$$ to compute
    1360 Jacobian sparsity patterns.
    1361 If you are using $cref RevSparseJac$$,
    1362 one of the versions of this
    1363 virtual function must be defined by the
    1364 $cref/atomic_user/atomic_ctor/atomic_user/$$ class.
    1365 $pre
    1366 
    1367 $$
    1368 For a fixed matrix $latex R \in B^{q \times m}$$,
    1369 the Jacobian of $latex R * f( x )$$ with respect to $latex x \in B^n$$ is
    1370 $latex \[
    1371         S(x) = R * f^{(1)} (x)
    1372 \] $$
    1373 Given a $cref/sparsity pattern/glossary/Sparsity Pattern/$$ for $latex R$$,
    1374 $code rev_sparse_jac$$ computes a sparsity pattern for $latex S(x)$$.
    1375 
    1376 $head Implementation$$
    1377 If you are using
    1378 $cref RevSparseJac$$ or $cref ForSparseHes$$,
    1379 this virtual function must be defined by the
    1380 $cref/atomic_user/atomic_ctor/atomic_user/$$ class.
    1381 
    1382 $subhead q$$
    1383 The argument $icode q$$ has prototype
    1384 $codei%
    1385      size_t %q%
    1386 %$$
    1387 It specifies the number of rows in
    1388 $latex R \in B^{q \times m}$$ and the Jacobian
    1389 $latex S(x) \in B^{q \times n}$$.
    1390 
    1391 $subhead rt$$
    1392 This argument has prototype
    1393 $codei%
    1394      const %atomic_sparsity%& %rt%
    1395 %$$
    1396 and is a
    1397 $cref/atomic_sparsity/atomic_option/atomic_sparsity/$$ pattern for
    1398 $latex R^\R{T} \in B^{m \times q}$$.
    1399 
    1400 $subhead st$$
    1401 This argument has prototype
    1402 $codei%
    1403         %atomic_sparsity%& %st%
    1404 %$$
    1405 The input value of its elements
    1406 are not specified (must not matter).
    1407 Upon return, $icode s$$ is a
    1408 $cref/atomic_sparsity/atomic_option/atomic_sparsity/$$ pattern for
    1409 $latex S(x)^\R{T} \in B^{n \times q}$$.
    1410 
    1411 $subhead x$$
    1412 $index deprecated$$
    1413 The argument has prototype
    1414 $codei%
    1415         const CppAD::vector<%Base%>& %x%
    1416 %$$
    1417 and size is equal to the $icode n$$.
    1418 This is the $cref Value$$ corresponding to the parameters in the
    1419 vector $cref/ax/atomic_afun/ax/$$ (when the atomic function was called).
    1420 To be specific, if
    1421 $codei%
    1422         if( Parameter(%ax%[%i%]) == true )
    1423                 %x%[%i%] = Value( %ax%[%i%] );
    1424         else
    1425                 %x%[%i%] = CppAD::numeric_limits<%Base%>::quiet_NaN();
    1426 %$$
    1427 The version of this function with out the $icode x$$ argument is deprecated;
    1428 i.e., you should include the argument even if you do not use it.
    1429 
    1430 $head ok$$
    1431 The return value $icode ok$$ has prototype
    1432 $codei%
    1433         bool %ok%
    1434 %$$
    1435 If it is $code true$$, the corresponding evaluation succeeded,
    1436 otherwise it failed.
    1437 
    1438 $children%
    1439         example/atomic/rev_sparse_jac.cpp
    1440 %$$
    1441 $head Examples$$
    1442 The file $cref atomic_rev_sparse_jac.cpp$$ contains an example and test
    1443 that uses this routine.
    1444 It returns true if the test passes and false if it fails.
    1445 
    1446 $end
    1447 -----------------------------------------------------------------------------
    1448 */
    1449 /*!
    1450 Link, after case split, from rev_jac_sweep to atomic_base
    1451 
    1452 \param q [in]
    1453 is the row dimension for the Jacobian sparsity partterns
    1454 
    1455 \param rt [out]
    1456 is the tansposed Jacobian sparsity pattern w.r.t to range variables y
    1457 
    1458 \param st [in]
    1459 is the tansposed Jacobian sparsity pattern for the argument variables x
    1460 
    1461 \param x
    1462 is the integer value for x arguments that are parameters.
    1463 */
    1464 virtual bool rev_sparse_jac(
    1465         size_t                                  q  ,
    1466         const vector< std::set<size_t> >&       rt ,
    1467               vector< std::set<size_t> >&       st ,
    1468         const vector<Base>&                     x  )
    1469 {       return false; }
    1470 virtual bool rev_sparse_jac(
    1471         size_t                                  q  ,
    1472         const vector<bool>&                     rt ,
    1473               vector<bool>&                     st ,
    1474         const vector<Base>&                     x  )
    1475 {       return false; }
    1476 virtual bool rev_sparse_jac(
    1477         size_t                                  q  ,
    1478         const vectorBool&                       rt ,
    1479               vectorBool&                       st ,
    1480         const vector<Base>&                     x  )
    1481 {       return false; }
    1482 // deprecated versions
    1483 virtual bool rev_sparse_jac(
    1484         size_t                                  q  ,
    1485         const vector< std::set<size_t> >&       rt ,
    1486               vector< std::set<size_t> >&       st )
    1487 {       return false; }
    1488 virtual bool rev_sparse_jac(
    1489         size_t                                  q  ,
    1490         const vector<bool>&                     rt ,
    1491               vector<bool>&                     st )
    1492 {       return false; }
    1493 virtual bool rev_sparse_jac(
    1494         size_t                                  q  ,
    1495         const vectorBool&                       rt ,
    1496               vectorBool&                       st )
    1497 {       return false; }
    1498 
    1499 /*!
    1500 Link, before case split, from rev_jac_sweep to atomic_base.
    1501 
    1502 \tparam InternalSparsity
    1503 Is the used internaly for sparsity calculations; i.e.,
    1504 sparse_pack or sparse_list.
    1505 
    1506 \param x
    1507 is parameter arguments to the function, other components are nan.
    1508 
    1509 \param x_index
    1510 is the variable index, on the tape, for the arguments to this function.
    1511 This size of x_index is n, the number of arguments to this function.
    1512 
    1513 \param y_index
    1514 is the variable index, on the tape, for the results for this function.
    1515 This size of y_index is m, the number of results for this function.
    1516 
    1517 \param var_sparsity
    1518 On input, for i = 0, ... , m-1, the sparsity pattern with index y_index[i],
    1519 is the sparsity for the i-th argument to this atomic function.
    1520 On output, for j = 0, ... , n-1, the sparsity pattern with index x_index[j],
    1521 the sparsity has been updated to remove y as a function of x.
    1522 */
    1523 template <class InternalSparsity>
    1524 void rev_sparse_jac(
    1525         const vector<Base>&        x            ,
    1526         const vector<size_t>&      x_index      ,
    1527         const vector<size_t>&      y_index      ,
    1528         InternalSparsity&          var_sparsity )
    1529 {
    1530         // initial results may be non-empty during reverse mode
    1531         size_t q           = var_sparsity.end();
    1532         bool   input_empty = false;
    1533         bool   zero_empty  = true;
    1534         bool   transpose   = false;
    1535         size_t n           = x_index.size();
    1536         bool   ok          = false;
    1537         size_t thread      = thread_alloc::thread_num();
    1538         allocate_work(thread);
    1539         //
    1540         std::string msg    = ": atomic_base.rev_sparse_jac: returned false";
    1541         if( sparsity_ == pack_sparsity_enum )
    1542         {       vectorBool& pack_rt ( work_[thread]->pack_r );
    1543                 vectorBool& pack_st ( work_[thread]->pack_s );
    1544                 local::get_internal_sparsity(
    1545                         transpose, y_index, var_sparsity, pack_rt
    1546                 );
    1547                 //
    1548                 pack_st.resize(n * q );
    1549                 ok = rev_sparse_jac(q, pack_rt, pack_st, x);
    1550                 if( ! ok )
    1551                         ok = rev_sparse_jac(q, pack_rt, pack_st);
    1552                 if( ! ok )
    1553                 {       msg = afun_name() + msg + " sparsity = pack_sparsity_enum";
    1554                         CPPAD_ASSERT_KNOWN(false, msg.c_str());
    1555                 }
    1556                 local::set_internal_sparsity(zero_empty, input_empty,
    1557                         transpose, x_index, var_sparsity, pack_st
    1558                 );
    1559         }
    1560         else if( sparsity_ == bool_sparsity_enum )
    1561         {       vector<bool>& bool_rt ( work_[thread]->bool_r );
    1562                 vector<bool>& bool_st ( work_[thread]->bool_s );
    1563                 local::get_internal_sparsity(
    1564                         transpose, y_index, var_sparsity, bool_rt
    1565                 );
    1566                 bool_st.resize(n * q );
    1567                 ok = rev_sparse_jac(q, bool_rt, bool_st, x);
    1568                 if( ! ok )
    1569                         ok = rev_sparse_jac(q, bool_rt, bool_st);
    1570                 if( ! ok )
    1571                 {       msg = afun_name() + msg + " sparsity = bool_sparsity_enum";
    1572                         CPPAD_ASSERT_KNOWN(false, msg.c_str());
    1573                 }
    1574                 local::set_internal_sparsity(zero_empty, input_empty,
    1575                         transpose, x_index, var_sparsity, bool_st
    1576                 );
    1577         }
    1578         else
    1579         {       CPPAD_ASSERT_UNKNOWN( sparsity_ == set_sparsity_enum );
    1580                 vector< std::set<size_t> >& set_rt ( work_[thread]->set_r );
    1581                 vector< std::set<size_t> >& set_st ( work_[thread]->set_s );
    1582                 local::get_internal_sparsity(
    1583                         transpose, y_index, var_sparsity, set_rt
    1584                 );
    1585                 set_st.resize(n);
    1586                 ok = rev_sparse_jac(q, set_rt, set_st, x);
    1587                 if( ! ok )
    1588                         ok = rev_sparse_jac(q, set_rt, set_st);
    1589                 if( ! ok )
    1590                 {       msg = afun_name() + msg + " sparsity = set_sparsity_enum";
    1591                         CPPAD_ASSERT_KNOWN(false, msg.c_str());
    1592                 }
    1593                 local::set_internal_sparsity(zero_empty, input_empty,
    1594                         transpose, x_index, var_sparsity, set_st
    1595                 );
    1596         }
    1597         return;
    1598 }
    1599 /*
    1600 -------------------------------------- ---------------------------------------
    1601 $begin atomic_for_sparse_hes$$
    1602 $spell
    1603         sq
    1604         mul.hpp
    1605         vx
    1606         afun
    1607         Jacobian
    1608         jac
    1609         CppAD
    1610         std
    1611         bool
    1612         hes
    1613         const
    1614 $$
    1615 
    1616 $section Atomic Forward Hessian Sparsity Patterns$$
    1617 
    1618 $head Syntax$$
    1619 $icode%ok% = %afun%.for_sparse_hes(%vx%, %r%, %s%, %h%, %x%)%$$
    1620 
    1621 $head Deprecated 2016-06-27$$
    1622 $icode%ok% = %afun%.for_sparse_hes(%vx%, %r%, %s%, %h%)%$$
    1623 
    1624 $head Purpose$$
    1625 This function is used by $cref ForSparseHes$$ to compute
    1626 Hessian sparsity patterns.
    1627 If you are using $cref ForSparseHes$$,
    1628 one of the versions of this
    1629 virtual function must be defined by the
    1630 $cref/atomic_user/atomic_ctor/atomic_user/$$ class.
    1631 $pre
    1632 
    1633 $$
    1634 Given a $cref/sparsity pattern/glossary/Sparsity Pattern/$$ for
    1635 a diagonal matrix $latex R \in B^{n \times n}$$, and
    1636 a row vector $latex S \in B^{1 \times m}$$,
    1637 this routine computes the sparsity pattern for
    1638 $latex \[
    1639         H(x) = R^\R{T} \cdot (S \cdot f)^{(2)}( x ) \cdot R
    1640 \] $$
    1641 
    1642 $head Implementation$$
    1643 If you are using and $cref ForSparseHes$$,
    1644 this virtual function must be defined by the
    1645 $cref/atomic_user/atomic_ctor/atomic_user/$$ class.
    1646 
    1647 $subhead vx$$
    1648 The argument $icode vx$$ has prototype
    1649 $codei%
    1650      const CppAD:vector<bool>& %vx%
    1651 %$$
    1652 $icode%vx%.size() == %n%$$, and
    1653 for $latex j = 0 , \ldots , n-1$$,
    1654 $icode%vx%[%j%]%$$ is true if and only if
    1655 $icode%ax%[%j%]%$$ is a $cref/variable/glossary/Variable/$$
    1656 or $cref/dynamic parameter/glossary/Parameter/Dynamic/$$
    1657 in the corresponding call to
    1658 $codei%
    1659         %afun%(%ax%, %ay%)
    1660 %$$
    1661 
    1662 $subhead r$$
    1663 This argument has prototype
    1664 $codei%
    1665      const CppAD:vector<bool>& %r%
    1666 %$$
    1667 and is a $cref/atomic_sparsity/atomic_option/atomic_sparsity/$$ pattern for
    1668 the diagonal of $latex R \in B^{n \times n}$$.
    1669 
    1670 $subhead s$$
    1671 The argument $icode s$$ has prototype
    1672 $codei%
    1673      const CppAD:vector<bool>& %s%
    1674 %$$
    1675 and its size is $icode m$$.
    1676 It is a sparsity pattern for $latex S \in B^{1 \times m}$$.
    1677 
    1678 $subhead h$$
    1679 This argument has prototype
    1680 $codei%
    1681      %atomic_sparsity%& %h%
    1682 %$$
    1683 The input value of its elements
    1684 are not specified (must not matter).
    1685 Upon return, $icode h$$ is a
    1686 $cref/atomic_sparsity/atomic_option/atomic_sparsity/$$ pattern for
    1687 $latex H(x) \in B^{n \times n}$$ which is defined above.
    1688 
    1689 $subhead x$$
    1690 $index deprecated$$
    1691 The argument has prototype
    1692 $codei%
    1693         const CppAD::vector<%Base%>& %x%
    1694 %$$
    1695 and size is equal to the $icode n$$.
    1696 This is the $cref Value$$ value corresponding to the parameters in the
    1697 vector $cref/ax/atomic_afun/ax/$$ (when the atomic function was called).
    1698 To be specific, if
    1699 $codei%
    1700         if( Parameter(%ax%[%i%]) == true )
    1701                 %x%[%i%] = Value( %ax%[%i%] );
    1702         else
    1703                 %x%[%i%] = CppAD::numeric_limits<%Base%>::quiet_NaN();
    1704 %$$
    1705 The version of this function with out the $icode x$$ argument is deprecated;
    1706 i.e., you should include the argument even if you do not use it.
    1707 
    1708 $children%
    1709         example/atomic/for_sparse_hes.cpp
    1710 %$$
    1711 $head Examples$$
    1712 The file $cref atomic_for_sparse_hes.cpp$$ contains an example and test
    1713 that uses this routine.
    1714 It returns true if the test passes and false if it fails.
    1715 
    1716 $end
    1717 -----------------------------------------------------------------------------
    1718 */
    1719 /*!
    1720 Link, after case split, from for_hes_sweep to atomic_base.
    1721 
    1722 \param vx [in]
    1723 which componens of x are variables.
    1724 
    1725 \param r [in]
    1726 is the forward Jacobian sparsity pattern w.r.t the argument vector x.
    1727 
    1728 \param s [in]
    1729 is the reverse Jacobian sparsity pattern w.r.t the result vector y.
    1730 
    1731 \param h [out]
    1732 is the Hessian sparsity pattern w.r.t the argument vector x.
    1733 
    1734 \param x
    1735 is the integer value of the x arguments that are parameters.
    1736 */
    1737 virtual bool for_sparse_hes(
    1738         const vector<bool>&             vx ,
    1739         const vector<bool>&             r  ,
    1740         const vector<bool>&             s  ,
    1741         vector< std::set<size_t> >&     h  ,
    1742         const vector<Base>&             x  )
    1743 {       return false; }
    1744 virtual bool for_sparse_hes(
    1745         const vector<bool>&             vx ,
    1746         const vector<bool>&             r  ,
    1747         const vector<bool>&             s  ,
    1748         vector<bool>&                   h  ,
    1749         const vector<Base>&             x  )
    1750 {       return false; }
    1751 virtual bool for_sparse_hes(
    1752         const vector<bool>&             vx ,
    1753         const vector<bool>&             r  ,
    1754         const vector<bool>&             s  ,
    1755         vectorBool&                     h  ,
    1756         const vector<Base>&             x  )
    1757 // deprecated
    1758 {       return false; }
    1759 virtual bool for_sparse_hes(
    1760         const vector<bool>&             vx ,
    1761         const vector<bool>&             r  ,
    1762         const vector<bool>&             s  ,
    1763         vector< std::set<size_t> >&     h  )
    1764 {       return false; }
    1765 // deprecated versions
    1766 virtual bool for_sparse_hes(
    1767         const vector<bool>&             vx ,
    1768         const vector<bool>&             r  ,
    1769         const vector<bool>&             s  ,
    1770         vector<bool>&                   h  )
    1771 {       return false; }
    1772 virtual bool for_sparse_hes(
    1773         const vector<bool>&             vx ,
    1774         const vector<bool>&             r  ,
    1775         const vector<bool>&             s  ,
    1776         vectorBool&                     h  )
    1777 {       return false; }
    1778 /*!
    1779 Link, before case split, from for_hes_sweep to atomic_base.
    1780 
    1781 \tparam InternalSparsity
    1782 Is the used internaly for sparsity calculations; i.e.,
    1783 sparse_pack or sparse_list.
    1784 
    1785 \param x
    1786 is parameter arguments to the function, other components are nan.
    1787 
    1788 \param x_index
    1789 is the variable index, on the tape, for the arguments to this function.
    1790 This size of x_index is n, the number of arguments to this function.
    1791 
    1792 \param y_index
    1793 is the variable index, on the tape, for the results for this function.
    1794 This size of y_index is m, the number of results for this function.
    1795 
    1796 \param for_jac_sparsity
    1797 On input, for j = 0, ... , n-1, the sparsity pattern with index x_index[j],
    1798 is the forward Jacobian sparsity for the j-th argument to this atomic function.
    1799 
    1800 \param rev_jac_sparsity
    1801 On input, for i = 0, ... , m-1, the sparsity pattern with index y_index[i],
    1802 is the reverse Jacobian sparsity for the i-th result to this atomic function.
    1803 This shows which components of the result affect the function we are
    1804 computing the Hessian of.
    1805 
    1806 \param for_hes_sparsity
    1807 This is the sparsity pattern for the Hessian. On input, the non-linear
    1808 terms in the atomic fuction have not been included. Upon return, they
    1809 have been included.
    1810 */
    1811 template <class InternalSparsity>
    1812 void for_sparse_hes(
    1813         const vector<Base>&        x                ,
    1814         const vector<size_t>&      x_index          ,
    1815         const vector<size_t>&      y_index          ,
    1816         const InternalSparsity&    for_jac_sparsity ,
    1817         const InternalSparsity&    rev_jac_sparsity ,
    1818         InternalSparsity&          for_hes_sparsity )
    1819 {       typedef typename InternalSparsity::const_iterator const_iterator;
    1820         CPPAD_ASSERT_UNKNOWN( rev_jac_sparsity.end() == 1 );
    1821         size_t n      = x_index.size();
    1822         size_t m      = y_index.size();
    1823         bool   ok     = false;
    1824         size_t thread = thread_alloc::thread_num();
    1825         allocate_work(thread);
    1826         //
    1827         // vx
    1828         vector<bool> vx(n);
    1829         for(size_t j = 0; j < n; j++)
    1830                 vx[j] = x_index[j] != 0;
    1831         //
    1832         // bool_r
    1833         vector<bool>& bool_r( work_[thread]->bool_r );
    1834         bool_r.resize(n);
    1835         for(size_t j = 0; j < n; j++)
    1836         {       // check if we must compute row and column j of h
    1837                 const_iterator itr(for_jac_sparsity, x_index[j]);
    1838                 size_t i = *itr;
    1839                 bool_r[j] = i < for_jac_sparsity.end();
    1840         }
    1841         //
    1842         // bool s
    1843         vector<bool>& bool_s( work_[thread]->bool_s );
    1844         bool_s.resize(m);
    1845         for(size_t i = 0; i < m; i++)
    1846         {       // check if row i of result is included in h
    1847                 bool_s[i] = rev_jac_sparsity.is_element(y_index[i], 0);
    1848         }
    1849         //
    1850         // h
    1851         vectorBool&                 pack_h( work_[thread]->pack_h );
    1852         vector<bool>&               bool_h( work_[thread]->bool_h );
    1853         vector< std::set<size_t> >& set_h(  work_[thread]->set_h );
    1854         //
    1855         // call user's version of atomic function
    1856         std::string msg    = ": atomic_base.for_sparse_hes: returned false";
    1857         if( sparsity_ == pack_sparsity_enum )
    1858         {       pack_h.resize(n * n);
    1859                 ok = for_sparse_hes(vx, bool_r, bool_s, pack_h, x);
    1860                 if( ! ok )
    1861                         ok = for_sparse_hes(vx, bool_r, bool_s, pack_h);
    1862                 if( ! ok )
    1863                 {       msg = afun_name() + msg + " sparsity = pack_sparsity_enum";
    1864                         CPPAD_ASSERT_KNOWN(false, msg.c_str());
    1865                 }
    1866         }
    1867         else if( sparsity_ == bool_sparsity_enum )
    1868         {       bool_h.resize(n * n);
    1869                 ok = for_sparse_hes(vx, bool_r, bool_s, bool_h, x);
    1870                 if( ! ok )
    1871                         ok = for_sparse_hes(vx, bool_r, bool_s, bool_h);
    1872                 if( ! ok )
    1873                 {       msg = afun_name() + msg + " sparsity = bool_sparsity_enum";
    1874                         CPPAD_ASSERT_KNOWN(false, msg.c_str());
    1875                 }
    1876         }
    1877         else
    1878         {       CPPAD_ASSERT_UNKNOWN( sparsity_ == set_sparsity_enum )
    1879                 set_h.resize(n);
    1880                 ok = for_sparse_hes(vx, bool_r, bool_s, set_h, x);
    1881                 if( ! ok )
    1882                         ok = for_sparse_hes(vx, bool_r, bool_s, set_h);
    1883                 if( ! ok )
    1884                 {       msg = afun_name() + msg + " sparsity = set_sparsity_enum";
    1885                         CPPAD_ASSERT_KNOWN(false, msg.c_str());
    1886                 }
    1887         }
    1888         CPPAD_ASSERT_UNKNOWN( ok );
    1889         //
    1890         // modify hessian in calling routine
    1891         for(size_t i = 0; i < n; i++)
    1892         {       for(size_t j = 0; j < n; j++)
    1893                 {       if( (x_index[i] > 0) & (x_index[j] > 0) )
    1894                         {       bool flag = false;
    1895                                 switch( sparsity_ )
    1896                                 {       case pack_sparsity_enum:
    1897                                         flag = pack_h[i * n + j];
    1898                                         break;
    1899                                         //
    1900                                         case bool_sparsity_enum:
    1901                                         flag = bool_h[i * n + j];
    1902                                         break;
    1903                                         //
    1904                                         case set_sparsity_enum:
    1905                                         flag = set_h[i].find(j) != set_h[i].end();
    1906                                         break;
    1907                                 }
    1908                                 if( flag )
    1909                                 {       const_iterator itr_i(for_jac_sparsity, x_index[i]);
    1910                                         size_t i_x = *itr_i;
    1911                                         while( i_x < for_jac_sparsity.end() )
    1912                                         {       for_hes_sparsity.binary_union(
    1913                                                         i_x, i_x, x_index[j], for_jac_sparsity
    1914                                                 );
    1915                                                 i_x = *(++itr_i);
    1916                                         }
    1917                                         const_iterator itr_j(for_jac_sparsity, x_index[j]);
    1918                                         size_t j_x = *itr_j;
    1919                                         while( j_x < for_jac_sparsity.end() )
    1920                                         {       for_hes_sparsity.binary_union(
    1921                                                         j_x, j_x, x_index[i], for_jac_sparsity
    1922                                                 );
    1923                                                 j_x = *(++itr_j);
    1924                                         }
    1925                                 }
    1926                         }
    1927                 }
    1928         }
    1929         return;
    1930 }
    1931 /*
    1932 -------------------------------------- ---------------------------------------
    1933 $begin atomic_rev_sparse_hes$$
    1934 $spell
    1935         sq
    1936         mul.hpp
    1937         vx
    1938         afun
    1939         Jacobian
    1940         jac
    1941         CppAD
    1942         std
    1943         bool
    1944         hes
    1945         const
    1946 $$
    1947 
    1948 $section Atomic Reverse Hessian Sparsity Patterns$$
    1949 
    1950 $head Syntax$$
    1951 $icode%ok% = %afun%.rev_sparse_hes(%vx%, %s%, %t%, %q%, %r%, %u%, %v%, %x%)%$$
    1952 
    1953 $head Deprecated 2016-06-27$$
    1954 $icode%ok% = %afun%.rev_sparse_hes(%vx%, %s%, %t%, %q%, %r%, %u%, %v%)%$$
    1955 
    1956 $head Purpose$$
    1957 This function is used by $cref RevSparseHes$$ to compute
    1958 Hessian sparsity patterns.
    1959 If you are using $cref RevSparseHes$$ to compute
    1960 one of the versions of this
    1961 virtual function muse be defined by the
    1962 $cref/atomic_user/atomic_ctor/atomic_user/$$ class.
    1963 $pre
    1964 
    1965 $$
    1966 There is an unspecified scalar valued function
    1967 $latex g : B^m \rightarrow B$$.
    1968 Given a $cref/sparsity pattern/glossary/Sparsity Pattern/$$ for
    1969 $latex R \in B^{n \times q}$$,
    1970 and information about the function $latex z = g(y)$$,
    1971 this routine computes the sparsity pattern for
    1972 $latex \[
    1973         V(x) = (g \circ f)^{(2)}( x ) R
    1974 \] $$
    1975 
    1976 $head Implementation$$
    1977 If you are using and $cref RevSparseHes$$,
    1978 this virtual function must be defined by the
    1979 $cref/atomic_user/atomic_ctor/atomic_user/$$ class.
    1980 
    1981 $subhead vx$$
    1982 The argument $icode vx$$ has prototype
    1983 $codei%
    1984      const CppAD:vector<bool>& %vx%
    1985 %$$
    1986 $icode%vx%.size() == %n%$$, and
    1987 for $latex j = 0 , \ldots , n-1$$,
    1988 $icode%vx%[%j%]%$$ is true if and only if
    1989 $icode%ax%[%j%]%$$ is a $cref/variable/glossary/Variable/$$
    1990 or $cref/dynamic parameter/glossary/Parameter/Dynamic/$$
    1991 in the corresponding call to
    1992 $codei%
    1993         %afun%(%ax%, %ay%)
    1994 %$$
    1995 
    1996 $subhead s$$
    1997 The argument $icode s$$ has prototype
    1998 $codei%
    1999      const CppAD:vector<bool>& %s%
    2000 %$$
    2001 and its size is $icode m$$.
    2002 It is a sparsity pattern for
    2003 $latex S(x) = g^{(1)} [ f(x) ] \in B^{1 \times m}$$.
    2004 
    2005 $subhead t$$
    2006 This argument has prototype
    2007 $codei%
    2008      CppAD:vector<bool>& %t%
    2009 %$$
    2010 and its size is $icode m$$.
    2011 The input values of its elements
    2012 are not specified (must not matter).
    2013 Upon return, $icode t$$ is a
    2014 sparsity pattern for
    2015 $latex T(x) \in B^{1 \times n}$$ where
    2016 $latex \[
    2017         T(x) = (g \circ f)^{(1)} (x) = S(x) * f^{(1)} (x)
    2018 \]$$
    2019 
    2020 $subhead q$$
    2021 The argument $icode q$$ has prototype
    2022 $codei%
    2023      size_t %q%
    2024 %$$
    2025 It specifies the number of columns in
    2026 $latex R \in B^{n \times q}$$,
    2027 $latex U(x) \in B^{m \times q}$$, and
    2028 $latex V(x) \in B^{n \times q}$$.
    2029 
    2030 $subhead r$$
    2031 This argument has prototype
    2032 $codei%
    2033      const %atomic_sparsity%& %r%
    2034 %$$
    2035 and is a $cref/atomic_sparsity/atomic_option/atomic_sparsity/$$ pattern for
    2036 $latex R \in B^{n \times q}$$.
    2037 
    2038 $head u$$
    2039 This argument has prototype
    2040 $codei%
    2041      const %atomic_sparsity%& %u%
    2042 %$$
    2043 and is a $cref/atomic_sparsity/atomic_option/atomic_sparsity/$$ pattern for
    2044 $latex U(x) \in B^{m \times q}$$ which is defined by
    2045 $latex \[
    2046 \begin{array}{rcl}
    2047 U(x)
    2048 & = &
    2049 \{ \partial_u \{ \partial_y g[ y + f^{(1)} (x) R u ] \}_{y=f(x)} \}_{u=0}
    2050 \\
    2051 & = &
    2052 \partial_u \{ g^{(1)} [ f(x) + f^{(1)} (x) R u ] \}_{u=0}
    2053 \\
    2054 & = &
    2055 g^{(2)} [ f(x) ] f^{(1)} (x) R
    2056 \end{array}
    2057 \] $$
    2058 
    2059 $subhead v$$
    2060 This argument has prototype
    2061 $codei%
    2062      %atomic_sparsity%& %v%
    2063 %$$
    2064 The input value of its elements
    2065 are not specified (must not matter).
    2066 Upon return, $icode v$$ is a
    2067 $cref/atomic_sparsity/atomic_option/atomic_sparsity/$$ pattern for
    2068 $latex V(x) \in B^{n \times q}$$ which is defined by
    2069 $latex \[
    2070 \begin{array}{rcl}
    2071 V(x)
    2072 & = &
    2073 \partial_u [ \partial_x (g \circ f) ( x + R u )  ]_{u=0}
    2074 \\
    2075 & = &
    2076 \partial_u [ (g \circ f)^{(1)}( x + R u )  ]_{u=0}
    2077 \\
    2078 & = &
    2079 (g \circ f)^{(2)}( x ) R
    2080 \\
    2081 & = &
    2082 f^{(1)} (x)^\R{T} g^{(2)} [ f(x) ] f^{(1)} (x)  R
    2083 +
    2084 \sum_{i=1}^m g_i^{(1)} [ f(x) ] \; f_i^{(2)} (x) R
    2085 \\
    2086 & = &
    2087 f^{(1)} (x)^\R{T} U(x)
    2088 +
    2089 \sum_{i=1}^m S_i (x) \; f_i^{(2)} (x) R
    2090 \end{array}
    2091 \] $$
    2092 
    2093 $subhead x$$
    2094 $index deprecated$$
    2095 The argument has prototype
    2096 $codei%
    2097         const CppAD::vector<%Base%>& %x%
    2098 %$$
    2099 and size is equal to the $icode n$$.
    2100 This is the $cref Value$$ value corresponding to the parameters in the
    2101 vector $cref/ax/atomic_afun/ax/$$ (when the atomic function was called).
    2102 To be specific, if
    2103 $codei%
    2104         if( Parameter(%ax%[%i%]) == true )
    2105                 %x%[%i%] = Value( %ax%[%i%] );
    2106         else
    2107                 %x%[%i%] = CppAD::numeric_limits<%Base%>::quiet_NaN();
    2108 %$$
    2109 The version of this function with out the $icode x$$ argument is deprecated;
    2110 i.e., you should include the argument even if you do not use it.
    2111 
    2112 $children%
    2113         example/atomic/rev_sparse_hes.cpp
    2114 %$$
    2115 $head Examples$$
    2116 The file $cref atomic_rev_sparse_hes.cpp$$ contains an example and test
    2117 that uses this routine.
    2118 It returns true if the test passes and false if it fails.
    2119 
    2120 $end
    2121 -----------------------------------------------------------------------------
    2122 */
    2123 /*!
    2124 Link from reverse Hessian sparsity sweep to base_atomic
    2125 
    2126 \param vx [in]
    2127 which componens of x are variables.
    2128 
    2129 \param s [in]
    2130 is the reverse Jacobian sparsity pattern w.r.t the result vector y.
    2131 
    2132 \param t [out]
    2133 is the reverse Jacobian sparsity pattern w.r.t the argument vector x.
    2134 
    2135 \param q [in]
    2136 is the column dimension for the sparsity partterns.
    2137 
    2138 \param r [in]
    2139 is the forward Jacobian sparsity pattern w.r.t the argument vector x
    2140 
    2141 \param u [in]
    2142 is the Hessian sparsity pattern w.r.t the result vector y.
    2143 
    2144 \param v [out]
    2145 is the Hessian sparsity pattern w.r.t the argument vector x.
    2146 
    2147 \param x [in]
    2148 is the integer value of the x arguments that are parameters.
    2149 */
    2150 virtual bool rev_sparse_hes(
    2151         const vector<bool>&                     vx ,
    2152         const vector<bool>&                     s  ,
    2153               vector<bool>&                     t  ,
    2154         size_t                                  q  ,
    2155         const vector< std::set<size_t> >&       r  ,
    2156         const vector< std::set<size_t> >&       u  ,
    2157               vector< std::set<size_t> >&       v  ,
    2158         const vector<Base>&                     x  )
    2159 {       return false; }
    2160 virtual bool rev_sparse_hes(
    2161         const vector<bool>&                     vx ,
    2162         const vector<bool>&                     s  ,
    2163               vector<bool>&                     t  ,
    2164         size_t                                  q  ,
    2165         const vector<bool>&                     r  ,
    2166         const vector<bool>&                     u  ,
    2167               vector<bool>&                     v  ,
    2168         const vector<Base>&                     x  )
    2169 {       return false; }
    2170 virtual bool rev_sparse_hes(
    2171         const vector<bool>&                     vx ,
    2172         const vector<bool>&                     s  ,
    2173               vector<bool>&                     t  ,
    2174         size_t                                  q  ,
    2175         const vectorBool&                       r  ,
    2176         const vectorBool&                       u  ,
    2177               vectorBool&                       v  ,
    2178         const vector<Base>&                     x  )
    2179 {       return false; }
    2180 // deprecated
    2181 virtual bool rev_sparse_hes(
    2182         const vector<bool>&                     vx ,
    2183         const vector<bool>&                     s  ,
    2184               vector<bool>&                     t  ,
    2185         size_t                                  q  ,
    2186         const vector< std::set<size_t> >&       r  ,
    2187         const vector< std::set<size_t> >&       u  ,
    2188               vector< std::set<size_t> >&       v  )
    2189 {       return false; }
    2190 virtual bool rev_sparse_hes(
    2191         const vector<bool>&                     vx ,
    2192         const vector<bool>&                     s  ,
    2193               vector<bool>&                     t  ,
    2194         size_t                                  q  ,
    2195         const vector<bool>&                     r  ,
    2196         const vector<bool>&                     u  ,
    2197               vector<bool>&                     v  )
    2198 {       return false; }
    2199 virtual bool rev_sparse_hes(
    2200         const vector<bool>&                     vx ,
    2201         const vector<bool>&                     s  ,
    2202               vector<bool>&                     t  ,
    2203         size_t                                  q  ,
    2204         const vectorBool&                       r  ,
    2205         const vectorBool&                       u  ,
    2206               vectorBool&                       v  )
    2207 {       return false; }
    2208 /*!
    2209 Link, before case split, from rev_hes_sweep to atomic_base.
    2210 
    2211 \tparam InternalSparsity
    2212 Is the used internaly for sparsity calculations; i.e.,
    2213 sparse_pack or sparse_list.
    2214 
    2215 \param x
    2216 is parameter arguments to the function, other components are nan.
    2217 
    2218 \param x_index
    2219 is the variable index, on the tape, for the arguments to this function.
    2220 This size of x_index is n, the number of arguments to this function.
    2221 
    2222 \param y_index
    2223 is the variable index, on the tape, for the results for this function.
    2224 This size of y_index is m, the number of results for this function.
    2225 
    2226 \param for_jac_sparsity
    2227 On input, for j = 0, ... , n-1, the sparsity pattern with index x_index[j],
    2228 is the forward Jacobian sparsity for the j-th argument to this atomic function.
    2229 
    2230 \param rev_jac_flag
    2231 This shows which variables affect the function we are
    2232 computing the Hessian of.
    2233 On input, for i = 0, ... , m-1, the rev_jac_flag[ y_index[i] ] is true
    2234 if the Jacobian of function (we are computing sparsity for) is no-zero.
    2235 Upon return, for j = 0, ... , n-1, rev_jac_flag [ x_index[j] ]
    2236 as been adjusted to accound removing this atomic function.
    2237 
    2238 \param rev_hes_sparsity
    2239 This is the sparsity pattern for the Hessian.
    2240 On input, for i = 0, ... , m-1, row y_index[i] is the reverse Hessian sparsity
    2241 with one of the partials with respect to to y_index[i].
    2242 */
    2243 template <class InternalSparsity>
    2244 void rev_sparse_hes(
    2245         const vector<Base>&        x                ,
    2246         const vector<size_t>&      x_index          ,
    2247         const vector<size_t>&      y_index          ,
    2248         const InternalSparsity&    for_jac_sparsity ,
    2249         bool*                      rev_jac_flag     ,
    2250         InternalSparsity&          rev_hes_sparsity )
    2251 {       CPPAD_ASSERT_UNKNOWN( for_jac_sparsity.end() == rev_hes_sparsity.end() );
    2252         size_t q           = rev_hes_sparsity.end();
    2253         size_t n           = x_index.size();
    2254         size_t m           = y_index.size();
    2255         bool   ok          = false;
    2256         size_t thread      = thread_alloc::thread_num();
    2257         allocate_work(thread);
    2258         bool   zero_empty  = true;
    2259         bool   input_empty = false;
    2260         bool   transpose   = false;
    2261         //
    2262         // vx
    2263         vector<bool> vx(n);
    2264         for(size_t j = 0; j < n; j++)
    2265                 vx[j] = x_index[j] != 0;
    2266         //
    2267         // note that s and t are vectors so transpose does not matter for bool case
    2268         vector<bool> bool_s( work_[thread]->bool_s );
    2269         vector<bool> bool_t( work_[thread]->bool_t );
    2270         //
    2271         bool_s.resize(m);
    2272         bool_t.resize(n);
    2273         //
    2274         for(size_t i = 0; i < m; i++)
    2275         {       if( y_index[i] > 0  )
    2276                         bool_s[i] = rev_jac_flag[ y_index[i] ];
    2277         }
    2278         //
    2279         std::string msg = ": atomic_base.rev_sparse_hes: returned false";
    2280         if( sparsity_ == pack_sparsity_enum )
    2281         {       vectorBool&  pack_r( work_[thread]->pack_r );
    2282                 vectorBool&  pack_u( work_[thread]->pack_u );
    2283                 vectorBool&  pack_v( work_[thread]->pack_h );
    2284                 //
    2285                 pack_v.resize(n * q);
    2286                 //
    2287                 local::get_internal_sparsity(
    2288                         transpose, x_index, for_jac_sparsity, pack_r
    2289                 );
    2290                 local::get_internal_sparsity(
    2291                         transpose, y_index, rev_hes_sparsity, pack_u
    2292                 );
    2293                 //
    2294                 ok = rev_sparse_hes(vx, bool_s, bool_t, q, pack_r, pack_u, pack_v, x);
    2295                 if( ! ok )
    2296                         ok = rev_sparse_hes(vx, bool_s, bool_t, q, pack_r, pack_u, pack_v);
    2297                 if( ! ok )
    2298                 {       msg = afun_name() + msg + " sparsity = pack_sparsity_enum";
    2299                         CPPAD_ASSERT_KNOWN(false, msg.c_str());
    2300                 }
    2301                 local::set_internal_sparsity(zero_empty, input_empty,
    2302                         transpose, x_index, rev_hes_sparsity, pack_v
    2303                 );
    2304         }
    2305         else if( sparsity_ == bool_sparsity_enum )
    2306         {       vector<bool>&  bool_r( work_[thread]->bool_r );
    2307                 vector<bool>&  bool_u( work_[thread]->bool_u );
    2308                 vector<bool>&  bool_v( work_[thread]->bool_h );
    2309                 //
    2310                 bool_v.resize(n * q);
    2311                 //
    2312                 local::get_internal_sparsity(
    2313                         transpose, x_index, for_jac_sparsity, bool_r
    2314                 );
    2315                 local::get_internal_sparsity(
    2316                         transpose, y_index, rev_hes_sparsity, bool_u
    2317                 );
    2318                 //
    2319                 ok = rev_sparse_hes(vx, bool_s, bool_t, q, bool_r, bool_u, bool_v, x);
    2320                 if( ! ok )
    2321                         ok = rev_sparse_hes(vx, bool_s, bool_t, q, bool_r, bool_u, bool_v);
    2322                 if( ! ok )
    2323                 {       msg = afun_name() + msg + " sparsity = bool_sparsity_enum";
    2324                         CPPAD_ASSERT_KNOWN(false, msg.c_str());
    2325                 }
    2326                 local::set_internal_sparsity(zero_empty, input_empty,
    2327                         transpose, x_index, rev_hes_sparsity, bool_v
    2328                 );
    2329         }
    2330         else
    2331         {       CPPAD_ASSERT_UNKNOWN( sparsity_ == set_sparsity_enum );
    2332                 vector< std::set<size_t> >&  set_r( work_[thread]->set_r );
    2333                 vector< std::set<size_t> >&  set_u( work_[thread]->set_u );
    2334                 vector< std::set<size_t> >&  set_v( work_[thread]->set_h );
    2335                 //
    2336                 set_v.resize(n);
    2337                 //
    2338                 local::get_internal_sparsity(
    2339                         transpose, x_index, for_jac_sparsity, set_r
    2340                 );
    2341                 local::get_internal_sparsity(
    2342                         transpose, y_index, rev_hes_sparsity, set_u
    2343                 );
    2344                 //
    2345                 ok = rev_sparse_hes(vx, bool_s, bool_t, q, set_r, set_u, set_v, x);
    2346                 if( ! ok )
    2347                         ok = rev_sparse_hes(vx, bool_s, bool_t, q, set_r, set_u, set_v);
    2348                 if( ! ok )
    2349                 {       msg = afun_name() + msg + " sparsity = set_sparsity_enum";
    2350                         CPPAD_ASSERT_KNOWN(false, msg.c_str());
    2351                 }
    2352                 local::set_internal_sparsity(zero_empty, input_empty,
    2353                         transpose, x_index, rev_hes_sparsity, set_v
    2354                 );
    2355         }
    2356         for(size_t j = 0; j < n; j++)
    2357         {       if( x_index[j] > 0  )
    2358                         rev_jac_flag[ x_index[j] ] |= bool_t[j];
    2359         }
    2360         return;
    2361 }
    2362 /*
    2363 ------------------------------------------------------------------------------
    2364 $begin atomic_base_clear$$
    2365 $spell
    2366         sq
    2367         alloc
    2368 $$
    2369 
    2370 $section Free Static Variables$$
    2371 $mindex clear$$
    2372 
    2373 $head Syntax$$
    2374 $codei%atomic_base<%Base%>::clear()%$$
    2375 
    2376 $head Purpose$$
    2377 Each $code atomic_base$$ objects holds onto work space in order to
    2378 avoid repeated memory allocation calls and thereby increase speed
    2379 (until it is deleted).
    2380 If an the $code atomic_base$$ object is global or static because,
    2381 the it does not get deleted.
    2382 This is a problem when using
    2383 $code thread_alloc$$ $cref/free_all/ta_free_all/$$
    2384 to check that all allocated memory has been freed.
    2385 Calling this $code clear$$ function will free all the
    2386 memory currently being held onto by the
    2387 $codei%atomic_base<%Base%>%$$ class.
    2388 
    2389 $head Future Use$$
    2390 If there is future use of an $code atomic_base$$ object,
    2391 after a call to $code clear$$,
    2392 the work space will be reallocated and held onto.
    2393 
    2394 $head Restriction$$
    2395 This routine cannot be called
    2396 while in $cref/parallel/ta_in_parallel/$$ execution mode.
    2397 
    2398 $end
    2399 ------------------------------------------------------------------------------
    2400 */
    2401 /*!
    2402 Free all thread_alloc static memory held by atomic_base (avoids reallocations).
    2403 (This does not include class_object() which is an std::vector.)
    2404 */
    2405 /// Free vector memory used by this class (work space)
    2406 static void clear(void)
    2407 {       CPPAD_ASSERT_KNOWN(
    2408                 ! thread_alloc::in_parallel() ,
    2409                 "cannot use atomic_base clear during parallel execution"
    2410         );
    2411         size_t i = class_object().size();
    2412         while(i--)
    2413         {       atomic_base* op = class_object()[i];
    2414                 if( op != CPPAD_NULL )
    2415                 {       for(size_t thread = 0; thread < CPPAD_MAX_NUM_THREADS; thread++)
    2416                                 op->free_work(thread);
    2417                 }
    2418         }
    2419         return;
    2420 }
    2421 // -------------------------------------------------------------------------
    2422 /*!
    2423 Set value of id (used by deprecated old_atomic class)
    2424 
    2425 This function is called just before calling any of the virtual function
    2426 and has the corresponding id of the corresponding virtual call.
    2427 */
    2428 virtual void set_old(size_t id)
    2429 { }
     559                return;
     560        }
     561        /// atomic_base function object corresponding to a certain index
     562        static atomic_base* class_object(size_t index)
     563        {       CPPAD_ASSERT_UNKNOWN( class_object().size() > index );
     564                return class_object()[index];
     565        }
     566        /// atomic_base function name corresponding to a certain index
     567        static const std::string& class_name(size_t index)
     568        {       CPPAD_ASSERT_UNKNOWN( class_name().size() > index );
     569                return class_name()[index];
     570        }
     571
     572        /*!
     573        Set value of id (used by deprecated old_atomic class)
     574
     575        This function is called just before calling any of the virtual function
     576        and has the corresponding id of the corresponding virtual call.
     577        */
     578        virtual void set_old(size_t id)
     579        { }
    2430580// ---------------------------------------------------------------------------
    2431581};
    2432582} // END_CPPAD_NAMESPACE
     583
     584// functitons implemented in cppad/core/atomic_base files
     585# include <cppad/core/atomic_base/ctor.hpp>
     586# include <cppad/core/atomic_base/option.hpp>
     587# include <cppad/core/atomic_base/afun.hpp>
     588# include <cppad/core/atomic_base/forward.hpp>
     589# include <cppad/core/atomic_base/reverse.hpp>
     590# include <cppad/core/atomic_base/for_sparse_jac.hpp>
     591# include <cppad/core/atomic_base/rev_sparse_jac.hpp>
     592# include <cppad/core/atomic_base/for_sparse_hes.hpp>
     593# include <cppad/core/atomic_base/rev_sparse_hes.hpp>
     594# include <cppad/core/atomic_base/clear.hpp>
     595
    2433596# endif
  • trunk/cppad/core/capacity_order.hpp

    r4010 r4024  
    153153*/
    154154
    155 template <typename Base>
    156 void ADFun<Base>::capacity_order(size_t c, size_t r)
     155template <typename Base, typename RecBase>
     156void ADFun<Base,RecBase>::capacity_order(size_t c, size_t r)
    157157{       // temporary indices
    158158        size_t i, k, ell;
     
    238238*/
    239239
    240 template <typename Base>
    241 void ADFun<Base>::capacity_order(size_t c)
     240template <typename Base, typename RecBase>
     241void ADFun<Base,RecBase>::capacity_order(size_t c)
    242242{       size_t r;
    243243        if( (c == 0) | (c == 1) )
  • trunk/cppad/core/check_for_nan.hpp

    r4023 r4024  
    156156new value for this flag.
    157157*/
    158 template <class Base>
    159 void ADFun<Base>::check_for_nan(bool value)
     158template <class Base, class RecBase>
     159void ADFun<Base,RecBase>::check_for_nan(bool value)
    160160{       check_for_nan_ = value; }
    161161
     
    166166current value of check_for_nan_.
    167167*/
    168 template <class Base>
    169 bool ADFun<Base>::check_for_nan(void) const
     168template <class Base, class RecBase>
     169bool ADFun<Base,RecBase>::check_for_nan(void) const
    170170{       return check_for_nan_; }
    171171
  • trunk/cppad/core/checkpoint.hpp

    r4023 r4024  
    2424$begin checkpoint$$
    2525$spell
     26        alloc
     27        inuse
    2628        sv
    2729        var
     
    5153
    5254$head See Also$$
    53 $cref reverse_checkpoint.cpp$$
     55$cref atomic_base$$, $cref reverse_checkpoint.cpp$$
    5456
    5557$head Purpose$$
     
    219221but this is not part of the CppAD specifications and my change.)
    220222
    221 $head clear$$
    222 The $code atomic_base$$ class holds onto static work space in order to
    223 increase speed by avoiding system memory allocation calls.
     223$head Memory$$
     224
     225$subhead Restriction$$
     226The $code clear$$ routine cannot be called
     227while in $cref/parallel/ta_in_parallel/$$ execution mode.
     228
     229$head Parallel Mode$$
     230The CppAD checkpoint function delays the caching of certain calculations
     231until they are needed.
     232In $cref/parallel model/ta_parallel_setup/$$,
     233this may result in $cref/thread_alloc::inuse(thread)/ta_inuse/$$
     234being non-zero even though the specified thread is no longer active.
     235This memory will be freed when the checkpoint object is deleted.
     236
     237$subhead clear$$
     238The $code atomic_base$$ class holds onto static work space
     239that is not connected to a particular object
     240(in order to increase speed by avoiding system memory allocation calls).
    224241This call makes to work space $cref/available/ta_available/$$ to
    225242for other uses by the same thread.
     
    227244user atomic functions for a specific value of $icode Base$$.
    228245
    229 $subhead Restriction$$
    230 The $code clear$$ routine cannot be called
    231 while in $cref/parallel/ta_in_parallel/$$ execution mode.
    232 
    233246$children%example/atomic/checkpoint.cpp
    234247        %example/atomic/mul_level.cpp
     
    243256$end
    244257*/
    245 
    246 # define CPPAD_TMP_MULTI_THREAD 0
    247 # ifdef _OPENMP
    248 # ifdef CPPAD_FOR_TMB
    249 # undef  CPPAD_TMP_MULTI_THREAD
    250 # define CPPAD_TMP_MULTI_THREAD 1
    251 # endif
    252 # endif
    253 
    254 // ============================================================================
    255 # if ! CPPAD_TMP_MULTI_THREAD
    256 // ============================================================================
    257258
    258259template <class Base>
     
    263264        typedef typename atomic_base<Base>::option_enum option_enum;
    264265        //
    265         /// AD function corresponding to this checkpoint object
    266         ADFun<Base> f_;
     266        // ------------------------------------------------------------------------
     267        // member_
     268        // ------------------------------------------------------------------------
     269        // same checkpoint object can be used by multiple threads
     270        struct member_struct {
     271                //
     272                /// AD function corresponding to this checkpoint object
     273                ADFun<Base>             f_;
     274                ADFun< AD<Base>, Base > af_;
     275                //
     276                /// sparsity for entire Jacobian f(x)^{(1)}
     277                /// does not change so can cache it
     278                local::sparse_list         jac_sparse_set_;
     279                vectorBool                 jac_sparse_bool_;
     280                //
     281                /// sparsity for sum_i f_i(x)^{(2)} does not change so can cache it
     282                local::sparse_list         hes_sparse_set_;
     283                vectorBool                 hes_sparse_bool_;
     284        };
     285        /// use pointers and allocate memory to avoid false sharing
     286        member_struct* member_[CPPAD_MAX_NUM_THREADS];
    267287        //
    268         /// sparsity for entire Jacobian f(x)^{(1)} does not change so can cache it
    269         local::sparse_list         jac_sparse_set_;
    270         vectorBool                 jac_sparse_bool_;
     288        /// allocate member_ for this thread
     289        void allocate_member(size_t thread)
     290        {       // thread zero is the master thread
     291                size_t master = 0;
     292                //
     293                if( member_[thread] == CPPAD_NULL )
     294                {       member_[thread] = new member_struct;
     295                        // The function is only recorded by the master thread,
     296                        // other threads have copy.
     297                        if( thread != master )
     298                                member_[thread]->f_ = member_[master]->f_;
     299                }
     300                return;
     301        }
    271302        //
    272         /// sparsity for sum_i f_i(x)^{(2)} does not change so can cache it
    273         local::sparse_list         hes_sparse_set_;
    274         vectorBool                 hes_sparse_bool_;
     303        /// free member_ for this thread
     304        void free_member(size_t thread)
     305        {       if( member_[thread] != CPPAD_NULL )
     306                {       delete member_[thread];
     307                        member_[thread] = CPPAD_NULL;
     308                }
     309                return;
     310        }
    275311        // ------------------------------------------------------------------------
    276312        option_enum sparsity(void)
     
    278314        // ------------------------------------------------------------------------
    279315        /// set jac_sparse_set_
    280         void set_jac_sparse_set(void)
    281         {       CPPAD_ASSERT_UNKNOWN( jac_sparse_set_.n_set() == 0 );
    282                 bool transpose  = false;
    283                 bool dependency = true;
    284                 size_t n = f_.Domain();
    285                 size_t m = f_.Range();
    286                 // Use the choice for forward / reverse that results in smaller
    287                 // size for the sparsity pattern of all variables in the tape.
    288                 if( n <= m )
    289                 {       local::sparse_list identity;
    290                         identity.resize(n, n);
    291                         for(size_t j = 0; j < n; j++)
    292                         {       // use add_element because only adding one element per set
    293                                 identity.add_element(j, j);
    294                         }
    295                         f_.ForSparseJacCheckpoint(
    296                                 n, identity, transpose, dependency, jac_sparse_set_
    297                         );
    298                         f_.size_forward_set(0);
    299                 }
    300                 else
    301                 {       local::sparse_list identity;
    302                         identity.resize(m, m);
    303                         for(size_t i = 0; i < m; i++)
    304                         {       // use add_element because only adding one element per set
    305                                 identity.add_element(i, i);
    306                         }
    307                         f_.RevSparseJacCheckpoint(
    308                                 m, identity, transpose, dependency, jac_sparse_set_
    309                         );
    310                 }
    311                 CPPAD_ASSERT_UNKNOWN( f_.size_forward_set() == 0 );
    312                 CPPAD_ASSERT_UNKNOWN( f_.size_forward_bool() == 0 );
    313         }
     316        void set_jac_sparse_set(void);
    314317        /// set jac_sparse_bool_
    315         void set_jac_sparse_bool(void)
    316         {       CPPAD_ASSERT_UNKNOWN( jac_sparse_bool_.size() == 0 );
    317                 bool transpose  = false;
    318                 bool dependency = true;
    319                 size_t n = f_.Domain();
    320                 size_t m = f_.Range();
    321                 // Use the choice for forward / reverse that results in smaller
    322                 // size for the sparsity pattern of all variables in the tape.
    323                 if( n <= m )
    324                 {       vectorBool identity(n * n);
    325                         for(size_t j = 0; j < n; j++)
    326                         {       for(size_t i = 0; i < n; i++)
    327                                         identity[ i * n + j ] = (i == j);
    328                         }
    329                         jac_sparse_bool_ = f_.ForSparseJac(
    330                                 n, identity, transpose, dependency
    331                         );
    332                         f_.size_forward_bool(0);
    333                 }
    334                 else
    335                 {       vectorBool identity(m * m);
    336                         for(size_t j = 0; j < m; j++)
    337                         {       for(size_t i = 0; i < m; i++)
    338                                         identity[ i * m + j ] = (i == j);
    339                         }
    340                         jac_sparse_bool_ = f_.RevSparseJac(
    341                                 m, identity, transpose, dependency
    342                         );
    343                 }
    344                 CPPAD_ASSERT_UNKNOWN( f_.size_forward_bool() == 0 );
    345                 CPPAD_ASSERT_UNKNOWN( f_.size_forward_set() == 0 );
    346         }
     318        void set_jac_sparse_bool(void);
    347319        // ------------------------------------------------------------------------
    348320        /// set hes_sparse_set_
    349         void set_hes_sparse_set(void)
    350         {       CPPAD_ASSERT_UNKNOWN( hes_sparse_set_.n_set() == 0 );
    351                 size_t n = f_.Domain();
    352                 size_t m = f_.Range();
    353                 //
    354                 // set version of sparsity for vector of all ones
    355                 vector<bool> all_one(m);
    356                 for(size_t i = 0; i < m; i++)
    357                         all_one[i] = true;
    358 
    359                 // set version of sparsity for n by n idendity matrix
    360                 local::sparse_list identity;
    361                 identity.resize(n, n);
    362                 for(size_t j = 0; j < n; j++)
    363                 {       // use add_element because only adding one element per set
    364                         identity.add_element(j, j);
    365                 }
    366 
    367                 // compute sparsity pattern for H(x) = sum_i f_i(x)^{(2)}
    368                 bool transpose  = false;
    369                 bool dependency = false;
    370                 f_.ForSparseJacCheckpoint(
    371                         n, identity, transpose, dependency, jac_sparse_set_
    372                 );
    373                 f_.RevSparseHesCheckpoint(
    374                         n, all_one, transpose, hes_sparse_set_
    375                 );
    376                 CPPAD_ASSERT_UNKNOWN( hes_sparse_set_.n_set() == n );
    377                 CPPAD_ASSERT_UNKNOWN( hes_sparse_set_.end()   == n );
    378                 //
    379                 // drop the forward sparsity results from f_
    380                 f_.size_forward_set(0);
    381         }
     321        void set_hes_sparse_set(void);
    382322        /// set hes_sparse_bool_
    383         void set_hes_sparse_bool(void)
    384         {       CPPAD_ASSERT_UNKNOWN( hes_sparse_bool_.size() == 0 );
    385                 size_t n = f_.Domain();
    386                 size_t m = f_.Range();
    387                 //
    388                 // set version of sparsity for vector of all ones
    389                 vectorBool all_one(m);
    390                 for(size_t i = 0; i < m; i++)
    391                         all_one[i] = true;
    392 
    393                 // set version of sparsity for n by n idendity matrix
    394                 vectorBool identity(n * n);
    395                 for(size_t j = 0; j < n; j++)
    396                 {       for(size_t i = 0; i < n; i++)
    397                                 identity[ i * n + j ] = (i == j);
    398                 }
    399 
    400                 // compute sparsity pattern for H(x) = sum_i f_i(x)^{(2)}
    401                 bool transpose  = false;
    402                 bool dependency = false;
    403                 f_.ForSparseJac(n, identity, transpose, dependency);
    404                 hes_sparse_bool_ = f_.RevSparseHes(n, all_one, transpose);
    405                 CPPAD_ASSERT_UNKNOWN( hes_sparse_bool_.size() == n * n );
    406                 //
    407                 // drop the forward sparsity results from f_
    408                 f_.size_forward_bool(0);
    409                 CPPAD_ASSERT_UNKNOWN( f_.size_forward_bool() == 0 );
    410                 CPPAD_ASSERT_UNKNOWN( f_.size_forward_set() == 0 );
    411         }
     323        void set_hes_sparse_bool(void);
    412324        // ------------------------------------------------------------------------
    413325        /*!
     
    421333                const sparsity_type&                    r  ,
    422334                      sparsity_type&                    s  ,
    423                 const vector<Base>&                     x  )
    424         {       // during user sparsity calculations
    425                 size_t m = f_.Range();
    426                 size_t n = f_.Domain();
    427                 if( jac_sparse_bool_.size() == 0 )
    428                         set_jac_sparse_bool();
    429                 if( jac_sparse_set_.n_set() != 0 )
    430                         jac_sparse_set_.resize(0, 0);
    431                 CPPAD_ASSERT_UNKNOWN( jac_sparse_bool_.size() == m * n );
    432                 CPPAD_ASSERT_UNKNOWN( jac_sparse_set_.n_set() == 0 );
    433                 CPPAD_ASSERT_UNKNOWN( r.size() == n * q );
    434                 CPPAD_ASSERT_UNKNOWN( s.size() == m * q );
    435                 //
    436                 bool ok = true;
    437                 for(size_t i = 0; i < m; i++)
    438                 {       for(size_t k = 0; k < q; k++)
    439                                 s[i * q + k] = false;
    440                 }
    441                 // sparsity for  s = jac_sparse_bool_ * r
    442                 for(size_t i = 0; i < m; i++)
    443                 {       for(size_t k = 0; k < q; k++)
    444                         {       // initialize sparsity for S(i,k)
    445                                 bool s_ik = false;
    446                                 // S(i,k) = sum_j J(i,j) * R(j,k)
    447                                 for(size_t j = 0; j < n; j++)
    448                                 {       bool J_ij = jac_sparse_bool_[ i * n + j];
    449                                         bool R_jk = r[j * q + k ];
    450                                         s_ik |= ( J_ij & R_jk );
    451                                 }
    452                                 s[i * q + k] = s_ik;
    453                         }
    454                 }
    455                 return ok;
    456         }
     335                const vector<Base>&                     x
     336        );
    457337        // ------------------------------------------------------------------------
    458338        /*!
     
    466346                const sparsity_type&                    rt ,
    467347                      sparsity_type&                    st ,
    468                 const vector<Base>&                     x  )
    469         {       // during user sparsity calculations
    470                 size_t m = f_.Range();
    471                 size_t n = f_.Domain();
    472                 if( jac_sparse_bool_.size() == 0 )
    473                         set_jac_sparse_bool();
    474                 if( jac_sparse_set_.n_set() != 0 )
    475                         jac_sparse_set_.resize(0, 0);
    476                 CPPAD_ASSERT_UNKNOWN( jac_sparse_bool_.size() == m * n );
    477                 CPPAD_ASSERT_UNKNOWN( jac_sparse_set_.n_set() == 0 );
    478                 CPPAD_ASSERT_UNKNOWN( rt.size() == m * q );
    479                 CPPAD_ASSERT_UNKNOWN( st.size() == n * q );
    480                 bool ok  = true;
    481                 //
    482                 // S = R * J where J is jacobian
    483                 for(size_t i = 0; i < q; i++)
    484                 {       for(size_t j = 0; j < n; j++)
    485                         {       // initialize sparsity for S(i,j)
    486                                 bool s_ij = false;
    487                                 // S(i,j) = sum_k R(i,k) * J(k,j)
    488                                 for(size_t k = 0; k < m; k++)
    489                                 {       // sparsity for R(i, k)
    490                                         bool R_ik = rt[ k * q + i ];
    491                                         bool J_kj = jac_sparse_bool_[ k * n + j ];
    492                                         s_ij     |= (R_ik & J_kj);
    493                                 }
    494                                 // set sparsity for S^T
    495                                 st[ j * q + i ] = s_ij;
    496                         }
    497                 }
    498                 return ok;
    499         }
     348                const vector<Base>&                     x
     349        );
    500350        /*!
    501351        Link from user_atomic to reverse sparse Hessian  bools
     
    512362                const sparsity_type&                    u  ,
    513363                      sparsity_type&                    v  ,
    514                 const vector<Base>&                     x  )
    515         {       size_t n = f_.Domain();
    516 # ifndef NDEBUG
    517                 size_t m = f_.Range();
    518 # endif
    519                 CPPAD_ASSERT_UNKNOWN( vx.size() == n );
    520                 CPPAD_ASSERT_UNKNOWN(  s.size() == m );
    521                 CPPAD_ASSERT_UNKNOWN(  t.size() == n );
    522                 CPPAD_ASSERT_UNKNOWN(  r.size() == n * q );
    523                 CPPAD_ASSERT_UNKNOWN(  u.size() == m * q );
    524                 CPPAD_ASSERT_UNKNOWN(  v.size() == n * q );
    525                 //
    526                 bool ok        = true;
    527 
    528                 // make sure hes_sparse_bool_ has been set
    529                 if( hes_sparse_bool_.size() == 0 )
    530                         set_hes_sparse_bool();
    531                 if( hes_sparse_set_.n_set() != 0 )
    532                         hes_sparse_set_.resize(0, 0);
    533                 CPPAD_ASSERT_UNKNOWN( hes_sparse_bool_.size() == n * n );
    534                 CPPAD_ASSERT_UNKNOWN( hes_sparse_set_.n_set() == 0 );
    535 
    536 
    537                 // compute sparsity pattern for T(x) = S(x) * f'(x)
    538                 t = f_.RevSparseJac(1, s);
    539 # ifndef NDEBUG
    540                 for(size_t j = 0; j < n; j++)
    541                         CPPAD_ASSERT_UNKNOWN( vx[j] || ! t[j] )
    542 # endif
    543 
    544                 // V(x) = f'(x)^T * g''(y) * f'(x) * R  +  g'(y) * f''(x) * R
    545                 // U(x) = g''(y) * f'(x) * R
    546                 // S(x) = g'(y)
    547 
    548                 // compute sparsity pattern for A(x) = f'(x)^T * U(x)
    549                 bool transpose = true;
    550                 sparsity_type a(n * q);
    551                 a = f_.RevSparseJac(q, u, transpose);
    552 
    553                 // Need sparsity pattern for H(x) = (S(x) * f(x))''(x) * R,
    554                 // but use less efficient sparsity for  f(x)''(x) * R so that
    555                 // hes_sparse_set_ can be used every time this is needed.
    556                 for(size_t i = 0; i < n; i++)
    557                 {       for(size_t k = 0; k < q; k++)
    558                         {       // initialize sparsity pattern for H(i,k)
    559                                 bool h_ik = false;
    560                                 // H(i,k) = sum_j f''(i,j) * R(j,k)
    561                                 for(size_t j = 0; j < n; j++)
    562                                 {       bool f_ij = hes_sparse_bool_[i * n + j];
    563                                         bool r_jk = r[j * q + k];
    564                                         h_ik     |= ( f_ij & r_jk );
    565                                 }
    566                                 // sparsity for H(i,k)
    567                                 v[i * q + k] = h_ik;
    568                         }
    569                 }
    570 
    571                 // compute sparsity pattern for V(x) = A(x) + H(x)
    572                 for(size_t i = 0; i < n; i++)
    573                 {       for(size_t k = 0; k < q; k++)
    574                                 // v[ i * q + k ] |= a[ i * q + k];
    575                                 v[ i * q + k ] = bool(v[ i * q + k]) | bool(a[ i * q + k]);
    576                 }
    577                 return ok;
    578         }
     364                const vector<Base>&                     x
     365        );
    579366public:
    580367        // ------------------------------------------------------------------------
     
    614401                                atomic_base<Base>::pack_sparsity_enum  ,
    615402                bool                           optimize = true
    616         ) : atomic_base<Base>(name, sparsity)
    617         {       CheckSimpleVector< CppAD::AD<Base> , ADVector>();
    618 
    619                 // make a copy of ax because Independent modifies AD information
    620                 ADVector x_tmp(ax);
    621                 // delcare x_tmp as the independent variables
    622                 Independent(x_tmp);
    623                 // record mapping from x_tmp to ay
    624                 algo(x_tmp, ay);
    625                 // create function f_ : x -> y
    626                 f_.Dependent(ay);
    627                 if( optimize )
    628                 {       // suppress checking for nan in f_ results
    629                         // (see optimize documentation for atomic functions)
    630                         f_.check_for_nan(false);
    631                         //
    632                         // now optimize
    633                         f_.optimize();
     403        );
     404        /// destructor
     405        ~checkpoint(void)
     406        {
     407# ifndef NDEBUG
     408                if( thread_alloc::in_parallel() )
     409                {       std::string msg = atomic_base<Base>::afun_name();
     410                        msg += ": checkpoint destructor called in parallel mode.";
     411                        CPPAD_ASSERT_KNOWN(false, msg.c_str() );
    634412                }
    635                 // now disable checking of comparison operations
    636                 // 2DO: add a debugging mode that checks for changes and aborts
    637                 f_.compare_change_count(0);
     413# endif
     414                for(size_t thread = 0; thread < CPPAD_MAX_NUM_THREADS; ++thread)
     415                        free_member(thread);
    638416        }
    639417        // ------------------------------------------------------------------------
    640418        /*!
    641         Implement the user call to <tt>atom_fun.size_var()</tt>.
     419        Implement the user call to atom_fun.size_var().
    642420        */
    643421        size_t size_var(void)
    644         {       return f_.size_var(); }
    645         // ------------------------------------------------------------------------
    646         /*!
    647         Implement the user call to <tt>atom_fun(ax, ay)</tt>.
     422        {   // make sure member_ is allocated for this thread
     423                size_t thread = thread_alloc::thread_num();
     424                allocate_member(thread);
     425                //
     426                return member_[thread]->f_.size_var();
     427        }
     428        // ------------------------------------------------------------------------
     429        /*!
     430        Implement the user call to atom_fun(ax, ay).
    648431
    649432        \tparam ADVector
    650         A simple vector class with elements of type <code>AD<Base></code>.
     433        A simple vector class with elements of type AD<Base>.
    651434
    652435        \param id
     
    655438        \param ax
    656439        is the argument vector for this call,
    657         <tt>ax.size()</tt> determines the number of arguments.
     440        ax.size() determines the number of arguments.
    658441
    659442        \param ay
    660443        is the result vector for this call,
    661         <tt>ay.size()</tt> determines the number of results.
     444        ay.size() determines the number of results.
    662445        */
    663446        template <class ADVector>
     
    676459        */
    677460        virtual bool forward(
    678                 size_t                    p ,
    679                 size_t                    q ,
    680                 const vector<bool>&      vx ,
    681                       vector<bool>&      vy ,
    682                 const vector<Base>&      tx ,
    683                       vector<Base>&      ty )
    684         {       size_t n = f_.Domain();
    685                 size_t m = f_.Range();
    686                 //
    687                 CPPAD_ASSERT_UNKNOWN( f_.size_var() > 0 );
    688                 CPPAD_ASSERT_UNKNOWN( tx.size() % (q+1) == 0 );
    689                 CPPAD_ASSERT_UNKNOWN( ty.size() % (q+1) == 0 );
    690                 CPPAD_ASSERT_UNKNOWN( n == tx.size() / (q+1) );
    691                 CPPAD_ASSERT_UNKNOWN( m == ty.size() / (q+1) );
    692                 bool ok  = true;
    693                 //
    694                 if( vx.size() == 0 )
    695                 {       // during user forward mode
    696                         if( jac_sparse_set_.n_set() != 0 )
    697                                 jac_sparse_set_.resize(0,0);
    698                         if( jac_sparse_bool_.size() != 0 )
    699                                 jac_sparse_bool_.clear();
    700                         //
    701                         if( hes_sparse_set_.n_set() != 0 )
    702                                 hes_sparse_set_.resize(0,0);
    703                         if( hes_sparse_bool_.size() != 0 )
    704                                 hes_sparse_bool_.clear();
    705                 }
    706                 if( vx.size() > 0 )
    707                 {       // need Jacobian sparsity pattern to determine variable relation
    708                         // during user recording using checkpoint functions
    709                         if( sparsity() == atomic_base<Base>::set_sparsity_enum )
    710                         {       if( jac_sparse_set_.n_set() == 0 )
    711                                         set_jac_sparse_set();
    712                                 CPPAD_ASSERT_UNKNOWN( jac_sparse_set_.n_set() == m );
    713                                 CPPAD_ASSERT_UNKNOWN( jac_sparse_set_.end()   == n );
    714                                 //
    715                                 for(size_t i = 0; i < m; i++)
    716                                 {       vy[i] = false;
    717                                         local::sparse_list::const_iterator set_itr(
    718                                                 jac_sparse_set_, i
    719                                         );
    720                                         size_t j = *set_itr;
    721                                         while(j < n )
    722                                         {       // y[i] depends on the value of x[j]
    723                                                 // cast avoid Microsoft warning (should not be needed)
    724                                                 vy[i] |= static_cast<bool>( vx[j] );
    725                                                 j = *(++set_itr);
    726                                         }
    727                                 }
    728                         }
    729                         else
    730                         {       if( jac_sparse_set_.n_set() != 0 )
    731                                         jac_sparse_set_.resize(0, 0);
    732                                 if( jac_sparse_bool_.size() == 0 )
    733                                         set_jac_sparse_bool();
    734                                 CPPAD_ASSERT_UNKNOWN( jac_sparse_set_.n_set() == 0 );
    735                                 CPPAD_ASSERT_UNKNOWN( jac_sparse_bool_.size() == m * n );
    736                                 //
    737                                 for(size_t i = 0; i < m; i++)
    738                                 {       vy[i] = false;
    739                                         for(size_t j = 0; j < n; j++)
    740                                         {       if( jac_sparse_bool_[ i * n + j ] )
    741                                                 {       // y[i] depends on the value of x[j]
    742                                                         // cast avoid Microsoft warning
    743                                                         vy[i] |= static_cast<bool>( vx[j] );
    744                                                 }
    745                                         }
    746                                 }
    747                         }
    748                 }
    749                 // compute forward results for orders zero through q
    750                 ty = f_.Forward(q, tx);
    751 
    752                 // no longer need the Taylor coefficients in f_
    753                 // (have to reconstruct them every time)
    754                 // Hold onto sparsity pattern because it is always good.
    755                 size_t c = 0;
    756                 size_t r = 0;
    757                 f_.capacity_order(c, r);
    758                 return ok;
    759         }
     461                size_t                      p  ,
     462                size_t                      q  ,
     463                const vector<bool>&         vx ,
     464                      vector<bool>&         vy ,
     465                const vector<Base>&         tx ,
     466                      vector<Base>&         ty
     467        );
     468        virtual bool forward(
     469                size_t                      p   ,
     470                size_t                      q   ,
     471                const vector<bool>&         vx  ,
     472                      vector<bool>&         vy  ,
     473                const vector< AD<Base> >&   atx ,
     474                      vector< AD<Base> >&   aty
     475        );
    760476        // ------------------------------------------------------------------------
    761477        /*!
     
    765481        */
    766482        virtual bool reverse(
    767                 size_t                    q  ,
    768                 const vector<Base>&       tx ,
    769                 const vector<Base>&       ty ,
    770                       vector<Base>&       px ,
    771                 const vector<Base>&       py )
    772         {
    773 # ifndef NDEBUG
    774                 size_t n = f_.Domain();
    775                 size_t m = f_.Range();
    776 # endif
    777                 CPPAD_ASSERT_UNKNOWN( n == tx.size() / (q+1) );
    778                 CPPAD_ASSERT_UNKNOWN( m == ty.size() / (q+1) );
    779                 CPPAD_ASSERT_UNKNOWN( f_.size_var() > 0 );
    780                 CPPAD_ASSERT_UNKNOWN( tx.size() % (q+1) == 0 );
    781                 CPPAD_ASSERT_UNKNOWN( ty.size() % (q+1) == 0 );
    782                 bool ok  = true;
    783 
    784                 // put proper forward mode coefficients in f_
    785 # ifdef NDEBUG
    786                 // compute forward results for orders zero through q
    787                 f_.Forward(q, tx);
    788 # else
    789                 CPPAD_ASSERT_UNKNOWN( px.size() == n * (q+1) );
    790                 CPPAD_ASSERT_UNKNOWN( py.size() == m * (q+1) );
    791                 size_t i, j, k;
    792                 //
    793                 // compute forward results for orders zero through q
    794                 vector<Base> check_ty = f_.Forward(q, tx);
    795                 for(i = 0; i < m; i++)
    796                 {       for(k = 0; k <= q; k++)
    797                         {       j = i * (q+1) + k;
    798                                 CPPAD_ASSERT_UNKNOWN( check_ty[j] == ty[j] );
    799                         }
    800                 }
    801 # endif
    802                 // now can run reverse mode
    803                 px = f_.Reverse(q+1, py);
    804 
    805                 // no longer need the Taylor coefficients in f_
    806                 // (have to reconstruct them every time)
    807                 size_t c = 0;
    808                 size_t r = 0;
    809                 f_.capacity_order(c, r);
    810                 return ok;
    811         }
     483                size_t                          q  ,
     484                const vector<Base>&             tx ,
     485                const vector<Base>&             ty ,
     486                      vector<Base>&             px ,
     487                const vector<Base>&             py
     488        );
     489        virtual bool reverse(
     490                size_t                          q  ,
     491                const vector< AD<Base> >&       atx ,
     492                const vector< AD<Base> >&       aty ,
     493                      vector< AD<Base> >&       apx ,
     494                const vector< AD<Base> >&       apy
     495        );
    812496        // ------------------------------------------------------------------------
    813497        /*!
     
    820504                const vectorBool&                       r  ,
    821505                      vectorBool&                       s  ,
    822                 const vector<Base>&                     x  )
    823         {       return for_sparse_jac< vectorBool >(q, r, s, x);
    824         }
     506                const vector<Base>&                     x
     507        );
    825508        /*!
    826509        Link from user_atomic to forward sparse Jacobian bool
     
    832515                const vector<bool>&                     r  ,
    833516                      vector<bool>&                     s  ,
    834                 const vector<Base>&                     x  )
    835         {       return for_sparse_jac< vector<bool> >(q, r, s, x);
    836         }
     517                const vector<Base>&                     x
     518        );
    837519        /*!
    838520        Link from user_atomic to forward sparse Jacobian sets
     
    844526                const vector< std::set<size_t> >&       r  ,
    845527                      vector< std::set<size_t> >&       s  ,
    846                 const vector<Base>&                     x  )
    847         {       // during user sparsity calculations
    848                 size_t m = f_.Range();
    849                 size_t n = f_.Domain();
    850                 if( jac_sparse_bool_.size() != 0 )
    851                         jac_sparse_bool_.clear();
    852                 if( jac_sparse_set_.n_set() == 0 )
    853                         set_jac_sparse_set();
    854                 CPPAD_ASSERT_UNKNOWN( jac_sparse_bool_.size() == 0 );
    855                 CPPAD_ASSERT_UNKNOWN( jac_sparse_set_.n_set() == m );
    856                 CPPAD_ASSERT_UNKNOWN( jac_sparse_set_.end()   == n );
    857                 CPPAD_ASSERT_UNKNOWN( r.size() == n );
    858                 CPPAD_ASSERT_UNKNOWN( s.size() == m );
    859 
    860                 bool ok = true;
    861                 for(size_t i = 0; i < m; i++)
    862                         s[i].clear();
    863 
    864                 // sparsity for  s = jac_sparse_set_ * r
    865                 for(size_t i = 0; i < m; i++)
    866                 {       // compute row i of the return pattern
    867                         local::sparse_list::const_iterator set_itr(
    868                                 jac_sparse_set_, i
    869                         );
    870                         size_t j = *set_itr;
    871                         while(j < n )
    872                         {       std::set<size_t>::const_iterator itr_j;
    873                                 const std::set<size_t>& r_j( r[j] );
    874                                 for(itr_j = r_j.begin(); itr_j != r_j.end(); itr_j++)
    875                                 {       size_t k = *itr_j;
    876                                         CPPAD_ASSERT_UNKNOWN( k < q );
    877                                         s[i].insert(k);
    878                                 }
    879                                 j = *(++set_itr);
    880                         }
    881                 }
    882 
    883                 return ok;
    884         }
     528                const vector<Base>&                     x
     529        );
    885530        // ------------------------------------------------------------------------
    886531        /*!
     
    893538                const vectorBool&                       rt ,
    894539                      vectorBool&                       st ,
    895                 const vector<Base>&                     x  )
    896         {       return rev_sparse_jac< vectorBool >(q, rt, st, x);
    897         }
     540                const vector<Base>&                     x
     541        );
    898542        /*!
    899543        Link from user_atomic to reverse sparse Jacobian bool
     
    905549                const vector<bool>&                     rt ,
    906550                      vector<bool>&                     st ,
    907                 const vector<Base>&                     x  )
    908         {       return rev_sparse_jac< vector<bool> >(q, rt, st, x);
    909         }
     551                const vector<Base>&                     x
     552        );
    910553        /*!
    911554        Link from user_atomic to reverse Jacobian sets
     
    917560                const vector< std::set<size_t> >&       rt ,
    918561                      vector< std::set<size_t> >&       st ,
    919                 const vector<Base>&                     x  )
    920         {       // during user sparsity calculations
    921                 size_t m = f_.Range();
    922                 size_t n = f_.Domain();
    923                 if( jac_sparse_bool_.size() != 0 )
    924                         jac_sparse_bool_.clear();
    925                 if( jac_sparse_set_.n_set() == 0 )
    926                         set_jac_sparse_set();
    927                 CPPAD_ASSERT_UNKNOWN( jac_sparse_bool_.size() == 0 );
    928                 CPPAD_ASSERT_UNKNOWN( jac_sparse_set_.n_set() == m );
    929                 CPPAD_ASSERT_UNKNOWN( jac_sparse_set_.end()   == n );
    930                 CPPAD_ASSERT_UNKNOWN( rt.size() == m );
    931                 CPPAD_ASSERT_UNKNOWN( st.size() == n );
    932                 //
    933                 bool ok  = true;
    934                 //
    935                 for(size_t j = 0; j < n; j++)
    936                         st[j].clear();
    937                 //
    938                 // sparsity for  s = r * jac_sparse_set_
    939                 // s^T = jac_sparse_set_^T * r^T
    940                 for(size_t i = 0; i < m; i++)
    941                 {       // i is the row index in r^T
    942                         std::set<size_t>::const_iterator itr_i;
    943                         const std::set<size_t>& r_i( rt[i] );
    944                         for(itr_i = r_i.begin(); itr_i != r_i.end(); itr_i++)
    945                         {       // k is the column index in r^T
    946                                 size_t k = *itr_i;
    947                                 CPPAD_ASSERT_UNKNOWN( k < q );
    948                                 //
    949                                 // i is column index in jac_sparse_set^T
    950                                 local::sparse_list::const_iterator set_itr(
    951                                         jac_sparse_set_, i
    952                                 );
    953                                 size_t j = *set_itr;
    954                                 while( j < n )
    955                                 {       // j is row index in jac_sparse_set^T
    956                                         st[j].insert(k);
    957                                         j = *(++set_itr);
    958                                 }
    959                         }
    960                 }
    961 
    962                 return ok;
    963         }
     562                const vector<Base>&                     x
     563        );
    964564        // ------------------------------------------------------------------------
    965565        /*!
     
    976576                const vectorBool&                       u  ,
    977577                      vectorBool&                       v  ,
    978                 const vector<Base>&                     x  )
    979         {       return rev_sparse_hes< vectorBool >(vx, s, t, q, r, u, v, x);
    980         }
     578                const vector<Base>&                     x
     579        );
    981580        /*!
    982581        Link from user_atomic to reverse sparse Hessian bool
     
    992591                const vector<bool>&                     u  ,
    993592                      vector<bool>&                     v  ,
    994                 const vector<Base>&                     x  )
    995         {       return rev_sparse_hes< vector<bool> >(vx, s, t, q, r, u, v, x);
    996         }
     593                const vector<Base>&                     x
     594        );
    997595        /*!
    998596        Link from user_atomic to reverse sparse Hessian sets
     
    1008606                const vector< std::set<size_t> >&       u  ,
    1009607                      vector< std::set<size_t> >&       v  ,
    1010                 const vector<Base>&                     x  )
    1011         {       size_t n = f_.Domain();
    1012 # ifndef NDEBUG
    1013                 size_t m = f_.Range();
     608                const vector<Base>&                     x
     609        );
     610};
     611
     612} // END_CPPAD_NAMESPACE
     613
     614// functions implemented in cppad/core/checkpoint files
     615# include <cppad/core/checkpoint/ctor.hpp>
     616# include <cppad/core/checkpoint/reverse.hpp>
     617# include <cppad/core/checkpoint/forward.hpp>
     618# include <cppad/core/checkpoint/rev_sparse_hes.hpp>
     619# include <cppad/core/checkpoint/rev_sparse_jac.hpp>
     620# include <cppad/core/checkpoint/for_sparse_jac.hpp>
     621# include <cppad/core/checkpoint/set_hes_sparse_bool.hpp>
     622# include <cppad/core/checkpoint/set_hes_sparse_set.hpp>
     623# include <cppad/core/checkpoint/set_jac_sparse_bool.hpp>
     624# include <cppad/core/checkpoint/set_jac_sparse_set.hpp>
     625
     626# undef CPPAD_MULTI_THREAD_TMB
    1014627# endif
    1015                 CPPAD_ASSERT_UNKNOWN( vx.size() == n );
    1016                 CPPAD_ASSERT_UNKNOWN(  s.size() == m );
    1017                 CPPAD_ASSERT_UNKNOWN(  t.size() == n );
    1018                 CPPAD_ASSERT_UNKNOWN(  r.size() == n );
    1019                 CPPAD_ASSERT_UNKNOWN(  u.size() == m );
    1020                 CPPAD_ASSERT_UNKNOWN(  v.size() == n );
    1021                 //
    1022                 bool ok        = true;
    1023 
    1024                 // make sure hes_sparse_set_ has been set
    1025                 if( hes_sparse_bool_.size() != 0 )
    1026                         hes_sparse_bool_.clear();
    1027                 if( hes_sparse_set_.n_set() == 0 )
    1028                         set_hes_sparse_set();
    1029                 CPPAD_ASSERT_UNKNOWN( hes_sparse_bool_.size() == 0 );
    1030                 CPPAD_ASSERT_UNKNOWN( hes_sparse_set_.n_set() == n );
    1031                 CPPAD_ASSERT_UNKNOWN( hes_sparse_set_.end()   == n );
    1032 
    1033                 // compute sparsity pattern for T(x) = S(x) * f'(x)
    1034                 t = f_.RevSparseJac(1, s);
    1035 # ifndef NDEBUG
    1036                 for(size_t j = 0; j < n; j++)
    1037                         CPPAD_ASSERT_UNKNOWN( vx[j] || ! t[j] )
    1038 # endif
    1039 
    1040                 // V(x) = f'(x)^T * g''(y) * f'(x) * R  +  g'(y) * f''(x) * R
    1041                 // U(x) = g''(y) * f'(x) * R
    1042                 // S(x) = g'(y)
    1043 
    1044                 // compute sparsity pattern for A(x) = f'(x)^T * U(x)
    1045                 // 2DO: change a to use INTERNAL_SPARSE_SET
    1046                 bool transpose = true;
    1047                 vector< std::set<size_t> > a(n);
    1048                 a = f_.RevSparseJac(q, u, transpose);
    1049 
    1050                 // Need sparsity pattern for H(x) = (S(x) * f(x))''(x) * R,
    1051                 // but use less efficient sparsity for  f(x)''(x) * R so that
    1052                 // hes_sparse_set_ can be used every time this is needed.
    1053                 for(size_t i = 0; i < n; i++)
    1054                 {       v[i].clear();
    1055                         local::sparse_list::const_iterator set_itr(
    1056                                 hes_sparse_set_, i
    1057                         );
    1058                         size_t j = *set_itr;
    1059                         while( j < n )
    1060                         {       std::set<size_t>::const_iterator itr_j;
    1061                                 const std::set<size_t>& r_j( r[j] );
    1062                                 for(itr_j = r_j.begin(); itr_j != r_j.end(); itr_j++)
    1063                                 {       size_t k = *itr_j;
    1064                                         v[i].insert(k);
    1065                                 }
    1066                                 j = *(++set_itr);
    1067                         }
    1068                 }
    1069                 // compute sparsity pattern for V(x) = A(x) + H(x)
    1070                 std::set<size_t>::const_iterator itr;
    1071                 for(size_t i = 0; i < n; i++)
    1072                 {       for(itr = a[i].begin(); itr != a[i].end(); itr++)
    1073                         {       size_t j = *itr;
    1074                                 CPPAD_ASSERT_UNKNOWN( j < q );
    1075                                 v[i].insert(j);
    1076                         }
    1077                 }
    1078 
    1079                 return ok;
    1080         }
    1081 };
    1082 // ============================================================================
    1083 # else
    1084 // ============================================================================
    1085 # define THREAD omp_get_thread_num()
    1086 template <class Base>
    1087 class checkpoint : public atomic_base<Base> {
    1088 // ---------------------------------------------------------------------------
    1089 private:
    1090         /// same as option_enum in base class
    1091         typedef typename atomic_base<Base>::option_enum option_enum;
    1092         //
    1093         /// AD function corresponding to this checkpoint object
    1094         vector< ADFun<Base> > f_;
    1095         //
    1096         /// sparsity for entire Jacobian f(x)^{(1)} does not change so can cache it
    1097         vector<local::sparse_list> jac_sparse_set_;
    1098         vector<vectorBool>         jac_sparse_bool_;
    1099         //
    1100         /// sparsity for sum_i f_i(x)^{(2)} does not change so can cache it
    1101         vector<local::sparse_list> hes_sparse_set_;
    1102         vector<vectorBool>         hes_sparse_bool_;
    1103         // ------------------------------------------------------------------------
    1104         option_enum sparsity(void)
    1105         {       return static_cast< atomic_base<Base>* >(this)->sparsity(); }
    1106         // ------------------------------------------------------------------------
    1107         /// set jac_sparse_set_
    1108         void set_jac_sparse_set(void)
    1109         {       CPPAD_ASSERT_UNKNOWN( jac_sparse_set_[THREAD].n_set() == 0 );
    1110                 bool transpose  = false;
    1111                 bool dependency = true;
    1112                 size_t n = f_[THREAD].Domain();
    1113                 size_t m = f_[THREAD].Range();
    1114                 // Use the choice for forward / reverse that results in smaller
    1115                 // size for the sparsity pattern of all variables in the tape.
    1116                 if( n <= m )
    1117                 {       local::sparse_list identity;
    1118                         identity.resize(n, n);
    1119                         for(size_t j = 0; j < n; j++)
    1120                         {       // use add_element because only adding one element per set
    1121                                 identity.add_element(j, j);
    1122                         }
    1123                         f_[THREAD].ForSparseJacCheckpoint(
    1124                                 n, identity, transpose, dependency, jac_sparse_set_[THREAD]
    1125                         );
    1126                         f_[THREAD].size_forward_set(0);
    1127                 }
    1128                 else
    1129                 {       local::sparse_list identity;
    1130                         identity.resize(m, m);
    1131                         for(size_t i = 0; i < m; i++)
    1132                         {       // use add_element because only adding one element per set
    1133                                 identity.add_element(i, i);
    1134                         }
    1135                         f_[THREAD].RevSparseJacCheckpoint(
    1136                                 m, identity, transpose, dependency, jac_sparse_set_[THREAD]
    1137                         );
    1138                 }
    1139                 CPPAD_ASSERT_UNKNOWN( f_[THREAD].size_forward_set() == 0 );
    1140                 CPPAD_ASSERT_UNKNOWN( f_[THREAD].size_forward_bool() == 0 );
    1141         }
    1142         /// set jac_sparse_bool_
    1143         void set_jac_sparse_bool(void)
    1144         {       CPPAD_ASSERT_UNKNOWN( jac_sparse_bool_[THREAD].size() == 0 );
    1145                 bool transpose  = false;
    1146                 bool dependency = true;
    1147                 size_t n = f_[THREAD].Domain();
    1148                 size_t m = f_[THREAD].Range();
    1149                 // Use the choice for forward / reverse that results in smaller
    1150                 // size for the sparsity pattern of all variables in the tape.
    1151                 if( n <= m )
    1152                 {       vectorBool identity(n * n);
    1153                         for(size_t j = 0; j < n; j++)
    1154                         {       for(size_t i = 0; i < n; i++)
    1155                                         identity[ i * n + j ] = (i == j);
    1156                         }
    1157                         jac_sparse_bool_[THREAD] = f_[THREAD].ForSparseJac(
    1158                                 n, identity, transpose, dependency
    1159                         );
    1160                         f_[THREAD].size_forward_bool(0);
    1161                 }
    1162                 else
    1163                 {       vectorBool identity(m * m);
    1164                         for(size_t j = 0; j < m; j++)
    1165                         {       for(size_t i = 0; i < m; i++)
    1166                                         identity[ i * m + j ] = (i == j);
    1167                         }
    1168                         jac_sparse_bool_[THREAD] = f_[THREAD].RevSparseJac(
    1169                                 m, identity, transpose, dependency
    1170                         );
    1171                 }
    1172                 CPPAD_ASSERT_UNKNOWN( f_[THREAD].size_forward_bool() == 0 );
    1173                 CPPAD_ASSERT_UNKNOWN( f_[THREAD].size_forward_set() == 0 );
    1174         }
    1175         // ------------------------------------------------------------------------
    1176         /// set hes_sparse_set_
    1177         void set_hes_sparse_set(void)
    1178         {       CPPAD_ASSERT_UNKNOWN( hes_sparse_set_[THREAD].n_set() == 0 );
    1179                 size_t n = f_[THREAD].Domain();
    1180                 size_t m = f_[THREAD].Range();
    1181                 //
    1182                 // set version of sparsity for vector of all ones
    1183                 vector<bool> all_one(m);
    1184                 for(size_t i = 0; i < m; i++)
    1185                         all_one[i] = true;
    1186 
    1187                 // set version of sparsity for n by n idendity matrix
    1188                 local::sparse_list identity;
    1189                 identity.resize(n, n);
    1190                 for(size_t j = 0; j < n; j++)
    1191                 {       // use add_element because only adding one element per set
    1192                         identity.add_element(j, j);
    1193                 }
    1194 
    1195                 // compute sparsity pattern for H(x) = sum_i f_i(x)^{(2)}
    1196                 bool transpose  = false;
    1197                 bool dependency = false;
    1198                 f_[THREAD].ForSparseJacCheckpoint(
    1199                         n, identity, transpose, dependency, jac_sparse_set_[THREAD]
    1200                 );
    1201                 f_[THREAD].RevSparseHesCheckpoint(
    1202                         n, all_one, transpose, hes_sparse_set_[THREAD]
    1203                 );
    1204                 CPPAD_ASSERT_UNKNOWN( hes_sparse_set_[THREAD].n_set() == n );
    1205                 CPPAD_ASSERT_UNKNOWN( hes_sparse_set_[THREAD].end()   == n );
    1206                 //
    1207                 // drop the forward sparsity results from f_
    1208                 f_[THREAD].size_forward_set(0);
    1209         }
    1210         /// set hes_sparse_bool_
    1211         void set_hes_sparse_bool(void)
    1212         {       CPPAD_ASSERT_UNKNOWN( hes_sparse_bool_[THREAD].size() == 0 );
    1213                 size_t n = f_[THREAD].Domain();
    1214                 size_t m = f_[THREAD].Range();
    1215                 //
    1216                 // set version of sparsity for vector of all ones
    1217                 vectorBool all_one(m);
    1218                 for(size_t i = 0; i < m; i++)
    1219                         all_one[i] = true;
    1220 
    1221                 // set version of sparsity for n by n idendity matrix
    1222                 vectorBool identity(n * n);
    1223                 for(size_t j = 0; j < n; j++)
    1224                 {       for(size_t i = 0; i < n; i++)
    1225                                 identity[ i * n + j ] = (i == j);
    1226                 }
    1227 
    1228                 // compute sparsity pattern for H(x) = sum_i f_i(x)^{(2)}
    1229                 bool transpose  = false;
    1230                 bool dependency = false;
    1231                 f_[THREAD].ForSparseJac(n, identity, transpose, dependency);
    1232                 hes_sparse_bool_[THREAD] = f_[THREAD].RevSparseHes(n, all_one, transpose);
    1233                 CPPAD_ASSERT_UNKNOWN( hes_sparse_bool_[THREAD].size() == n * n );
    1234                 //
    1235                 // drop the forward sparsity results from f_
    1236                 f_[THREAD].size_forward_bool(0);
    1237                 CPPAD_ASSERT_UNKNOWN( f_[THREAD].size_forward_bool() == 0 );
    1238                 CPPAD_ASSERT_UNKNOWN( f_[THREAD].size_forward_set() == 0 );
    1239         }
    1240         // ------------------------------------------------------------------------
    1241         /*!
    1242         Link from user_atomic to forward sparse Jacobian pack and bool
    1243 
    1244         \copydetails atomic_base::for_sparse_jac
    1245         */
    1246         template <class sparsity_type>
    1247         bool for_sparse_jac(
    1248                 size_t                                  q  ,
    1249                 const sparsity_type&                    r  ,
    1250                       sparsity_type&                    s  ,
    1251                 const vector<Base>&                     x  )
    1252         {       // during user sparsity calculations
    1253                 size_t m = f_[THREAD].Range();
    1254                 size_t n = f_[THREAD].Domain();
    1255                 if( jac_sparse_bool_[THREAD].size() == 0 )
    1256                         set_jac_sparse_bool();
    1257                 if( jac_sparse_set_[THREAD].n_set() != 0 )
    1258                         jac_sparse_set_[THREAD].resize(0, 0);
    1259                 CPPAD_ASSERT_UNKNOWN( jac_sparse_bool_[THREAD].size() == m * n );
    1260                 CPPAD_ASSERT_UNKNOWN( jac_sparse_set_[THREAD].n_set() == 0 );
    1261                 CPPAD_ASSERT_UNKNOWN( r.size() == n * q );
    1262                 CPPAD_ASSERT_UNKNOWN( s.size() == m * q );
    1263                 //
    1264                 bool ok = true;
    1265                 for(size_t i = 0; i < m; i++)
    1266                 {       for(size_t k = 0; k < q; k++)
    1267                                 s[i * q + k] = false;
    1268                 }
    1269                 // sparsity for  s = jac_sparse_bool_ * r
    1270                 for(size_t i = 0; i < m; i++)
    1271                 {       for(size_t k = 0; k < q; k++)
    1272                         {       // initialize sparsity for S(i,k)
    1273                                 bool s_ik = false;
    1274                                 // S(i,k) = sum_j J(i,j) * R(j,k)
    1275                                 for(size_t j = 0; j < n; j++)
    1276                                 {       bool J_ij = jac_sparse_bool_[THREAD][ i * n + j];
    1277                                         bool R_jk = r[j * q + k ];
    1278                                         s_ik |= ( J_ij & R_jk );
    1279                                 }
    1280                                 s[i * q + k] = s_ik;
    1281                         }
    1282                 }
    1283                 return ok;
    1284         }
    1285         // ------------------------------------------------------------------------
    1286         /*!
    1287         Link from user_atomic to reverse sparse Jacobian pack and bool
    1288 
    1289         \copydetails atomic_base::rev_sparse_jac
    1290         */
    1291         template <class sparsity_type>
    1292         bool rev_sparse_jac(
    1293                 size_t                                  q  ,
    1294                 const sparsity_type&                    rt ,
    1295                       sparsity_type&                    st ,
    1296                 const vector<Base>&                     x  )
    1297         {       // during user sparsity calculations
    1298                 size_t m = f_[THREAD].Range();
    1299                 size_t n = f_[THREAD].Domain();
    1300                 if( jac_sparse_bool_[THREAD].size() == 0 )
    1301                         set_jac_sparse_bool();
    1302                 if( jac_sparse_set_[THREAD].n_set() != 0 )
    1303                         jac_sparse_set_[THREAD].resize(0, 0);
    1304                 CPPAD_ASSERT_UNKNOWN( jac_sparse_bool_[THREAD].size() == m * n );
    1305                 CPPAD_ASSERT_UNKNOWN( jac_sparse_set_[THREAD].n_set() == 0 );
    1306                 CPPAD_ASSERT_UNKNOWN( rt.size() == m * q );
    1307                 CPPAD_ASSERT_UNKNOWN( st.size() == n * q );
    1308                 bool ok  = true;
    1309                 //
    1310                 // S = R * J where J is jacobian
    1311                 for(size_t i = 0; i < q; i++)
    1312                 {       for(size_t j = 0; j < n; j++)
    1313                         {       // initialize sparsity for S(i,j)
    1314                                 bool s_ij = false;
    1315                                 // S(i,j) = sum_k R(i,k) * J(k,j)
    1316                                 for(size_t k = 0; k < m; k++)
    1317                                 {       // sparsity for R(i, k)
    1318                                         bool R_ik = rt[ k * q + i ];
    1319                                         bool J_kj = jac_sparse_bool_[THREAD][ k * n + j ];
    1320                                         s_ij     |= (R_ik & J_kj);
    1321                                 }
    1322                                 // set sparsity for S^T
    1323                                 st[ j * q + i ] = s_ij;
    1324                         }
    1325                 }
    1326                 return ok;
    1327         }
    1328         /*!
    1329         Link from user_atomic to reverse sparse Hessian  bools
    1330 
    1331         \copydetails atomic_base::rev_sparse_hes
    1332         */
    1333         template <class sparsity_type>
    1334         bool rev_sparse_hes(
    1335                 const vector<bool>&                     vx ,
    1336                 const vector<bool>&                     s  ,
    1337                       vector<bool>&                     t  ,
    1338                 size_t                                  q  ,
    1339                 const sparsity_type&                    r  ,
    1340                 const sparsity_type&                    u  ,
    1341                       sparsity_type&                    v  ,
    1342                 const vector<Base>&                     x  )
    1343         {       size_t n = f_[THREAD].Domain();
    1344 # ifndef NDEBUG
    1345                 size_t m = f_[THREAD].Range();
    1346 # endif
    1347                 CPPAD_ASSERT_UNKNOWN( vx.size() == n );
    1348                 CPPAD_ASSERT_UNKNOWN(  s.size() == m );
    1349                 CPPAD_ASSERT_UNKNOWN(  t.size() == n );
    1350                 CPPAD_ASSERT_UNKNOWN(  r.size() == n * q );
    1351                 CPPAD_ASSERT_UNKNOWN(  u.size() == m * q );
    1352                 CPPAD_ASSERT_UNKNOWN(  v.size() == n * q );
    1353                 //
    1354                 bool ok        = true;
    1355 
    1356                 // make sure hes_sparse_bool_ has been set
    1357                 if( hes_sparse_bool_[THREAD].size() == 0 )
    1358                         set_hes_sparse_bool();
    1359                 if( hes_sparse_set_[THREAD].n_set() != 0 )
    1360                         hes_sparse_set_[THREAD].resize(0, 0);
    1361                 CPPAD_ASSERT_UNKNOWN( hes_sparse_bool_[THREAD].size() == n * n );
    1362                 CPPAD_ASSERT_UNKNOWN( hes_sparse_set_[THREAD].n_set() == 0 );
    1363 
    1364 
    1365                 // compute sparsity pattern for T(x) = S(x) * f'(x)
    1366                 t = f_[THREAD].RevSparseJac(1, s);
    1367 # ifndef NDEBUG
    1368                 for(size_t j = 0; j < n; j++)
    1369                         CPPAD_ASSERT_UNKNOWN( vx[j] || ! t[j] )
    1370 # endif
    1371 
    1372                 // V(x) = f'(x)^T * g''(y) * f'(x) * R  +  g'(y) * f''(x) * R
    1373                 // U(x) = g''(y) * f'(x) * R
    1374                 // S(x) = g'(y)
    1375 
    1376                 // compute sparsity pattern for A(x) = f'(x)^T * U(x)
    1377                 bool transpose = true;
    1378                 sparsity_type a(n * q);
    1379                 a = f_[THREAD].RevSparseJac(q, u, transpose);
    1380 
    1381                 // Need sparsity pattern for H(x) = (S(x) * f(x))''(x) * R,
    1382                 // but use less efficient sparsity for  f(x)''(x) * R so that
    1383                 // hes_sparse_set_ can be used every time this is needed.
    1384                 for(size_t i = 0; i < n; i++)
    1385                 {       for(size_t k = 0; k < q; k++)
    1386                         {       // initialize sparsity pattern for H(i,k)
    1387                                 bool h_ik = false;
    1388                                 // H(i,k) = sum_j f''(i,j) * R(j,k)
    1389                                 for(size_t j = 0; j < n; j++)
    1390                                 {       bool f_ij = hes_sparse_bool_[THREAD][i * n + j];
    1391                                         bool r_jk = r[j * q + k];
    1392                                         h_ik     |= ( f_ij & r_jk );
    1393                                 }
    1394                                 // sparsity for H(i,k)
    1395                                 v[i * q + k] = h_ik;
    1396                         }
    1397                 }
    1398 
    1399                 // compute sparsity pattern for V(x) = A(x) + H(x)
    1400                 for(size_t i = 0; i < n; i++)
    1401                 {       for(size_t k = 0; k < q; k++)
    1402                                 // v[ i * q + k ] |= a[ i * q + k];
    1403                                 v[ i * q + k ] = bool(v[ i * q + k]) | bool(a[ i * q + k]);
    1404                 }
    1405                 return ok;
    1406         }
    1407 public:
    1408         // ------------------------------------------------------------------------
    1409         /*!
    1410         Constructor of a checkpoint object
    1411 
    1412         \param name [in]
    1413         is the user's name for the AD version of this atomic operation.
    1414 
    1415         \param algo [in/out]
    1416         user routine that compute AD function values
    1417         (not const because state may change during evaluation).
    1418 
    1419         \param ax [in]
    1420         argument value where algo operation sequence is taped.
    1421 
    1422         \param ay [out]
    1423         function value at specified argument value.
    1424 
    1425         \param sparsity [in]
    1426         what type of sparsity patterns are computed by this function,
    1427         pack_sparsity_enum bool_sparsity_enum, or set_sparsity_enum.
    1428         The default value is unspecified.
    1429 
    1430         \param optimize [in]
    1431         should the operation sequence corresponding to the algo be optimized.
    1432         The default value is true, but it is
    1433         sometimes useful to use false for debugging purposes.
    1434         */
    1435         template <class Algo, class ADVector>
    1436         checkpoint(
    1437                 const char*                    name            ,
    1438                 Algo&                          algo            ,
    1439                 const ADVector&                ax              ,
    1440                 ADVector&                      ay              ,
    1441                 option_enum                    sparsity =
    1442                                 atomic_base<Base>::pack_sparsity_enum  ,
    1443                 bool                           optimize = true
    1444         ) :
    1445         atomic_base<Base>(name  , sparsity)        ,
    1446         f_( size_t( omp_get_max_threads() ) )                ,
    1447         jac_sparse_set_( size_t( omp_get_max_threads() ) )   ,
    1448         jac_sparse_bool_( size_t( omp_get_max_threads() ) )  ,
    1449         hes_sparse_set_( size_t( omp_get_max_threads() ) )   ,
    1450         hes_sparse_bool_( size_t( omp_get_max_threads() ) )
    1451         {
    1452                 CheckSimpleVector< CppAD::AD<Base> , ADVector>();
    1453 
    1454                 // make a copy of ax because Independent modifies AD information
    1455                 ADVector x_tmp(ax);
    1456                 // delcare x_tmp as the independent variables
    1457                 Independent(x_tmp);
    1458                 // record mapping from x_tmp to ay
    1459                 algo(x_tmp, ay);
    1460                 // create function f_ : x -> y
    1461                 f_[0].Dependent(ay);
    1462                 if( optimize )
    1463                 {       // suppress checking for nan in f_ results
    1464                         // (see optimize documentation for atomic functions)
    1465                         f_[0].check_for_nan(false);
    1466                         //
    1467                         // now optimize
    1468                         f_[0].optimize();
    1469                 }
    1470                 // now disable checking of comparison operations
    1471                 // 2DO: add a debugging mode that checks for changes and aborts
    1472                 f_[0].compare_change_count(0);
    1473                 //
    1474                 // copy for use during multi-threading
    1475                 for(int i = 1; i < omp_get_max_threads(); ++i)
    1476                         f_[i] = f_[0];
    1477         }
    1478         // ------------------------------------------------------------------------
    1479         /*!
    1480         Implement the user call to <tt>atom_fun.size_var()</tt>.
    1481         */
    1482         size_t size_var(void)
    1483         {       return f_[THREAD].size_var(); }
    1484         // ------------------------------------------------------------------------
    1485         /*!
    1486         Implement the user call to <tt>atom_fun(ax, ay)</tt>.
    1487 
    1488         \tparam ADVector
    1489         A simple vector class with elements of type <code>AD<Base></code>.
    1490 
    1491         \param id
    1492         optional parameter which must be zero if present.
    1493 
    1494         \param ax
    1495         is the argument vector for this call,
    1496         <tt>ax.size()</tt> determines the number of arguments.
    1497 
    1498         \param ay
    1499         is the result vector for this call,
    1500         <tt>ay.size()</tt> determines the number of results.
    1501         */
    1502         template <class ADVector>
    1503         void operator()(const ADVector& ax, ADVector& ay, size_t id = 0)
    1504         {       CPPAD_ASSERT_KNOWN(
    1505                         id == 0,
    1506                         "checkpoint: id is non-zero in atom_fun(ax, ay, id)"
    1507                 );
    1508                 this->atomic_base<Base>::operator()(ax, ay, id);
    1509         }
    1510         // ------------------------------------------------------------------------
    1511         /*!
    1512         Link from user_atomic to forward mode
    1513 
    1514         \copydetails atomic_base::forward
    1515         */
    1516         virtual bool forward(
    1517                 size_t                    p ,
    1518                 size_t                    q ,
    1519                 const vector<bool>&      vx ,
    1520                       vector<bool>&      vy ,
    1521                 const vector<Base>&      tx ,
    1522                       vector<Base>&      ty )
    1523         {       size_t n = f_[THREAD].Domain();
    1524                 size_t m = f_[THREAD].Range();
    1525                 //
    1526                 CPPAD_ASSERT_UNKNOWN( f_[THREAD].size_var() > 0 );
    1527                 CPPAD_ASSERT_UNKNOWN( tx.size() % (q+1) == 0 );
    1528                 CPPAD_ASSERT_UNKNOWN( ty.size() % (q+1) == 0 );
    1529                 CPPAD_ASSERT_UNKNOWN( n == tx.size() / (q+1) );
    1530                 CPPAD_ASSERT_UNKNOWN( m == ty.size() / (q+1) );
    1531                 bool ok  = true;
    1532                 //
    1533                 if( vx.size() == 0 )
    1534                 {       // during user forward mode
    1535                         if( jac_sparse_set_[THREAD].n_set() != 0 )
    1536                                 jac_sparse_set_[THREAD].resize(0,0);
    1537                         if( jac_sparse_bool_[THREAD].size() != 0 )
    1538                                 jac_sparse_bool_[THREAD].clear();
    1539                         //
    1540                         if( hes_sparse_set_[THREAD].n_set() != 0 )
    1541                                 hes_sparse_set_[THREAD].resize(0,0);
    1542                         if( hes_sparse_bool_[THREAD].size() != 0 )
    1543                                 hes_sparse_bool_[THREAD].clear();
    1544                 }
    1545                 if( vx.size() > 0 )
    1546                 {       // need Jacobian sparsity pattern to determine variable relation
    1547                         // during user recording using checkpoint functions
    1548                         if( sparsity() == atomic_base<Base>::set_sparsity_enum )
    1549                         {       if( jac_sparse_set_[THREAD].n_set() == 0 )
    1550                                         set_jac_sparse_set();
    1551                                 CPPAD_ASSERT_UNKNOWN( jac_sparse_set_[THREAD].n_set() == m );
    1552                                 CPPAD_ASSERT_UNKNOWN( jac_sparse_set_[THREAD].end()   == n );
    1553                                 //
    1554                                 for(size_t i = 0; i < m; i++)
    1555                                 {       vy[i] = false;
    1556                                         local::sparse_list::const_iterator set_itr(
    1557                                                 jac_sparse_set_[THREAD], i
    1558                                         );
    1559                                         size_t j = *set_itr;
    1560                                         while(j < n )
    1561                                         {       // y[i] depends on the value of x[j]
    1562                                                 // cast avoid Microsoft warning (should not be needed)
    1563                                                 vy[i] |= static_cast<bool>( vx[j] );
    1564                                                 j = *(++set_itr);
    1565                                         }
    1566                                 }
    1567                         }
    1568                         else
    1569                         {       if( jac_sparse_set_[THREAD].n_set() != 0 )
    1570                                         jac_sparse_set_[THREAD].resize(0, 0);
    1571                                 if( jac_sparse_bool_[THREAD].size() == 0 )
    1572                                         set_jac_sparse_bool();
    1573                                 CPPAD_ASSERT_UNKNOWN( jac_sparse_set_[THREAD].n_set() == 0 );
    1574                                 CPPAD_ASSERT_UNKNOWN( jac_sparse_bool_[THREAD].size() == m * n );
    1575                                 //
    1576                                 for(size_t i = 0; i < m; i++)
    1577                                 {       vy[i] = false;
    1578                                         for(size_t j = 0; j < n; j++)
    1579                                         {       if( jac_sparse_bool_[THREAD][ i * n + j ] )
    1580                                                 {       // y[i] depends on the value of x[j]
    1581                                                         // cast avoid Microsoft warning
    1582                                                         vy[i] |= static_cast<bool>( vx[j] );
    1583                                                 }
    1584                                         }
    1585                                 }
    1586                         }
    1587                 }
    1588                 // compute forward results for orders zero through q
    1589                 ty = f_[THREAD].Forward(q, tx);
    1590 
    1591                 // no longer need the Taylor coefficients in f_
    1592                 // (have to reconstruct them every time)
    1593                 // Hold onto sparsity pattern because it is always good.
    1594                 size_t c = 0;
    1595                 size_t r = 0;
    1596                 f_[THREAD].capacity_order(c, r);
    1597                 return ok;
    1598         }
    1599         // ------------------------------------------------------------------------
    1600         /*!
    1601         Link from user_atomic to reverse mode
    1602 
    1603         \copydetails atomic_base::reverse
    1604         */
    1605         virtual bool reverse(
    1606                 size_t                    q  ,
    1607                 const vector<Base>&       tx ,
    1608                 const vector<Base>&       ty ,
    1609                       vector<Base>&       px ,
    1610                 const vector<Base>&       py )
    1611         {
    1612 # ifndef NDEBUG
    1613                 size_t n = f_[THREAD].Domain();
    1614                 size_t m = f_[THREAD].Range();
    1615 # endif
    1616                 CPPAD_ASSERT_UNKNOWN( n == tx.size() / (q+1) );
    1617                 CPPAD_ASSERT_UNKNOWN( m == ty.size() / (q+1) );
    1618                 CPPAD_ASSERT_UNKNOWN( f_[THREAD].size_var() > 0 );
    1619                 CPPAD_ASSERT_UNKNOWN( tx.size() % (q+1) == 0 );
    1620                 CPPAD_ASSERT_UNKNOWN( ty.size() % (q+1) == 0 );
    1621                 bool ok  = true;
    1622 
    1623                 // put proper forward mode coefficients in f_
    1624 # ifdef NDEBUG
    1625                 // compute forward results for orders zero through q
    1626                 f_[THREAD].Forward(q, tx);
    1627 # else
    1628                 CPPAD_ASSERT_UNKNOWN( px.size() == n * (q+1) );
    1629                 CPPAD_ASSERT_UNKNOWN( py.size() == m * (q+1) );
    1630                 size_t i, j, k;
    1631                 //
    1632                 // compute forward results for orders zero through q
    1633                 vector<Base> check_ty = f_[THREAD].Forward(q, tx);
    1634                 for(i = 0; i < m; i++)
    1635                 {       for(k = 0; k <= q; k++)
    1636                         {       j = i * (q+1) + k;
    1637                                 CPPAD_ASSERT_UNKNOWN( check_ty[j] == ty[j] );
    1638                         }
    1639                 }
    1640 # endif
    1641                 // now can run reverse mode
    1642                 px = f_[THREAD].Reverse(q+1, py);
    1643 
    1644                 // no longer need the Taylor coefficients in f_
    1645                 // (have to reconstruct them every time)
    1646                 size_t c = 0;
    1647                 size_t r = 0;
    1648                 f_[THREAD].capacity_order(c, r);
    1649                 return ok;
    1650         }
    1651         // ------------------------------------------------------------------------
    1652         /*!
    1653         Link from user_atomic to forward sparse Jacobian pack
    1654 
    1655         \copydetails atomic_base::for_sparse_jac
    1656         */
    1657         virtual bool for_sparse_jac(
    1658                 size_t                                  q  ,
    1659                 const vectorBool&                       r  ,
    1660                       vectorBool&                       s  ,
    1661                 const vector<Base>&                     x  )
    1662         {       return for_sparse_jac< vectorBool >(q, r, s, x);
    1663         }
    1664         /*!
    1665         Link from user_atomic to forward sparse Jacobian bool
    1666 
    1667         \copydetails atomic_base::for_sparse_jac
    1668         */
    1669         virtual bool for_sparse_jac(
    1670                 size_t                                  q  ,
    1671                 const vector<bool>&                     r  ,
    1672                       vector<bool>&                     s  ,
    1673                 const vector<Base>&                     x  )
    1674         {       return for_sparse_jac< vector<bool> >(q, r, s, x);
    1675         }
    1676         /*!
    1677         Link from user_atomic to forward sparse Jacobian sets
    1678 
    1679         \copydetails atomic_base::for_sparse_jac
    1680         */
    1681         virtual bool for_sparse_jac(
    1682                 size_t                                  q  ,
    1683                 const vector< std::set<size_t> >&       r  ,
    1684                       vector< std::set<size_t> >&       s  ,
    1685                 const vector<Base>&                     x  )
    1686         {       // during user sparsity calculations
    1687                 size_t m = f_[THREAD].Range();
    1688                 size_t n = f_[THREAD].Domain();
    1689                 if( jac_sparse_bool_[THREAD].size() != 0 )
    1690                         jac_sparse_bool_[THREAD].clear();
    1691                 if( jac_sparse_set_[THREAD].n_set() == 0 )
    1692                         set_jac_sparse_set();
    1693                 CPPAD_ASSERT_UNKNOWN( jac_sparse_bool_[THREAD].size() == 0 );
    1694                 CPPAD_ASSERT_UNKNOWN( jac_sparse_set_[THREAD].n_set() == m );
    1695                 CPPAD_ASSERT_UNKNOWN( jac_sparse_set_[THREAD].end()   == n );
    1696                 CPPAD_ASSERT_UNKNOWN( r.size() == n );
    1697                 CPPAD_ASSERT_UNKNOWN( s.size() == m );
    1698 
    1699                 bool ok = true;
    1700                 for(size_t i = 0; i < m; i++)
    1701                         s[i].clear();
    1702 
    1703                 // sparsity for  s = jac_sparse_set_ * r
    1704                 for(size_t i = 0; i < m; i++)
    1705                 {       // compute row i of the return pattern
    1706                         local::sparse_list::const_iterator set_itr(
    1707                                 jac_sparse_set_[THREAD], i
    1708                         );
    1709                         size_t j = *set_itr;
    1710                         while(j < n )
    1711                         {       std::set<size_t>::const_iterator itr_j;
    1712                                 const std::set<size_t>& r_j( r[j] );
    1713                                 for(itr_j = r_j.begin(); itr_j != r_j.end(); itr_j++)
    1714                                 {       size_t k = *itr_j;
    1715                                         CPPAD_ASSERT_UNKNOWN( k < q );
    1716                                         s[i].insert(k);
    1717                                 }
    1718                                 j = *(++set_itr);
    1719                         }
    1720                 }
    1721 
    1722                 return ok;
    1723         }
    1724         // ------------------------------------------------------------------------
    1725         /*!
    1726         Link from user_atomic to reverse sparse Jacobian pack
    1727 
    1728         \copydetails atomic_base::rev_sparse_jac
    1729         */
    1730         virtual bool rev_sparse_jac(
    1731                 size_t                                  q  ,
    1732                 const vectorBool&                       rt ,
    1733                       vectorBool&                       st ,
    1734                 const vector<Base>&                     x  )
    1735         {       return rev_sparse_jac< vectorBool >(q, rt, st, x);
    1736         }
    1737         /*!
    1738         Link from user_atomic to reverse sparse Jacobian bool
    1739 
    1740         \copydetails atomic_base::rev_sparse_jac
    1741         */
    1742         virtual bool rev_sparse_jac(
    1743                 size_t                                  q  ,
    1744                 const vector<bool>&                     rt ,
    1745                       vector<bool>&                     st ,
    1746                 const vector<Base>&                     x  )
    1747         {       return rev_sparse_jac< vector<bool> >(q, rt, st, x);
    1748         }
    1749         /*!
    1750         Link from user_atomic to reverse Jacobian sets
    1751 
    1752         \copydetails atomic_base::rev_sparse_jac
    1753         */
    1754         virtual bool rev_sparse_jac(
    1755                 size_t                                  q  ,
    1756                 const vector< std::set<size_t> >&       rt ,
    1757                       vector< std::set<size_t> >&       st ,
    1758                 const vector<Base>&                     x  )
    1759         {       // during user sparsity calculations
    1760                 size_t m = f_[THREAD].Range();
    1761                 size_t n = f_[THREAD].Domain();
    1762                 if( jac_sparse_bool_[THREAD].size() != 0 )
    1763                         jac_sparse_bool_[THREAD].clear();
    1764                 if( jac_sparse_set_[THREAD].n_set() == 0 )
    1765                         set_jac_sparse_set();
    1766                 CPPAD_ASSERT_UNKNOWN( jac_sparse_bool_[THREAD].size() == 0 );
    1767                 CPPAD_ASSERT_UNKNOWN( jac_sparse_set_[THREAD].n_set() == m );
    1768                 CPPAD_ASSERT_UNKNOWN( jac_sparse_set_[THREAD].end()   == n );
    1769                 CPPAD_ASSERT_UNKNOWN( rt.size() == m );
    1770                 CPPAD_ASSERT_UNKNOWN( st.size() == n );
    1771                 //
    1772                 bool ok  = true;
    1773                 //
    1774                 for(size_t j = 0; j < n; j++)
    1775                         st[j].clear();
    1776                 //
    1777                 // sparsity for  s = r * jac_sparse_set_
    1778                 // s^T = jac_sparse_set_^T * r^T
    1779                 for(size_t i = 0; i < m; i++)
    1780                 {       // i is the row index in r^T
    1781                         std::set<size_t>::const_iterator itr_i;
    1782                         const std::set<size_t>& r_i( rt[i] );
    1783                         for(itr_i = r_i.begin(); itr_i != r_i.end(); itr_i++)
    1784                         {       // k is the column index in r^T
    1785                                 size_t k = *itr_i;
    1786                                 CPPAD_ASSERT_UNKNOWN( k < q );
    1787                                 //
    1788                                 // i is column index in jac_sparse_set^T
    1789                                 local::sparse_list::const_iterator set_itr(
    1790                                         jac_sparse_set_[THREAD], i
    1791                                 );
    1792                                 size_t j = *set_itr;
    1793                                 while( j < n )
    1794                                 {       // j is row index in jac_sparse_set^T
    1795                                         st[j].insert(k);
    1796                                         j = *(++set_itr);
    1797                                 }
    1798                         }
    1799                 }
    1800 
    1801                 return ok;
    1802         }
    1803         // ------------------------------------------------------------------------
    1804         /*!
    1805         Link from user_atomic to reverse sparse Hessian pack
    1806 
    1807         \copydetails atomic_base::rev_sparse_hes
    1808         */
    1809         virtual bool rev_sparse_hes(
    1810                 const vector<bool>&                     vx ,
    1811                 const vector<bool>&                     s  ,
    1812                       vector<bool>&                     t  ,
    1813                 size_t                                  q  ,
    1814                 const vectorBool&                       r  ,
    1815                 const vectorBool&                       u  ,
    1816                       vectorBool&                       v  ,
    1817                 const vector<Base>&                     x  )
    1818         {       return rev_sparse_hes< vectorBool >(vx, s, t, q, r, u, v, x);
    1819         }
    1820         /*!
    1821         Link from user_atomic to reverse sparse Hessian bool
    1822 
    1823         \copydetails atomic_base::rev_sparse_hes
    1824         */
    1825         virtual bool rev_sparse_hes(
    1826                 const vector<bool>&                     vx ,
    1827                 const vector<bool>&                     s  ,
    1828                       vector<bool>&                     t  ,
    1829                 size_t                                  q  ,
    1830                 const vector<bool>&                     r  ,
    1831                 const vector<bool>&                     u  ,
    1832                       vector<bool>&                     v  ,
    1833                 const vector<Base>&                     x  )
    1834         {       return rev_sparse_hes< vector<bool> >(vx, s, t, q, r, u, v, x);
    1835         }
    1836         /*!
    1837         Link from user_atomic to reverse sparse Hessian sets
    1838 
    1839         \copydetails atomic_base::rev_sparse_hes
    1840         */
    1841         virtual bool rev_sparse_hes(
    1842                 const vector<bool>&                     vx ,
    1843                 const vector<bool>&                     s  ,
    1844                       vector<bool>&                     t  ,
    1845                 size_t                                  q  ,
    1846                 const vector< std::set<size_t> >&       r  ,
    1847                 const vector< std::set<size_t> >&       u  ,
    1848                       vector< std::set<size_t> >&       v  ,
    1849                 const vector<Base>&                     x  )
    1850         {       size_t n = f_[THREAD].Domain();
    1851 # ifndef NDEBUG
    1852                 size_t m = f_[THREAD].Range();
    1853 # endif
    1854                 CPPAD_ASSERT_UNKNOWN( vx.size() == n );
    1855                 CPPAD_ASSERT_UNKNOWN(  s.size() == m );
    1856                 CPPAD_ASSERT_UNKNOWN(  t.size() == n );
    1857                 CPPAD_ASSERT_UNKNOWN(  r.size() == n );
    1858                 CPPAD_ASSERT_UNKNOWN(  u.size() == m );
    1859                 CPPAD_ASSERT_UNKNOWN(  v.size() == n );
    1860                 //
    1861                 bool ok        = true;
    1862 
    1863                 // make sure hes_sparse_set_ has been set
    1864                 if( hes_sparse_bool_[THREAD].size() != 0 )
    1865                         hes_sparse_bool_[THREAD].clear();
    1866                 if( hes_sparse_set_[THREAD].n_set() == 0 )
    1867                         set_hes_sparse_set();
    1868                 CPPAD_ASSERT_UNKNOWN( hes_sparse_bool_[THREAD].size() == 0 );
    1869                 CPPAD_ASSERT_UNKNOWN( hes_sparse_set_[THREAD].n_set() == n );
    1870                 CPPAD_ASSERT_UNKNOWN( hes_sparse_set_[THREAD].end()   == n );
    1871 
    1872                 // compute sparsity pattern for T(x) = S(x) * f'(x)
    1873                 t = f_[THREAD].RevSparseJac(1, s);
    1874 # ifndef NDEBUG
    1875                 for(size_t j = 0; j < n; j++)
    1876                         CPPAD_ASSERT_UNKNOWN( vx[j] || ! t[j] )
    1877 # endif
    1878 
    1879                 // V(x) = f'(x)^T * g''(y) * f'(x) * R  +  g'(y) * f''(x) * R
    1880                 // U(x) = g''(y) * f'(x) * R
    1881                 // S(x) = g'(y)
    1882 
    1883                 // compute sparsity pattern for A(x) = f'(x)^T * U(x)
    1884                 // 2DO: change a to use INTERNAL_SPARSE_SET
    1885                 bool transpose = true;
    1886                 vector< std::set<size_t> > a(n);
    1887                 a = f_[THREAD].RevSparseJac(q, u, transpose);
    1888 
    1889                 // Need sparsity pattern for H(x) = (S(x) * f(x))''(x) * R,
    1890                 // but use less efficient sparsity for  f(x)''(x) * R so that
    1891                 // hes_sparse_set_ can be used every time this is needed.
    1892                 for(size_t i = 0; i < n; i++)
    1893                 {       v[i].clear();
    1894                         local::sparse_list::const_iterator set_itr(
    1895                                 hes_sparse_set_[THREAD], i
    1896                         );
    1897                         size_t j = *set_itr;
    1898                         while( j < n )
    1899                         {       std::set<size_t>::const_iterator itr_j;
    1900                                 const std::set<size_t>& r_j( r[j] );
    1901                                 for(itr_j = r_j.begin(); itr_j != r_j.end(); itr_j++)
    1902                                 {       size_t k = *itr_j;
    1903                                         v[i].insert(k);
    1904                                 }
    1905                                 j = *(++set_itr);
    1906                         }
    1907                 }
    1908                 // compute sparsity pattern for V(x) = A(x) + H(x)
    1909                 std::set<size_t>::const_iterator itr;
    1910                 for(size_t i = 0; i < n; i++)
    1911                 {       for(itr = a[i].begin(); itr != a[i].end(); itr++)
    1