Changeset 2859


Ignore:
Timestamp:
May 28, 2013 2:03:21 AM (7 years ago)
Author:
bradbell
Message:

merge in changes from branches/atomic; see bin/svn_merge.sh

Location:
trunk
Files:
8 deleted
96 edited
11 copied

Legend:

Unmodified
Added
Removed
  • trunk

  • trunk/AUTHORS

    r2792 r2859  
    22             ===========================================
    33
    4 To date, 2013-04-28, Bradley M. Bell is the sole author of CppAD.
     4To date, 2013-05-27, Bradley M. Bell is the sole author of CppAD.
    55While Bradley M. Bell worked for the University of Washington during
    66the development of CppAD, the following are also true:
  • trunk/CMakeLists.txt

    r2792 r2859  
    1717
    1818# cppad_version is used by set_version.sh to get the version number.
    19 SET(cppad_version      "20130428" )
     19SET(cppad_version      "20130527" )
    2020SET(cppad_url          "http://www.coin-or.org/CppAD" )
    2121SET(cppad_description  "Differentiation of C++ Algorithms" )
  • trunk/bin/check_include_def.sh

    r2506 r2859  
    22# $Id$
    33# -----------------------------------------------------------------------------
    4 # CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
     4# CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-13 Bradley M. Bell
    55#
    66# CppAD is distributed under multiple licenses. This distribution is under
     
    2727        cppad/speed/*.hpp \
    2828        example/*.hpp \
     29        example/atomic/*.hpp \
    2930        multi_thread/*.hpp \
    3031        | sed \
     
    4142        cppad/speed/*.hpp \
    4243        example/*.hpp \
     44        example/atomic/*.hpp \
    4345        multi_thread/*.hpp \
    4446        | sed -e 's|.*/||' \
  • trunk/bin/check_include_file.sh

    r2751 r2859  
    2929do
    3030        dir_list=`find . -name "*$ext" | sed -e '/junk\.[^.]*$/d' \
    31                 -e 's|^\./||' -e '/^work/d' -e 's|/[^/]*$||' | sort -u` 
     31                -e 's|^\./||' -e '/^build/d' -e 's|/[^/]*$||' | sort -u` 
    3232        for dir in $dir_list
    3333        do
  • trunk/bin/doxyfile.sh

    r2621 r2859  
    22# $Id$
    33# -----------------------------------------------------------------------------
    4 # CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
     4# CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-13 Bradley M. Bell
    55#
    66# CppAD is distributed under multiple licenses. This distribution is under
     
    2626error_file="$2"
    2727output_directory="$3"
     28# 2DO: change EXTRACT_ALL to NO so get warnings for undocumented functions.
    2829echo "create bin/doxyfile.$$"
    2930cat << EOF > bin/doxyfile.$$
  • trunk/bin/search.sh

    r2770 r2859  
    3838        cppad_ipopt/test
    3939        example
     40        example/atomic
     41        example/ipopt_solve
    4042        introduction/exp_apx
    4143        introduction/get_started
  • trunk/bin/svn_merge.sh

    r2506 r2859  
    3333#
    3434# Name of the directory where the changes have been committed
    35 from_branch='branches/sparse'
     35from_branch='branches/atomic'
    3636#
    3737# Version of the repository corresponding to from_branch just before changes
    38 Start=2474
     38Start=2797
    3939#
    4040# Version of the repository corresponding to from_branch after the changes
    41 End=2479
     41End=2858
    4242#
    4343# the svn merge command
  • trunk/cmake/add_cppad_cxx_flags.cmake

    r2772 r2859  
    3232                ENDIF( ${target_name} MATCHES ".*_${package}$" )
    3333        ENDFOREACH(package)
    34         SET_TARGET_PROPERTIES(
    35                 ${target_name} PROPERTIES COMPILE_FLAGS "${flags}"
    36         )
     34        IF( flags )
     35                SET_TARGET_PROPERTIES(
     36                        ${target_name} PROPERTIES COMPILE_FLAGS "${flags}"
     37                )
     38        ELSE( flags )
     39                SET_TARGET_PROPERTIES(
     40                        ${target_name} PROPERTIES COMPILE_FLAGS ""
     41                )
     42        ENDIF( flags )
    3743ENDMACRO(add_cppad_cxx_flags)
  • trunk/configure

    r2792 r2859  
    11#! /bin/sh
    22# Guess values for system-dependent variables and create Makefiles.
    3 # Generated by GNU Autoconf 2.68 for cppad 20130428.
     3# Generated by GNU Autoconf 2.68 for cppad 20130527.
    44#
    55# Report bugs to <cppad@list.coin-or.org>.
     
    561561PACKAGE_NAME='cppad'
    562562PACKAGE_TARNAME='cppad'
    563 PACKAGE_VERSION='20130428'
    564 PACKAGE_STRING='cppad 20130428'
     563PACKAGE_VERSION='20130527'
     564PACKAGE_STRING='cppad 20130527'
    565565PACKAGE_BUGREPORT='cppad@list.coin-or.org'
    566566PACKAGE_URL=''
     
    13691369  # This message is too long to be a string in the A/UX 3.1 sh.
    13701370  cat <<_ACEOF
    1371 \`configure' configures cppad 20130428 to adapt to many kinds of systems.
     1371\`configure' configures cppad 20130527 to adapt to many kinds of systems.
    13721372
    13731373Usage: $0 [OPTION]... [VAR=VALUE]...
     
    14351435if test -n "$ac_init_help"; then
    14361436  case $ac_init_help in
    1437      short | recursive ) echo "Configuration of cppad 20130428:";;
     1437     short | recursive ) echo "Configuration of cppad 20130527:";;
    14381438   esac
    14391439  cat <<\_ACEOF
     
    15631563if $ac_init_version; then
    15641564  cat <<\_ACEOF
    1565 cppad configure 20130428
     1565cppad configure 20130527
    15661566generated by GNU Autoconf 2.68
    15671567
     
    21872187running configure, to aid debugging if configure makes a mistake.
    21882188
    2189 It was created by cppad $as_me 20130428, which was
     2189It was created by cppad $as_me 20130527, which was
    21902190generated by GNU Autoconf 2.68.  Invocation command line was
    21912191
     
    50345034# Define the identity of the package.
    50355035 PACKAGE='cppad'
    5036  VERSION='20130428'
     5036 VERSION='20130527'
    50375037
    50385038
     
    74367436ipopt_prefix=${IPOPT_DIR}
    74377437
    7438 ac_config_files="$ac_config_files cppad/configure.hpp cppad_ipopt/example/test.sh cppad_ipopt/speed/test.sh cppad_ipopt/test/test.sh example/ipopt_solve/test.sh example/test_one.sh pkgconfig/cppad.pc pkgconfig/cppad-uninstalled.pc test_more/test_one.sh makefile compare_c/makefile cppad_ipopt/src/makefile cppad_ipopt/example/makefile cppad_ipopt/speed/makefile cppad_ipopt/test/makefile example/makefile example/ipopt_solve/makefile introduction/get_started/makefile introduction/exp_apx/makefile multi_thread/makefile multi_thread/test_multi/makefile print_for/makefile speed/adolc/makefile speed/cppad/makefile speed/double/makefile speed/example/makefile speed/fadbad/makefile speed/profile/makefile speed/profile/gprof.sed speed/sacado/makefile speed/src/makefile test_more/makefile"
     7438ac_config_files="$ac_config_files cppad/configure.hpp cppad_ipopt/example/test.sh cppad_ipopt/speed/test.sh cppad_ipopt/test/test.sh example/ipopt_solve/test.sh example/test_one.sh pkgconfig/cppad.pc pkgconfig/cppad-uninstalled.pc test_more/test_one.sh makefile compare_c/makefile cppad_ipopt/src/makefile cppad_ipopt/example/makefile cppad_ipopt/speed/makefile cppad_ipopt/test/makefile example/makefile example/atomic/makefile example/ipopt_solve/makefile introduction/get_started/makefile introduction/exp_apx/makefile multi_thread/makefile multi_thread/test_multi/makefile print_for/makefile speed/adolc/makefile speed/cppad/makefile speed/double/makefile speed/example/makefile speed/fadbad/makefile speed/profile/makefile speed/profile/gprof.sed speed/sacado/makefile speed/src/makefile test_more/makefile"
    74397439
    74407440
     
    80868086# values after options handling.
    80878087ac_log="
    8088 This file was extended by cppad $as_me 20130428, which was
     8088This file was extended by cppad $as_me 20130527, which was
    80898089generated by GNU Autoconf 2.68.  Invocation command line was
    80908090
     
    81438143ac_cs_config="`$as_echo "$ac_configure_args" | sed 's/^ //; s/[\\""\`\$]/\\\\&/g'`"
    81448144ac_cs_version="\\
    8145 cppad config.status 20130428
     8145cppad config.status 20130527
    81468146configured by $0, generated by GNU Autoconf 2.68,
    81478147  with options \\"\$ac_cs_config\\"
     
    82868286    "cppad_ipopt/test/makefile") CONFIG_FILES="$CONFIG_FILES cppad_ipopt/test/makefile" ;;
    82878287    "example/makefile") CONFIG_FILES="$CONFIG_FILES example/makefile" ;;
     8288    "example/atomic/makefile") CONFIG_FILES="$CONFIG_FILES example/atomic/makefile" ;;
    82888289    "example/ipopt_solve/makefile") CONFIG_FILES="$CONFIG_FILES example/ipopt_solve/makefile" ;;
    82898290    "introduction/get_started/makefile") CONFIG_FILES="$CONFIG_FILES introduction/get_started/makefile" ;;
  • trunk/configure.ac

    r2792 r2859  
    1313dnl Process this file with autoconf to produce a configure script.
    1414dnl   package   version              bug-report
    15 AC_INIT(cppad, 20130428, cppad@list.coin-or.org)
     15AC_INIT(cppad, 20130527, cppad@list.coin-or.org)
    1616AM_SILENT_RULES([yes])
    1717
     
    726726        cppad_ipopt/test/makefile
    727727        example/makefile
     728        example/atomic/makefile
    728729        example/ipopt_solve/makefile
    729730        introduction/get_started/makefile
  • trunk/cppad/local/abs_op.hpp

    r2625 r2859  
    44
    55/* --------------------------------------------------------------------------
    6 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
     6CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-13 Bradley M. Bell
    77
    88CppAD is distributed under multiple licenses. This distribution is under
     
    3535template <class Base>
    3636inline void forward_abs_op(
    37         size_t j           ,
     37        size_t q           ,
     38        size_t p           ,
    3839        size_t i_z         ,
    3940        size_t i_x         ,
     
    4546        CPPAD_ASSERT_UNKNOWN( NumRes(AbsOp) == 1 );
    4647        CPPAD_ASSERT_UNKNOWN( i_x < i_z );
    47         CPPAD_ASSERT_UNKNOWN( j < nc_taylor );
     48        CPPAD_ASSERT_UNKNOWN( p < nc_taylor );
     49        CPPAD_ASSERT_UNKNOWN( q <= p );
    4850
    4951        // Taylor coefficients corresponding to argument and result
     
    5153        Base* z = taylor + i_z * nc_taylor;
    5254
    53         z[j] = sign(x[0]) * x[j];
     55        for(size_t j = q; j <= p; j++)
     56                z[j] = sign(x[0]) * x[j];
    5457}
    5558
  • trunk/cppad/local/acos_op.hpp

    r2794 r2859  
    4242template <class Base>
    4343inline void forward_acos_op(
    44         size_t j           ,
     44        size_t q           ,
     45        size_t p           ,
    4546        size_t i_z         ,
    4647        size_t i_x         ,
     
    5253        CPPAD_ASSERT_UNKNOWN( NumRes(AcosOp) == 2 );
    5354        CPPAD_ASSERT_UNKNOWN( i_x + 1 < i_z );
    54         CPPAD_ASSERT_UNKNOWN( j < nc_taylor );
     55        CPPAD_ASSERT_UNKNOWN( p < nc_taylor );
     56        CPPAD_ASSERT_UNKNOWN( q <= p );
    5557
    5658        // Taylor coefficients corresponding to argument and result
     
    6163        size_t k;
    6264        Base uj;
    63         if( j == 0 )
     65        if( q == 0 )
    6466        {       z[0] = acos( x[0] );
    6567                uj   = Base(1) - x[0] * x[0];
    6668                b[0] = sqrt( uj );
     69                q++;
    6770        }
    68         else
     71        for(size_t j = q; j <= p; j++)
    6972        {       uj = 0.;
    7073                for(k = 0; k <= j; k++)
  • trunk/cppad/local/ad.hpp

    r2692 r2859  
    44
    55/* --------------------------------------------------------------------------
    6 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
     6CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-13 Bradley M. Bell
    77
    88CppAD is distributed under multiple licenses. This distribution is under
     
    101101        friend class ADTape<Base>;
    102102        friend class ADFun<Base>;
     103        friend class atomic_base<Base>;
    103104        friend class discrete<Base>;
    104         friend class user_atomic<Base>;
    105105        friend class VecAD<Base>;
    106106        friend class VecAD_reference<Base>;
  • trunk/cppad/local/ad_fun.hpp

    r2625 r2859  
    44
    55/* --------------------------------------------------------------------------
    6 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
     6CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-13 Bradley M. Bell
    77
    88CppAD is distributed under multiple licenses. This distribution is under
     
    127127        void ForSparseJacCase(
    128128                bool               set_type  ,
     129                bool               transpose ,
    129130                size_t             q         ,
    130131                const VectorSet&   r         , 
     
    136137        void ForSparseJacCase(
    137138                const std::set<size_t>&  set_type  ,
     139                bool                     transpose ,
    138140                size_t                   q         ,
    139141                const VectorSet&         r         , 
     
    146148        void RevSparseJacCase(
    147149                bool               set_type  ,
     150                bool               transpose ,
    148151                size_t             p         ,
    149152                const VectorSet&   s         , 
     
    155158        void RevSparseJacCase(
    156159                const std::set<size_t>&  set_type  ,
     160                bool                     transpose ,
    157161                size_t                   p         ,
    158162                const VectorSet&         s         , 
     
    165169        void RevSparseHesCase(
    166170                bool               set_type  ,
     171                bool               transpose ,
    167172                size_t             q         ,
    168173                const VectorSet&   s         , 
     
    174179        void RevSparseHesCase(
    175180                const std::set<size_t>&  set_type  ,
     181                bool                     transpose ,
    176182                size_t                   q         ,
    177183                const VectorSet&         s         , 
     
    341347        template <typename VectorSet>
    342348        VectorSet ForSparseJac(
    343                 size_t q, const VectorSet &r
     349                size_t q, const VectorSet &r, bool transpose = false
    344350        );
    345351        // reverse mode Jacobian sparsity
     
    347353        template <typename VectorSet>
    348354        VectorSet RevSparseJac(
    349                 size_t q, const VectorSet &s
     355                size_t q, const VectorSet &s, bool transpose = false
    350356        );
    351357        // reverse mode Hessian sparsity
     
    353359        template <typename VectorSet>
    354360        VectorSet RevSparseHes(
    355                 size_t q, const VectorSet &s
     361                size_t q, const VectorSet &s, bool transpose = false
    356362        );
    357363
  • trunk/cppad/local/ad_tape.hpp

    r2506 r2859  
    44
    55/* --------------------------------------------------------------------------
    6 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
     6CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-13 Bradley M. Bell
    77
    88CppAD is distributed under multiple licenses. This distribution is under
     
    3131        friend class AD<Base>;
    3232        friend class ADFun<Base>;
     33        friend class atomic_base<Base>;
    3334        friend class discrete<Base>;
    34         friend class user_atomic<Base>;
    3535        friend class VecAD<Base>;
    3636        friend class VecAD_reference<Base>;
  • trunk/cppad/local/ad_valued.hpp

    r2506 r2859  
    44
    55/* --------------------------------------------------------------------------
    6 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-11 Bradley M. Bell
     6CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-13 Bradley M. Bell
    77
    88CppAD is distributed under multiple licenses. This distribution is under
     
    2424$section AD Valued Operations and Functions$$
    2525
     26$comment atomic.omh includes atomic_base.omh which atomic_base.hpp$$
    2627$childtable%
    2728        cppad/local/arithmetic.hpp%
     
    3031        cppad/local/cond_exp.hpp%
    3132        cppad/local/discrete.hpp%
    32         cppad/local/user_atomic.hpp
     33        omh/atomic.omh
    3334%$$
    3435
     
    4445# include <cppad/local/math_other.hpp>
    4546# include <cppad/local/discrete.hpp>
    46 # include <cppad/local/user_atomic.hpp>
     47# include <cppad/local/atomic_base.hpp>
     48# include <cppad/local/checkpoint.hpp>
     49# include <cppad/local/old_atomic.hpp>
    4750
    4851# endif
  • trunk/cppad/local/add_op.hpp

    r2625 r2859  
    44
    55/* --------------------------------------------------------------------------
    6 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
     6CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-13 Bradley M. Bell
    77
    88CppAD is distributed under multiple licenses. This distribution is under
     
    3939template <class Base>
    4040inline void forward_addvv_op(
    41         size_t        d           ,
     41        size_t        q           ,
     42        size_t        p           ,
    4243        size_t        i_z         ,
    4344        const addr_t* arg         ,
     
    5152        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < i_z );
    5253        CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < i_z );
    53         CPPAD_ASSERT_UNKNOWN( d < nc_taylor );
     54        CPPAD_ASSERT_UNKNOWN( p < nc_taylor );
     55        CPPAD_ASSERT_UNKNOWN( q <= p  );
    5456
    5557        // Taylor coefficients corresponding to arguments and result
     
    5860        Base* z = taylor + i_z    * nc_taylor;
    5961
    60         z[d] = x[d] + y[d];
     62        for(size_t j = q; j <= p; j++)
     63                z[j] = x[j] + y[j];
    6164}
    6265
     
    160163template <class Base>
    161164inline void forward_addpv_op(
    162         size_t        d           ,
     165        size_t        q           ,
     166        size_t        p           ,
    163167        size_t        i_z         ,
    164168        const addr_t* arg         ,
     
    171175        CPPAD_ASSERT_UNKNOWN( NumRes(AddpvOp) == 1 );
    172176        CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < i_z );
    173         CPPAD_ASSERT_UNKNOWN( d < nc_taylor );
    174 
    175         // Taylor coefficients corresponding to arguments and result
    176         Base* y = taylor + arg[1] * nc_taylor;
    177         Base* z = taylor + i_z    * nc_taylor;
    178 
    179 # if CPPAD_USE_FORWARD0SWEEP
    180         CPPAD_ASSERT_UNKNOWN( d > 0 );
    181         z[d] = y[d];
    182 # else
    183         // Paraemter value
    184         Base x = parameter[ arg[0] ];
    185         if( d == 0 )
    186                 z[d] = x + y[d];
    187         else    z[d] = y[d];
    188 # endif
     177        CPPAD_ASSERT_UNKNOWN( p < nc_taylor );
     178        CPPAD_ASSERT_UNKNOWN( q <= p );
     179
     180        // Taylor coefficients corresponding to arguments and result
     181        Base* y = taylor + arg[1] * nc_taylor;
     182        Base* z = taylor + i_z    * nc_taylor;
     183
     184        if( q == 0 )
     185        {       // Paraemter value
     186                Base x = parameter[ arg[0] ];
     187                z[0] = x + y[0];
     188                q++;
     189        }
     190        for(size_t j = q; j <= p; j++)
     191                z[j] = y[j];
    189192}
    190193/*!
  • trunk/cppad/local/asin_op.hpp

    r2794 r2859  
    4242template <class Base>
    4343inline void forward_asin_op(
    44         size_t j           ,
     44        size_t q           ,
     45        size_t p           ,
    4546        size_t i_z         ,
    4647        size_t i_x         ,
     
    5253        CPPAD_ASSERT_UNKNOWN( NumRes(AsinOp) == 2 );
    5354        CPPAD_ASSERT_UNKNOWN( i_x + 1 < i_z );
    54         CPPAD_ASSERT_UNKNOWN( j < nc_taylor );
     55        CPPAD_ASSERT_UNKNOWN( p < nc_taylor );
     56        CPPAD_ASSERT_UNKNOWN( q <= p );
    5557
    5658        // Taylor coefficients corresponding to argument and result
     
    6163        size_t k;
    6264        Base uj;
    63         if( j == 0 )
     65        if( q == 0 )
    6466        {       z[0] = asin( x[0] );
    6567                uj   = Base(1) - x[0] * x[0];
    6668                b[0] = sqrt( uj );
     69                q++;
    6770        }
    68         else
     71        for(size_t j = q; j <= p; j++)
    6972        {       uj = 0.;
    7073                for(k = 0; k <= j; k++)
  • trunk/cppad/local/atan_op.hpp

    r2625 r2859  
    44
    55/* --------------------------------------------------------------------------
    6 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
     6CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-13 Bradley M. Bell
    77
    88CppAD is distributed under multiple licenses. This distribution is under
     
    4242template <class Base>
    4343inline void forward_atan_op(
    44         size_t j           ,
     44        size_t q           ,
     45        size_t p           ,
    4546        size_t i_z         ,
    4647        size_t i_x         ,
     
    5253        CPPAD_ASSERT_UNKNOWN( NumRes(AtanOp) == 2 );
    5354        CPPAD_ASSERT_UNKNOWN( i_x + 1 < i_z );
    54         CPPAD_ASSERT_UNKNOWN( j < nc_taylor );
     55        CPPAD_ASSERT_UNKNOWN( p < nc_taylor );
     56        CPPAD_ASSERT_UNKNOWN( q <= p );
    5557
    5658        // Taylor coefficients corresponding to argument and result
     
    6062
    6163        size_t k;
    62         if( j == 0 )
    63         {       z[j] = atan( x[0] );
    64                 b[j] = Base(1) + x[0] * x[0];
     64        if( q == 0 )
     65        {       z[0] = atan( x[0] );
     66                b[0] = Base(1) + x[0] * x[0];
     67                q++;
    6568        }
    66         else
     69        for(size_t j = q; j <= p; j++)
    6770        {
    6871                b[j] = Base(2) * x[0] * x[j];
  • trunk/cppad/local/bender_quad.hpp

    r2506 r2859  
    44
    55/* --------------------------------------------------------------------------
    6 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
     6CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-13 Bradley M. Bell
    77
    88CppAD is distributed under multiple licenses. This distribution is under
     
    3737
    3838$head Syntax$$
    39 $icode%
    40 %# include <cppad/cppad.hpp>
     39$codei%
     40# include <cppad/cppad.hpp>
    4141BenderQuad(%x%, %y%, %fun%, %g%, %gx%, %gxx%)%$$ 
    4242
  • trunk/cppad/local/cond_op.hpp

    r2625 r2859  
    33# define CPPAD_COND_OP_INCLUDED
    44/* --------------------------------------------------------------------------
    5 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
     5CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-13 Bradley M. Bell
    66
    77CppAD is distributed under multiple licenses. This distribution is under
     
    2626\copydetails conditional_exp_op
    2727
    28 \param d
    29 is the order of the Taylor coefficient of z that we are computing.
     28\param q
     29is the lowest order of the Taylor coefficient of z that we are computing.
     30
     31\param p
     32is the highest order of the Taylor coefficient of z that we are computing.
    3033
    3134\param taylor
    3235\b Input:
    33 For j = 0, 1, 2, 3 and k = 0 , ... , \a d,
     36For j = 0, 1, 2, 3 and k = 0 , ... , p,
    3437if y_j is a variable then
    35 \a taylor [ \a arg[2+j] * nc_taylor + k ]
     38<code>taylor [ arg[2+j] * nc_taylor + k ]</code>
    3639is the k-th order Taylor coefficient corresponding to y_j.
    3740\n
    38 \b Input: \a taylor [ \a i_z * \a nc_taylor + k ]
    39 for k = 0 , ... , \a d - 1
     41\b Input: <code>taylor [ i_z * nc_taylor + k ]</code>
     42for k = 0 , ... , q-1,
    4043is the k-th order Taylor coefficient corresponding to z.
    4144\n
    42 \b Output: \a taylor [ \a i_z * \a nc_taylor + \a d ]
    43 is the \a d-th order Taylor coefficient corresponding to z.
     45\b Output: <code>taylor [ i_z * nc_taylor + k ]</code>
     46for k = q , ... , p,
     47is the k-th order Taylor coefficient corresponding to z.
    4448
    4549*/
    4650template <class Base>
    4751inline void forward_cond_op(
    48         size_t         d           ,
     52        size_t         q           ,
     53        size_t         p           ,
    4954        size_t         i_z         ,
    5055        const addr_t*  arg         ,
     
    5560{       Base y_0, y_1, y_2, y_3;
    5661        Base zero(0);
    57         Base* z;
     62        Base* z = taylor + i_z * nc_taylor;
    5863
    5964        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < static_cast<size_t> (CompareNe) );
     
    7883                y_1 = parameter[ arg[3] ];
    7984        }
    80 # if CPPAD_USE_FORWARD0SWEEP
    81         CPPAD_ASSERT_UNKNOWN( d > 0 );
    82 # else
    83         if( d == 0 )
     85        if( q == 0 )
    8486        {       if( arg[1] & 4 )
    8587                {       CPPAD_ASSERT_UNKNOWN( size_t(arg[4]) < i_z );
     
    98100                        y_3 = parameter[ arg[5] ];
    99101                }
    100         }
    101         else
    102 # endif
     102                z[0] = CondExpOp(
     103                        CompareOp( arg[0] ),
     104                        y_0,
     105                        y_1,
     106                        y_2,
     107                        y_3
     108                );
     109                q++;
     110        }
     111        for(size_t d = q; d <= p; d++)
    103112        {       if( arg[1] & 4 )
    104113                {       CPPAD_ASSERT_UNKNOWN( size_t(arg[4]) < i_z );
     
    111120                }
    112121                else    y_3 = zero;
    113         }
    114         z = taylor + i_z * nc_taylor;
    115         z[d] = CondExpOp(
    116                 CompareOp( arg[0] ),
    117                 y_0,
    118                 y_1,
    119                 y_2,
    120                 y_3
    121         );
     122                z[d] = CondExpOp(
     123                        CompareOp( arg[0] ),
     124                        y_0,
     125                        y_1,
     126                        y_2,
     127                        y_3
     128                );
     129        }
    122130        return;
    123131}
  • trunk/cppad/local/cos_op.hpp

    r2625 r2859  
    44
    55/* --------------------------------------------------------------------------
    6 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
     6CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-13 Bradley M. Bell
    77
    88CppAD is distributed under multiple licenses. This distribution is under
     
    4242template <class Base>
    4343inline void forward_cos_op(
    44         size_t j           ,
     44        size_t q           ,
     45        size_t p           ,
    4546        size_t i_z         ,
    4647        size_t i_x         ,
     
    5253        CPPAD_ASSERT_UNKNOWN( NumRes(CosOp) == 2 );
    5354        CPPAD_ASSERT_UNKNOWN( i_x + 1 < i_z );
    54         CPPAD_ASSERT_UNKNOWN( j < nc_taylor );
     55        CPPAD_ASSERT_UNKNOWN( p < nc_taylor );
     56        CPPAD_ASSERT_UNKNOWN( q <= p );
    5557
    5658        // Taylor coefficients corresponding to argument and result
     
    6365        // forward_sin_op, forward_cos_op, forward_sinh_op, forward_cosh_op.
    6466        size_t k;
    65         if( j == 0 )
    66         {       s[j] = sin( x[0] );
    67                 c[j] = cos( x[0] );
     67        if( q == 0 )
     68        {       s[0] = sin( x[0] );
     69                c[0] = cos( x[0] );
     70                q++;
    6871        }
    69         else
     72        for(size_t j = q; j <= p; j++)
    7073        {
    7174                s[j] = Base(0);
  • trunk/cppad/local/cosh_op.hpp

    r2625 r2859  
    44
    55/* --------------------------------------------------------------------------
    6 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
     6CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-13 Bradley M. Bell
    77
    88CppAD is distributed under multiple licenses. This distribution is under
     
    4242template <class Base>
    4343inline void forward_cosh_op(
    44         size_t j           ,
     44        size_t q           ,
     45        size_t p           ,
    4546        size_t i_z         ,
    4647        size_t i_x         ,
     
    5253        CPPAD_ASSERT_UNKNOWN( NumRes(CoshOp) == 2 );
    5354        CPPAD_ASSERT_UNKNOWN( i_x + 1 < i_z );
    54         CPPAD_ASSERT_UNKNOWN( j < nc_taylor );
     55        CPPAD_ASSERT_UNKNOWN( p < nc_taylor );
     56        CPPAD_ASSERT_UNKNOWN( q <= p );
    5557
    5658        // Taylor coefficients corresponding to argument and result
     
    6264        // forward_sin_op, forward_cos_op, forward_sinh_op, forward_cosh_op.
    6365        size_t k;
    64         if( j == 0 )
    65         {       s[j] = sinh( x[0] );
    66                 c[j] = cosh( x[0] );
     66        if( q == 0 )
     67        {       s[0] = sinh( x[0] );
     68                c[0] = cosh( x[0] );
     69                q++;
    6770        }
    68         else
     71        for(size_t j = q; j <= p; j++)
    6972        {
    7073                s[j] = Base(0);
  • trunk/cppad/local/csum_op.hpp

    r2625 r2859  
    44
    55/* --------------------------------------------------------------------------
    6 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
     6CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-13 Bradley M. Bell
    77
    88CppAD is distributed under multiple licenses. This distribution is under
     
    2727This operation is
    2828\verbatim
    29         z = p + x(1) + ... + x(m) - y(1) - ... - y(n).
     29        z = s + x(1) + ... + x(m) - y(1) - ... - y(n).
    3030\endverbatim
    3131
     
    3535\a Base.
    3636
    37 \param d
    38 order of the Taylor coefficient that we are computing.
     37\param q
     38lowest order of the Taylor coefficient that we are computing.
     39
     40\param p
     41highest order of the Taylor coefficient that we are computing.
    3942
    4043\param i_z
     
    5255\n
    5356<tt>parameter[ arg[2] ]</tt>
    54 is the parameter value \c p in this cummunative summation.
     57is the parameter value \c s in this cummunative summation.
    5558\n
    5659<tt>arg[2+i]</tt>
    57 for <tt>i = 1 , ... , m</tt> is the value <tt>x(i)</tt>.
     60for <tt>i = 1 , ... , m</tt> is the variable index of <tt>x(i)</tt>.
    5861\n
    5962<tt>arg[2+arg[0]+i]</tt>
    60 for <tt>i = 1 , ... , n</tt> is the value <tt>y(i)</tt>.
     63for <tt>i = 1 , ... , n</tt> is the variable index of <tt>y(i)</tt>.
    6164
    6265\param num_par
     
    7275\b Input: <tt>taylor [ arg[2+i] * nc_taylor + k ]</tt>
    7376for <tt>i = 1 , ... , m</tt>
    74 and <tt>k = 0 , ... , d</tt>
     77and <tt>k = 0 , ... , p</tt>
    7578is the k-th order Taylor coefficient corresponding to <tt>x(i)</tt>
    7679\n
    7780\b Input: <tt>taylor [ arg[2+m+i] * nc_taylor + k ]</tt>
    7881for <tt>i = 1 , ... , n</tt>
    79 and <tt>k = 0 , ... , d</tt>
     82and <tt>k = 0 , ... , p</tt>
    8083is the k-th order Taylor coefficient corresponding to <tt>y(i)</tt>
    8184\n
    8285\b Input: <tt>taylor [ i_z * nc_taylor + k ]</tt>
    83 for k = 0 , ... , \a d - 1
     86for k = 0 , ... , q,
    8487is the k-th order Taylor coefficient corresponding to z.
    8588\n
    86 \b Output: <tt>taylor [ i_z * nc_taylor + d ]</tt>
    87 is the \a d-th order Taylor coefficient corresponding to z.
     89\b Output: <tt>taylor [ i_z * nc_taylor + k ]</tt>
     90for k = q , ... , p,
     91is the \a k-th order Taylor coefficient corresponding to z.
    8892*/
    8993template <class Base>
    9094inline void forward_csum_op(
    91         size_t        d           ,
     95        size_t        q           ,
     96        size_t        p           ,
    9297        size_t        i_z         ,
    9398        const addr_t* arg         ,
     
    97102        Base*         taylor      )
    98103{       Base zero(0);
     104        size_t i, j, k;
    99105
    100106        // check assumptions
    101107        CPPAD_ASSERT_UNKNOWN( NumRes(CSumOp) == 1 );
    102         CPPAD_ASSERT_UNKNOWN( d < nc_taylor );
     108        CPPAD_ASSERT_UNKNOWN( p < nc_taylor );
     109        CPPAD_ASSERT_UNKNOWN( q <= p );
    103110        CPPAD_ASSERT_UNKNOWN( size_t(arg[2]) < num_par );
    104111        CPPAD_ASSERT_UNKNOWN(
     
    108115        // Taylor coefficients corresponding to result
    109116        Base* z = taylor + i_z    * nc_taylor;
    110         if( d == 0 )
    111                 z[d] = parameter[ arg[2] ];
    112         else    z[d] = zero;
     117        for(k = q; k <= p; k++)
     118                z[q] = zero;
     119        if( q == 0 )
     120                z[q] = parameter[ arg[2] ];
    113121        Base* x;
    114         size_t i, j;
    115122        i = arg[0];
    116123        j = 2;
     
    118125        {       CPPAD_ASSERT_UNKNOWN( size_t(arg[j+1]) < i_z );
    119126                x     = taylor + arg[++j] * nc_taylor;
    120                 z[d] += x[d];
     127                for(k = q; k <= p; k++)
     128                        z[k] += x[k];
    121129        }       
    122130        i = arg[1];
     
    124132        {       CPPAD_ASSERT_UNKNOWN( size_t(arg[j+1]) < i_z );
    125133                x     = taylor + arg[++j] * nc_taylor;
    126                 z[d] -= x[d];
     134                for(k = q; k <= p; k++)
     135                        z[k] -= x[k];
    127136        }       
    128137}
  • trunk/cppad/local/declare_ad.hpp

    r2692 r2859  
    44
    55/* --------------------------------------------------------------------------
    6 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
     6CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-13 Bradley M. Bell
    77
    88CppAD is distributed under multiple licenses. This distribution is under
     
    3535        template <class Base> class ADFun;
    3636        template <class Base> class ADTape;
     37        template <class Base> class atomic_base;
    3738        template <class Base> class discrete;
    3839        template <class Base> class player;
    3940        template <class Base> class recorder;
    40         template <class Base> class user_atomic;
    4141        template <class Base> class VecAD;
    4242        template <class Base> class VecAD_reference;
  • trunk/cppad/local/div_op.hpp

    r2794 r2859  
    3939template <class Base>
    4040inline void forward_divvv_op(
    41         size_t        d           ,
     41        size_t        q           ,
     42        size_t        p           ,
    4243        size_t        i_z         ,
    4344        const addr_t* arg         ,
     
    5152        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < i_z );
    5253        CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < i_z );
    53         CPPAD_ASSERT_UNKNOWN( d < nc_taylor );
     54        CPPAD_ASSERT_UNKNOWN( p < nc_taylor );
     55        CPPAD_ASSERT_UNKNOWN( q <= p );
    5456
    5557        // Taylor coefficients corresponding to arguments and result
     
    6264        // so do not make it an error.
    6365        size_t k;
    64         z[d] = x[d];
    65         for(k = 1; k <= d; k++)
    66                 z[d] -= z[d-k] * y[k];
    67         z[d] /= y[0];
     66        for(size_t d = q; d <= p; d++)
     67        {       z[d] = x[d];
     68                for(k = 1; k <= d; k++)
     69                        z[d] -= z[d-k] * y[k];
     70                z[d] /= y[0];
     71        }
    6872}
    6973
     
    183187template <class Base>
    184188inline void forward_divpv_op(
    185         size_t        d           ,
     189        size_t        q           ,
     190        size_t        p           ,
    186191        size_t        i_z         ,
    187192        const addr_t* arg         ,
     
    194199        CPPAD_ASSERT_UNKNOWN( NumRes(DivpvOp) == 1 );
    195200        CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < i_z );
    196         CPPAD_ASSERT_UNKNOWN( d < nc_taylor );
     201        CPPAD_ASSERT_UNKNOWN( p < nc_taylor );
     202        CPPAD_ASSERT_UNKNOWN( q <= p );
    197203
    198204        // Taylor coefficients corresponding to arguments and result
     
    206212        // so do not make it an error.
    207213        size_t k;
    208 # if USE_CPPAD_FORWARD0SWEEP
    209         z[d] = Base(0);
    210 # else
    211         if( d == 0 )
    212                 z[d] = x;
    213         else    z[d] = Base(0);
    214 # endif
    215         for(k = 1; k <= d; k++)
    216                 z[d] -= z[d-k] * y[k];
    217         z[d] /= y[0];
    218 
     214        if( q == 0 )
     215        {       z[0] = x / y[0];
     216                q++;
     217        }
     218        for(size_t d = q; d <= p; d++)
     219        {       z[d] = Base(0);
     220                for(k = 1; k <= d; k++)
     221                        z[d] -= z[d-k] * y[k];
     222                z[d] /= y[0];
     223        }
    219224}
    220225
     
    330335template <class Base>
    331336inline void forward_divvp_op(
    332         size_t        d           ,
     337        size_t        q           ,
     338        size_t        p           ,
    333339        size_t        i_z         ,
    334340        const addr_t* arg         ,
     
    341347        CPPAD_ASSERT_UNKNOWN( NumRes(DivvpOp) == 1 );
    342348        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < i_z );
    343         CPPAD_ASSERT_UNKNOWN( d < nc_taylor );
     349        CPPAD_ASSERT_UNKNOWN( p < nc_taylor );
     350        CPPAD_ASSERT_UNKNOWN( q <= p );
    344351
    345352        // Taylor coefficients corresponding to arguments and result
     
    352359        // Using CondExp and multiple levels of AD, it can make sense
    353360        // to divide by zero so do not make it an error.
    354         z[d] = x[d] / y;
     361        for(size_t d = q; d <= p; d++)
     362                z[d] = x[d] / y;
    355363}
    356364
  • trunk/cppad/local/exp_op.hpp

    r2625 r2859  
    44
    55/* --------------------------------------------------------------------------
    6 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
     6CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-13 Bradley M. Bell
    77
    88CppAD is distributed under multiple licenses. This distribution is under
     
    3636template <class Base>
    3737inline void forward_exp_op(
    38         size_t j           ,
     38        size_t q           ,
     39        size_t p           ,
    3940        size_t i_z         ,
    4041        size_t i_x         ,
     
    4647        CPPAD_ASSERT_UNKNOWN( NumRes(ExpOp) == 1 );
    4748        CPPAD_ASSERT_UNKNOWN( i_x < i_z );
    48         CPPAD_ASSERT_UNKNOWN( j < nc_taylor );
     49        CPPAD_ASSERT_UNKNOWN( p < nc_taylor );
     50        CPPAD_ASSERT_UNKNOWN( q <= p );
    4951
    5052        // Taylor coefficients corresponding to argument and result
     
    5355
    5456        size_t k;
    55         if( j == 0 )
    56                 z[0] = exp( x[0] );
    57         else
     57        if( q == 0 )
     58        {       z[0] = exp( x[0] );
     59                q++;
     60        }
     61        for(size_t j = q; j <= p; j++)
    5862        {
    5963                z[j] = x[1] * z[j-1];
     
    129133        Base* pz       = partial + i_z * nc_partial;
    130134
    131         // lopp through orders in reverse
     135        // loop through orders in reverse
    132136        size_t j, k;
    133137        j = d;
  • trunk/cppad/local/for_jac_sweep.hpp

    r2794 r2859  
    1212Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
    1313-------------------------------------------------------------------------- */
     14
    1415# include <set>
    1516# include <cppad/local/pod_vector.hpp>
     
    3031*/
    3132# define CPPAD_FOR_JAC_SWEEP_TRACE 0
     33
     34/*
     35\def CPPAD_ATOMIC_CALL
     36This avoids warnings when NDEBUG is defined and user_ok is not used.
     37If \c NDEBUG is defined, this resolves to
     38\code
     39        user_atom->for_sparse_jac
     40\endcode
     41otherwise, it respolves to
     42\code
     43        user_ok = user_atom->for_sparse_jac
     44\endcode
     45This maco is undefined at the end of this file to facillitate is
     46use with a different definition in other files.
     47*/
     48# ifdef NDEBUG
     49# define CPPAD_ATOMIC_CALL user_atom->for_sparse_jac
     50# else
     51# define CPPAD_ATOMIC_CALL user_ok = user_atom->for_sparse_jac
     52# endif
    3253
    3354/*!
     
    130151        }
    131152
     153        // --------------------------------------------------------------
    132154        // work space used by UserOp.
     155        //
    133156        typedef std::set<size_t> size_set;
    134         const size_t user_q = limit; // maximum element plus one
    135157        size_set::iterator set_itr;  // iterator for a standard set
    136158        size_set::iterator set_end;  // end of iterator sequence
    137         vector< size_set > user_r;   // sparsity pattern for the argument x
    138         vector< size_set > user_s;   // sparisty pattern for the result y
    139         size_t user_index = 0;       // indentifier for this user_atomic operation
     159        vector< size_set > set_r;    // set sparsity pattern for the argument x
     160        vector< size_set > set_s;    // set sparisty pattern for the result y
     161        //
     162        vector<bool>       bool_r;   // bool sparsity pattern for the argument x
     163        vector<bool>       bool_s;   // bool sparisty pattern for the result y
     164        //
     165        const size_t user_q = limit; // maximum element plus one
     166        size_t user_index = 0;       // indentifier for this atomic operation
    140167        size_t user_id    = 0;       // user identifier for this call to operator
    141168        size_t user_i     = 0;       // index in result vector
     
    143170        size_t user_m     = 0;       // size of result vector
    144171        size_t user_n     = 0;       // size of arugment vector
     172        //
     173        atomic_base<Base>* user_atom = CPPAD_NULL; // user's atomic op calculator
     174        bool               user_bool = false;      // use bool or set sparsity ?
     175# ifndef NDEBUG
     176        bool               user_ok   = false;      // atomic op return value
     177# endif
     178        //
    145179        // next expected operator in a UserOp sequence
    146180        enum { user_start, user_arg, user_ret, user_end } user_state = user_start;
     181        // --------------------------------------------------------------
    147182
    148183# if CPPAD_FOR_JAC_SWEEP_TRACE
     
    512547                                user_n     = arg[2];
    513548                                user_m     = arg[3];
    514                                 if(user_r.size() < user_n )
    515                                         user_r.resize(user_n);
    516                                 if(user_s.size() < user_m )
    517                                         user_s.resize(user_m);
     549                                user_atom  = atomic_base<Base>::list(user_index);
     550                                user_bool  = user_atom->sparsity() ==
     551                                                        atomic_base<Base>::bool_sparsity_enum;
     552                                if( user_bool )
     553                                {       if( bool_r.size() != user_n * user_q )
     554                                                bool_r.resize( user_n * user_q );
     555                                        if( bool_s.size() != user_m * user_q )
     556                                                bool_s.resize( user_m * user_q );
     557                                        for(i = 0; i < user_n; i++)
     558                                                for(j = 0; j < user_q; j++)
     559                                                        bool_r[ i * user_q + j] = false;
     560                                }
     561                                else
     562                                {       if(set_r.size() != user_n )
     563                                                set_r.resize(user_n);
     564                                        if(set_s.size() != user_m )
     565                                                set_s.resize(user_m);
     566                                        for(i = 0; i < user_n; i++)
     567                                                set_r[i].clear();
     568                                }
    518569                                user_j     = 0;
    519570                                user_i     = 0;
     
    526577                                CPPAD_ASSERT_UNKNOWN( user_n     == size_t(arg[2]) );
    527578                                CPPAD_ASSERT_UNKNOWN( user_m     == size_t(arg[3]) );
     579# ifndef NDEBUG
     580                                if( ! user_ok )
     581                                {       std::string msg = user_atom->afun_name()
     582                                        + ": atomic_base.for_sparse_jac: returned false";
     583                                        CPPAD_ASSERT_KNOWN(false, msg.c_str() );
     584                                }
     585# endif
    528586                                user_state = user_start;
    529587                        }
     
    535593                        CPPAD_ASSERT_UNKNOWN( user_j < user_n );
    536594                        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_par );
    537                         // parameters have an empty sparsity pattern
    538                         user_r[user_j].clear();
     595                        // set row user_j to empty sparsity pattern
    539596                        ++user_j;
    540597                        if( user_j == user_n )
    541598                        {       // call users function for this operation
    542                                 user_atomic<Base>::for_jac_sparse(user_index, user_id,
    543                                         user_n, user_m, user_q, user_r, user_s
     599                                user_atom->set_id(user_id);
     600                                if( user_bool )
     601                                        CPPAD_ATOMIC_CALL(
     602                                                user_q, bool_r, bool_s
     603                                );
     604                                else
     605                                        CPPAD_ATOMIC_CALL(
     606                                                user_q, set_r, set_s
    544607                                );
    545608                                user_state = user_ret;
     
    552615                        CPPAD_ASSERT_UNKNOWN( user_j < user_n );
    553616                        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) <= i_var );
    554                         // set user_r[user_j] to sparsity pattern for variable arg[0]
    555                         user_r[user_j].clear();
     617                        // set row user_j to sparsity pattern for variable arg[0]
    556618                        var_sparsity.begin(arg[0]);
    557619                        i = var_sparsity.next_element();
    558620                        while( i < user_q )
    559                         {       user_r[user_j].insert(i);
     621                        {       if( user_bool )
     622                                        bool_r[user_j * user_q + i] = true;
     623                                else
     624                                        set_r[user_j].insert(i);
    560625                                i = var_sparsity.next_element();
    561626                        }
     
    563628                        if( user_j == user_n )
    564629                        {       // call users function for this operation
    565                                 user_atomic<Base>::for_jac_sparse(user_index, user_id,
    566                                         user_n, user_m, user_q, user_r, user_s
     630                                user_atom->set_id(user_id);
     631                                if( user_bool )
     632                                        CPPAD_ATOMIC_CALL(
     633                                                user_q, bool_r, bool_s
     634                                );
     635                                else
     636                                        CPPAD_ATOMIC_CALL(
     637                                                user_q, set_r, set_s
    567638                                );
    568639                                user_state = user_ret;
     
    585656                        // It might be faster if we add set union to var_sparsity
    586657                        // where one of the sets is not in var_sparsity
    587                         set_itr = user_s[user_i].begin();
    588                         set_end = user_s[user_i].end();
    589                         while( set_itr != set_end )
    590                                 var_sparsity.add_element(i_var, *set_itr++);
     658                        if( user_bool )
     659                        {       for(j = 0; j < user_q; j++)
     660                                        if( bool_s[ user_i * user_q + j ] )
     661                                                var_sparsity.add_element(i_var, j);
     662                        }
     663                        else
     664                        {       set_itr = set_s[user_i].begin();
     665                                set_end = set_s[user_i].end();
     666                                while( set_itr != set_end )
     667                                        var_sparsity.add_element(i_var, *set_itr++);
     668                        }
    591669                        user_i++;
    592670                        if( user_i == user_m )
     
    634712// preprocessor symbols that are local to this file
    635713# undef CPPAD_FOR_JAC_SWEEP_TRACE
     714# undef CPPAD_ATOMIC_CALL
    636715
    637716# endif
  • trunk/cppad/local/for_sparse_jac.hpp

    r2625 r2859  
    44
    55/* --------------------------------------------------------------------------
    6 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
     6CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-13 Bradley M. Bell
    77
    88CppAD is distributed under multiple licenses. This distribution is under
     
    3535
    3636$head Syntax$$
    37 $icode%s% = %f%.ForSparseJac(%q%, %r%)%$$
     37$icode%s% = %f%.ForSparseJac(%q%, %r%)
     38%$$
     39$icode%s% = %f%.ForSparseJac(%q%, %r%, %transpose%)%$$
    3840
    3941$head Purpose$$
     
    117119$latex S(x) \in B^{m \times q}$$.
    118120
     121$head transpose$$
     122The argument $icode transpose$$ has prototype
     123$codei%
     124        bool %transpose%
     125%$$
     126The default value $code false$$ is used when $icode transpose$$ is not present.
     127
     128
    119129$head r$$
    120130The argument $icode r$$ has prototype
     
    122132        const %VectorSet%& %r%
    123133%$$
    124 (see $cref/VectorSet/ForSparseJac/VectorSet/$$ below).
    125 If it has elements of type $code bool$$,
     134see $cref/VectorSet/ForSparseJac/VectorSet/$$ below.
     135
     136$subhead transpose false$$
     137If $icode r$$ has elements of type $code bool$$,
    126138its size is $latex n * q$$.
    127139If it has elements of type $code std::set<size_t>$$,
     
    130142It specifies a
    131143$cref/sparsity pattern/glossary/Sparsity Pattern/$$
    132 for the matrix $icode R$$.
     144for the matrix $icode R \in B^{n \times q}$$.
     145
     146$subhead transpose true$$
     147If $icode r$$ has elements of type $code bool$$,
     148its size is $latex q * n$$.
     149If it has elements of type $code std::set<size_t>$$,
     150its size is $latex q$$ and all the set elements must be between
     151zero and $icode%n%-1%$$ inclusive.
     152It specifies a
     153$cref/sparsity pattern/glossary/Sparsity Pattern/$$
     154for the matrix $icode R^\R{T} \in B^{q \times n}$$.
    133155
    134156$head s$$
     
    137159        %VectorSet% %s%
    138160%$$
    139 (see $cref/VectorSet/ForSparseJac/VectorSet/$$ below).
    140 If it has elements of type $code bool$$,
     161see $cref/VectorSet/ForSparseJac/VectorSet/$$ below.
     162
     163$subhead transpose false$$
     164If $icode s$$ has elements of type $code bool$$,
    141165its size is $latex m * q$$.
    142166If it has elements of type $code std::set<size_t>$$,
     
    145169It specifies a
    146170$cref/sparsity pattern/glossary/Sparsity Pattern/$$
    147 for the matrix $latex S(x)$$.
     171for the matrix $latex S(x) \in B^{m \times q}$$.
     172
     173$subhead transpose true$$
     174If $icode s$$ has elements of type $code bool$$,
     175its size is $latex q * m$$.
     176If it has elements of type $code std::set<size_t>$$,
     177its size is $latex q$$ and all its set elements are between
     178zero and $icode%m%-1%$$ inclusive.
     179It specifies a
     180$cref/sparsity pattern/glossary/Sparsity Pattern/$$
     181for the matrix $latex S(x)^\R{T} \in B^{q \times m}$$.
    148182
    149183$head VectorSet$$
     
    198232\tparam VectorSet
    199233is a simple vector class with elements of type \c bool.
     234
     235\param transpose
     236are the sparsity patterns transposed.
    200237
    201238\param q
     
    241278template <class Base, class VectorSet>
    242279void ForSparseJacBool(
     280        bool                   transpose        ,
    243281        size_t                 q                ,
    244282        const VectorSet&       r                ,
     
    257295        size_t n = ind_taddr.size();
    258296
    259         CPPAD_ASSERT_UNKNOWN(q > 0 );
    260         CPPAD_ASSERT_UNKNOWN( size_t(r.size()) == n * q );
    261 
    262         // allocate memory for the requested sparsity calculation
     297        CPPAD_ASSERT_KNOWN(
     298                q > 0,
     299                "ForSparseJac: q is not greater than zero"
     300        );
     301        CPPAD_ASSERT_KNOWN(
     302                size_t(r.size()) == n * q,
     303                "ForSparseJac: size of r is not equal to\n"
     304                "q times domain dimension for ADFun object."
     305        );
     306
     307        // allocate memory for the requested sparsity calculation result
    263308        for_jac_sparsity.resize(total_num_var, q);
    264309
     
    270315
    271316                // set bits that are true
    272                 for(j = 0; j < q; j++) if( r[ i * q + j ] )
    273                         for_jac_sparsity.add_element( ind_taddr[i], j);
     317                if( transpose )
     318                {       for(j = 0; j < q; j++) if( r[ j * n + i ] )
     319                                for_jac_sparsity.add_element( ind_taddr[i], j);
     320                }
     321                else
     322                {       for(j = 0; j < q; j++) if( r[ i * q + j ] )
     323                                for_jac_sparsity.add_element( ind_taddr[i], j);
     324                }
    274325        }
    275326
     
    288339
    289340                // extract the result from for_jac_sparsity
    290                 for(j = 0; j < q; j++)
    291                         s[ i * q + j ] = false;
     341                if( transpose )
     342                {       for(j = 0; j < q; j++)
     343                                s[ j * m + i ] = false;
     344                }
     345                else
     346                {       for(j = 0; j < q; j++)
     347                                s[ i * q + j ] = false;
     348                }
    292349                CPPAD_ASSERT_UNKNOWN( for_jac_sparsity.end() == q );
    293350                for_jac_sparsity.begin( dep_taddr[i] );
    294351                j = for_jac_sparsity.next_element();
    295352                while( j < q )
    296                 {       s[i * q + j ] = true;
     353                {       if( transpose )
     354                                s[j * m + i] = true;
     355                        else    s[i * q + j] = true;
    297356                        j = for_jac_sparsity.next_element();
    298357                }
     
    314373is a simple vector class with elements of type \c std::set<size_t>.
    315374
     375\param transpose
     376see \c SparseJacBool.
     377
    316378\param q
    317379see \c SparseJacBool.
     
    341403template <class Base, class VectorSet>
    342404void ForSparseJacSet(
     405        bool                        transpose        ,
    343406        size_t                      q                ,
    344407        const VectorSet&            r                ,
     
    358421        size_t n = ind_taddr.size();
    359422
    360         CPPAD_ASSERT_UNKNOWN(q > 0 );
    361         CPPAD_ASSERT_UNKNOWN( size_t(r.size()) == n );
     423        CPPAD_ASSERT_KNOWN(
     424                q > 0,
     425                "RevSparseJac: q is not greater than zero"
     426        );
     427        CPPAD_ASSERT_KNOWN(
     428                size_t(r.size()) == n || transpose,
     429                "RevSparseJac: size of r is not equal to n and transpose is false."
     430        );
     431        CPPAD_ASSERT_KNOWN(
     432                size_t(r.size()) == q || ! transpose,
     433                "RevSparseJac: size of r is not equal to q and transpose is true."
     434        );
    362435
    363436        // allocate memory for the requested sparsity calculation
     
    365438
    366439        // set values corresponding to independent variables
    367         for(i = 0; i < n; i++)
    368         {       CPPAD_ASSERT_UNKNOWN( ind_taddr[i] < total_num_var );
    369                 // ind_taddr[i] is operator taddr for i-th independent variable
    370                 CPPAD_ASSERT_UNKNOWN( play.GetOp( ind_taddr[i] ) == InvOp );
    371 
    372                 // add the elements that are present
    373                 itr = r[i].begin();
    374                 while( itr != r[i].end() )
    375                 {       j = *itr++;
    376                         CPPAD_ASSERT_KNOWN(
    377                                 j < q,
    378                                 "ForSparseJac: an element of the set r[i] "
    379                                 "has value greater than or equal q."
    380                         );
    381                         for_jac_sparsity.add_element( ind_taddr[i], j);
     440        if( transpose )
     441        {       for(i = 0; i < q; i++)
     442                {       // add the elements that are present
     443                        itr = r[i].begin();
     444                        while( itr != r[i].end() )
     445                        {       j = *itr++;
     446                                CPPAD_ASSERT_KNOWN(
     447                                j < n,
     448                                "ForSparseJac: transpose is true and element of the set\n"
     449                                "r[j] has value greater than or equal n."
     450                                );
     451                                CPPAD_ASSERT_UNKNOWN( ind_taddr[j] < total_num_var );
     452                                // operator for j-th independent variable
     453                                CPPAD_ASSERT_UNKNOWN( play.GetOp( ind_taddr[j] ) == InvOp );
     454                                for_jac_sparsity.add_element( ind_taddr[j], i);
     455                        }
    382456                }
    383457        }
    384 
     458        else
     459        {       for(i = 0; i < n; i++)
     460                {       CPPAD_ASSERT_UNKNOWN( ind_taddr[i] < total_num_var );
     461                        // ind_taddr[i] is operator taddr for i-th independent variable
     462                        CPPAD_ASSERT_UNKNOWN( play.GetOp( ind_taddr[i] ) == InvOp );
     463
     464                        // add the elements that are present
     465                        itr = r[i].begin();
     466                        while( itr != r[i].end() )
     467                        {       j = *itr++;
     468                                CPPAD_ASSERT_KNOWN(
     469                                        j < q,
     470                                        "ForSparseJac: an element of the set r[i] "
     471                                        "has value greater than or equal q."
     472                                );
     473                                for_jac_sparsity.add_element( ind_taddr[i], j);
     474                        }
     475                }
     476        }
    385477        // evaluate the sparsity patterns
    386478        ForJacSweep(
     
    392484
    393485        // return values corresponding to dependent variables
    394         CPPAD_ASSERT_UNKNOWN( size_t(s.size()) == m );
     486        CPPAD_ASSERT_UNKNOWN( size_t(s.size()) == m || transpose );
     487        CPPAD_ASSERT_UNKNOWN( size_t(s.size()) == q || ! transpose );
    395488        for(i = 0; i < m; i++)
    396489        {       CPPAD_ASSERT_UNKNOWN( dep_taddr[i] < total_num_var );
    397                 CPPAD_ASSERT_UNKNOWN( s[i].empty() );
    398490
    399491                // extract results from for_jac_sparsity
     
    403495                j = for_jac_sparsity.next_element();
    404496                while( j < q )
    405                 {       s[i].insert(j);
     497                {       if( transpose )
     498                                s[j].insert(i);
     499                        else    s[i].insert(j);
    406500                        j = for_jac_sparsity.next_element();
    407501                }
     
    419513code depending on the value of \c VectorSet::value_type.
    420514
     515\param transpose
     516See \c ForSparseJac(q, r).
     517
    421518\param q
    422519See \c ForSparseJac(q, r).
     
    433530void ADFun<Base>::ForSparseJacCase(
    434531        bool                set_type      ,
     532        bool                transpose     ,
    435533        size_t              q             ,
    436534        const VectorSet&    r             ,
     
    444542        s.resize( m * q );
    445543
    446         CPPAD_ASSERT_KNOWN(
    447                 size_t(r.size()) == Domain() * q,
    448                 "ForSparseJac: using vectors of bools and r (second argument) length"
    449                 "\nis not equal to q (first argument) times domain dimension."
    450         );
    451 
    452         // store results in r and for_jac_sparse_pack_
     544        // store results in s and for_jac_sparse_pack_
    453545        ForSparseJacBool(
     546                transpose        ,
    454547                q                ,
    455548                r                ,
     
    475568code depending on the value of \c VectorSet::value_type.
    476569
     570\param transpose
     571See \c ForSparseJac(q, r).
     572
    477573\param q
    478574See \c ForSparseJac(q, r).
     
    488584void ADFun<Base>::ForSparseJacCase(
    489585        const std::set<size_t>&    set_type      ,
     586        bool                       transpose     ,
    490587        size_t                     q             ,
    491588        const VectorSet&           r             ,
    492589        VectorSet&                 s             )
    493 {
     590{       size_t m = Range();
    494591
    495592        // check VectorSet is Simple Vector class with sets for elements
     
    498595        );
    499596
    500         CPPAD_ASSERT_KNOWN(
    501                 size_t(r.size()) == Domain() ,
    502                 "ForSparseJac: using vectors of sets and r (second argument) length"
    503                 "\nis not equal to the domain dimension for ADFun object."
    504         );
    505 
    506597        // dimension size of result vector
    507         size_t m = Range();
    508         s.resize( m );
     598        if( transpose )
     599                s.resize(q);
     600        else    s.resize( m );
    509601
    510602        // store results in r and for_jac_sparse_pack_
    511603        CppAD::ForSparseJacSet(
     604                transpose        ,
    512605                q                ,
    513606                r                ,
     
    543636is a sparsity pattern for the matrix \f$ R \f$.
    544637
     638\param transpose
     639are sparsity patterns for \f$ R \f$ and \f$ S(x) \f$ transposed.
     640
    545641\return
    546 If \c VectorSet::value_type is \c bool,
    547 the return value \c s is a vector with size \c m*q
    548 where \c m is the number of dependent variables
    549 corresponding to the operation sequence stored in \c f.
    550 If \c VectorSet::value_type is \c std::set<size_t>,
    551 the return value \c s is a vector of sets with size \c m
    552 and with all its elements between zero and \a q - 1.
    553 The value of \a s is the sparsity pattern for the matrix
     642The value of \c transpose is false (true),
     643the return value is a sparsity pattern for \f$ S(x) \f$ (\f$ S(x)^T \f$) where
    554644\f[
    555645        S(x) = F^{(1)} (x) * R
     
    557647where \f$ F \f$ is the function corresponding to the operation sequence
    558648and \a x is any argument value.
     649If \c VectorSet::value_type is \c bool,
     650the return value has size \f$ m * q \f$ (\f$ q * m \f$).
     651where \c m is the number of dependent variables
     652corresponding to the operation sequence stored in \c f.
     653If \c VectorSet::value_type is \c std::set<size_t>,
     654the return value has size \f$ m \f$ ( \f$ q \f$ )
     655and with all its elements between zero and
     656\f$ q - 1 \f$ ( \f$ m - 1 \f$).
    559657
    560658\par Side Effects
     
    586684VectorSet ADFun<Base>::ForSparseJac(
    587685        size_t             q             ,
    588         const VectorSet&   r             )
     686        const VectorSet&   r             ,
     687        bool               transpose     )
    589688{       VectorSet s;
    590689        typedef typename VectorSet::value_type Set_type;
    591 
    592         CPPAD_ASSERT_KNOWN(
    593                 q > 0,
    594                 "ForSparseJac: q (first arugment) is not greater than zero"
    595         );
    596690
    597691        // free all memory currently in sparsity patterns
     
    601695        ForSparseJacCase(
    602696                Set_type()  ,
     697                transpose   ,
    603698                q           ,
    604699                r           ,
  • trunk/cppad/local/forward.hpp

    r2625 r2859  
    44
    55/* --------------------------------------------------------------------------
    6 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
     6CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-13 Bradley M. Bell
    77
    88CppAD is distributed under multiple licenses. This distribution is under
     
    2222        omh/forward.omh%
    2323        cppad/local/cap_taylor.hpp%
    24         example/forward.cpp
     24        example/forward.cpp%
     25        example/forward_mul.cpp
    2526%$$
    2627
     
    5152
    5253\param p
    53 is the order for this forward mode computation; i.e., the number
    54 of Taylor coefficient is <code>p+1</code>.
     54is the hightest order for this forward mode computation; i.e.,
     55after this calculation there will be <code>p+1</code>
     56Taylor coefficients per variables.
    5557
    5658\param x_p
    57 Is the \c p th order Taylor coefficient vector for the independent variables.
     59contains Taylor coefficients for the independent variables.
     60The size of \a x_p must either be \c n or <code>n*(p+1)</code>,
     61We define <code>q = p + 1 - x_p.size() / n</code>.
     62The Taylor coefficients of order k, for
     63k = q, ... , p are calculated.
    5864
    5965\param s
     
    6874        std::ostream& s             )
    6975{       // temporary indices
    70         size_t i, j;
     76        size_t i, j, k;
    7177
    7278        // number of independent variables
     
    8086
    8187        CPPAD_ASSERT_KNOWN(
    82                 size_t(x_p.size()) == n,
    83                 "Second argument to Forward does not have length equal to\n"
    84                 "the dimension of the domain for the corresponding ADFun."
     88                size_t(x_p.size()) == n || size_t(x_p.size()) == n*(p+1),
     89                "Forward: x_p.size() is not equal n or n*(p+1)"
    8590        );
     91        size_t n_order = size_t(x_p.size()) / n;
    8692        CPPAD_ASSERT_KNOWN(
    87                 p <= taylor_per_var_,
    88                 "The number of taylor_ coefficient currently stored\n"
    89                 "in this ADFun object is less than p."
     93                p <= taylor_per_var_ || n_order == p + 1,
     94                "The number of Taylor coefficient currently stored in this ADFun\n"
     95                "is less than p and x_p.size() != n*(p+1)."
    9096        ); 
    9197
     
    103109
    104110                // It is also variable taddr for j-th independent variable
    105                 taylor_[ind_taddr_[j] * taylor_col_dim_ + p] = x_p[j];
     111                if( n_order ==  1 )
     112                        taylor_[ind_taddr_[j] * taylor_col_dim_ + p] = x_p[j];
     113                else for(k = 0; k < n_order; k++)
     114                        taylor_[ind_taddr_[j] * taylor_col_dim_ + k] =
     115                                x_p[j * n_order + k];
    106116        }
    107117
    108118        // evaluate the derivatives
     119        size_t q = (p + 1) - n_order;
    109120        if( p == 0 )
    110121        {
     
    114125                );
    115126# else
    116                 compare_change_ = forward_sweep(s, true,
     127                compare_change_ = forward_sweep(s, true, q,
    117128                        p, n, total_num_var_, &play_, taylor_col_dim_, taylor_.data()
    118129                );
    119130# endif
    120131        }
    121         else forward_sweep(s, false,
    122                 p, n, total_num_var_, &play_, taylor_col_dim_, taylor_.data()
    123         );
     132        else if( q == 0 )
     133        {       compare_change_ = forward_sweep(s, true, q,
     134                        p, n, total_num_var_, &play_, taylor_col_dim_, taylor_.data()
     135                );
     136        }
     137        else
     138        {       forward_sweep(s, true, q,
     139                        p, n, total_num_var_, &play_, taylor_col_dim_, taylor_.data()
     140                );
     141        }
    124142
    125         // return the p-th order taylor_ coefficients for dependent variables
    126         Vector y_p(m);
    127         for(i = 0; i < m; i++)
    128         {       CPPAD_ASSERT_UNKNOWN( dep_taddr_[i] < total_num_var_ );
    129                 y_p[i] = taylor_[dep_taddr_[i] * taylor_col_dim_ + p];
     143        // return Taylor coefficients for dependent variables
     144        Vector y_p;
     145        if( n_order == 1 )
     146        {       y_p.resize(m);
     147                for(i = 0; i < m; i++)
     148                {       CPPAD_ASSERT_UNKNOWN( dep_taddr_[i] < total_num_var_ );
     149                        y_p[i] = taylor_[dep_taddr_[i] * taylor_col_dim_ + p];
     150                }
     151        }
     152        else
     153        {       y_p.resize(m * n_order );
     154                for(i = 0; i < m; i++) 
     155                {       for(k = 0; k < n_order; k++)
     156                                y_p[ i * n_order + k] =
     157                                        taylor_[ dep_taddr_[i] * taylor_col_dim_ + k ];
     158                }
    130159        }
    131160# ifndef NDEBUG
    132         if( hasnan(y_p) )
    133         {       if( p == 0 )
    134                 {       CPPAD_ASSERT_KNOWN(false,
    135                                 "y = f.Forward(0, x): has a nan in y."
    136                         ); 
    137                 }
    138                 else
    139                 {       CPPAD_ASSERT_KNOWN(false,
    140                                 "y_p = f.Forward(p, x_p): has a nan in y_p for p > 0, "
    141                                 "but not for p = 0."
    142                         );
     161        bool ok = true;
     162        if( p == 0 && n_order == 1 )
     163                ok = ! hasnan(y_p);
     164        else if( n_order != 1 )
     165        {       for(i = 0; i < m; i++)
     166                        ok &= ! isnan( y_p[ i * n_order + 0 ] );
     167        }
     168        CPPAD_ASSERT_KNOWN(ok,
     169                "y_p = f.Forward(p, x): has a zero order Taylor coefficient "
     170                "with the value nan."
     171        ); 
     172        if( p != 0 && n_order == 1 )
     173                ok = ! hasnan(y_p);
     174        else if( n_order != 1 )
     175        {       for(i = 0; i < m; i++)
     176                {       for(k = 1; k < n_order; k++)
     177                                ok &= ! isnan( y_p[ i * n_order + k ] );
    143178                }
    144179        }
     180        CPPAD_ASSERT_KNOWN(ok,
     181                "y_p = f.Forward(p, x): has a non-zero order Taylor coefficient\n"
     182                "with the value nan (but zero order coefficients are not nan)."
     183        );
    145184# endif
    146185
  • trunk/cppad/local/forward0sweep.hpp

    r2794 r2859  
    2121Compute zero order forward mode Taylor coefficients.
    2222*/
     23
     24/*
     25\def CPPAD_ATOMIC_CALL
     26This avoids warnings when NDEBUG is defined and user_ok is not used.
     27If \c NDEBUG is defined, this resolves to
     28\code
     29        user_atom->forward
     30\endcode
     31otherwise, it respolves to
     32\code
     33        user_ok = user_atom->forward
     34\endcode
     35This maco is undefined at the end of this file to facillitate is
     36use with a different definition in other files.
     37*/
     38# ifdef NDEBUG
     39# define CPPAD_ATOMIC_CALL user_atom->forward
     40# else
     41# define CPPAD_ATOMIC_CALL user_ok = user_atom->forward
     42# endif
    2343
    2444/*!
     
    141161
    142162        // work space used by UserOp.
    143         const size_t user_k = 0;     // order of this forward mode calculation
     163        const size_t user_q = 0;     // lowest order
     164        const size_t user_p = 0;     // highest order
     165        vector<bool> user_vx;        // empty vecotor
     166        vector<bool> user_vy;        // empty vecotor
    144167        vector<Base> user_tx;        // argument vector Taylor coefficients
    145168        vector<Base> user_ty;        // result vector Taylor coefficients
    146         size_t user_index = 0;       // indentifier for this user_atomic operation
     169        size_t user_index = 0;       // indentifier for this atomic operation
    147170        size_t user_id    = 0;       // user identifier for this call to operator
    148171        size_t user_i     = 0;       // index in result vector
     
    150173        size_t user_m     = 0;       // size of result vector
    151174        size_t user_n     = 0;       // size of arugment vector
     175        //
     176        atomic_base<Base>* user_atom = CPPAD_NULL; // user's atomic op calculator
     177# ifndef NDEBUG
     178        bool               user_ok   = false;      // atomic op return value
     179# endif
     180        //
    152181        // next expected operator in a UserOp sequence
    153182        enum { user_start, user_arg, user_ret, user_end } user_state = user_start;
     
    235264                        Rec->forward_csum(op, arg, i_op, i_var);
    236265                        forward_csum_op(
    237                                 0, i_var, arg, num_par, parameter, J, Taylor
     266                                0, 0, i_var, arg, num_par, parameter, J, Taylor
    238267                        );
    239268                        break;
     
    510539                                user_n     = arg[2];
    511540                                user_m     = arg[3];
    512                                 if(user_tx.size() < user_n)
     541                                user_atom  = atomic_base<Base>::list(user_index);
     542                                if(user_tx.size() != user_n)
    513543                                        user_tx.resize(user_n);
    514                                 if(user_ty.size() < user_m)
     544                                if(user_ty.size() != user_m)
    515545                                        user_ty.resize(user_m);
    516546                                user_j     = 0;
     
    524554                                CPPAD_ASSERT_UNKNOWN( user_n     == size_t(arg[2]) );
    525555                                CPPAD_ASSERT_UNKNOWN( user_m     == size_t(arg[3]) );
     556# ifndef NDEBUG
     557                                if( ! user_ok )
     558                                {       std::string msg = user_atom->afun_name()
     559                                        + ": atomic_base.forward: returned false";
     560                                        CPPAD_ASSERT_KNOWN(false, msg.c_str() );
     561                                }
     562# endif
    526563                                user_state = user_start;
    527564                        }
     
    536573                        if( user_j == user_n )
    537574                        {       // call users function for this operation
    538                                 user_atomic<Base>::forward(user_index, user_id,
    539                                         user_k, user_n, user_m, user_tx, user_ty
     575                                user_atom->set_id(user_id);
     576                                CPPAD_ATOMIC_CALL(user_q, user_p,
     577                                        user_vx, user_vy, user_tx, user_ty
    540578                                );
    541579                                user_state = user_ret;
     
    551589                        if( user_j == user_n )
    552590                        {       // call users function for this operation
    553                                 user_atomic<Base>::forward(user_index, user_id,
    554                                         user_k, user_n, user_m, user_tx, user_ty
     591                                user_atom->set_id(user_id);
     592                                CPPAD_ATOMIC_CALL(user_q, user_p,
     593                                        user_vx, user_vy, user_tx, user_ty
    555594                                );
    556595                                user_state = user_ret;
     
    612651// preprocessor symbols that are local to this file
    613652# undef CPPAD_FORWARD0SWEEP_TRACE
     653# undef CPPAD_ATOMIC_CALL
    614654
    615655# endif
  • trunk/cppad/local/forward_sweep.hpp

    r2794 r2859  
    2222*/
    2323
     24/*
     25\def CPPAD_ATOMIC_CALL
     26This avoids warnings when NDEBUG is defined and user_ok is not used.
     27If \c NDEBUG is defined, this resolves to
     28\code
     29        user_atom->forward
     30\endcode
     31otherwise, it respolves to
     32\code
     33        user_ok = user_atom->forward
     34\endcode
     35This maco is undefined at the end of this file to facillitate is
     36use with a different definition in other files.
     37*/
     38# ifdef NDEBUG
     39# define CPPAD_ATOMIC_CALL user_atom->forward
     40# else
     41# define CPPAD_ATOMIC_CALL user_ok = user_atom->forward
     42# endif
     43
    2444/*!
    2545\def CPPAD_FORWARD_SWEEP_TRACE
     
    4464\param print
    4565If \a print is false,
    46 suppress the output that is otherwise generated by the PriOp instructions
    47 (must be false when \c d is nonzero).
     66suppress the output that is otherwise generated by the PriOp instructions.
     67
     68\param q
     69is the lowest order of the Taylor coefficients
     70that are computed during this call.
    4871
    4972\param p
    50 is the order of the Taylor coefficients that are computed during this call.
     73is the highest order of the Taylor coefficients
     74that are computed during this call.
    5175
    5276\param n
     
    110134        std::ostream&         s_out,
    111135        const bool            print,
     136        const size_t          q,
    112137        const size_t          p,
    113138        const size_t          n,
    114139        const size_t          numvar,
    115140        player<Base>         *Rec,
    116         const size_t          J,
     141        const size_t         J,
    117142        Base                 *Taylor
    118143)
    119144{       CPPAD_ASSERT_UNKNOWN( J >= p + 1 );
     145        CPPAD_ASSERT_UNKNOWN( q <= p );
    120146
    121147        // op code for current instruction
     
    128154        size_t        i_var;
    129155
    130         CPPAD_ASSERT_UNKNOWN( p == 0 || ! print );
    131 # if CPPAD_USE_FORWARD0SWEEP
    132         CPPAD_ASSERT_UNKNOWN( p > 0 );
    133 # else
    134         addr_t*         non_const_arg;
    135 # endif
     156        // arg (not as a constant)
     157        addr_t*         non_const_arg = CPPAD_NULL;
     158
     159        // arg (as a constant)
    136160        const addr_t*   arg = CPPAD_NULL;
    137161
     
    155179        }
    156180
    157         // work space used by UserOp.
    158         const size_t user_k  = p;    // order of this forward mode calculation
    159         const size_t user_k1 = p+1;  // number of orders for this calculation
     181        // Work space used by UserOp. Note User assumes q = p.
     182        const size_t user_p1 = p+1;  // number of orders for this user calculation
     183        vector<bool> user_vx;        // empty vector
     184        vector<bool> user_vy;        // empty vector
    160185        vector<Base> user_tx;        // argument vector Taylor coefficients
    161186        vector<Base> user_ty;        // result vector Taylor coefficients
    162187        vector<size_t> user_iy;      // variable indices for results vector
    163         size_t user_index = 0;       // indentifier for this user_atomic operation
     188        size_t user_index = 0;       // indentifier for this atomic operation
    164189        size_t user_id    = 0;       // user identifier for this call to operator
    165190        size_t user_i     = 0;       // index in result vector
     
    167192        size_t user_m     = 0;       // size of result vector
    168193        size_t user_n     = 0;       // size of arugment vector
     194        //
     195        atomic_base<Base>* user_atom = CPPAD_NULL; // user's atomic op calculator
     196# ifndef NDEBUG
     197        bool               user_ok   = false;      // atomic op return value
     198# endif
     199        //
    169200        // next expected operator in a UserOp sequence
    170201        enum { user_start, user_arg, user_ret, user_end } user_state = user_start;
     
    181212                parameter = Rec->GetPar();
    182213
    183 # if ! CPPAD_USE_FORWARD0SWEEP
    184214        // length of the text vector (used by CppAD assert macros)
    185215        const size_t num_text = Rec->num_rec_text();
     
    189219        if( num_text > 0 )
    190220                text = Rec->GetTxt(0);
    191 # endif
    192221
    193222        // skip the BeginOp at the beginning of the recording
     
    208237                {
    209238                        case AbsOp:
    210                         forward_abs_op(p, i_var, arg[0], J, Taylor);
     239                        forward_abs_op(q, p, i_var, arg[0], J, Taylor);
    211240                        break;
    212241                        // -------------------------------------------------
    213242
    214243                        case AddvvOp:
    215                         forward_addvv_op(p, i_var, arg, parameter, J, Taylor);
     244                        forward_addvv_op(q, p, i_var, arg, parameter, J, Taylor);
    216245                        break;
    217246                        // -------------------------------------------------
     
    219248                        case AddpvOp:
    220249                        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_par );
    221                         forward_addpv_op(p, i_var, arg, parameter, J, Taylor);
     250                        forward_addpv_op(q, p, i_var, arg, parameter, J, Taylor);
    222251                        break;
    223252                        // -------------------------------------------------
     
    226255                        // sqrt(1 - x * x), acos(x)
    227256                        CPPAD_ASSERT_UNKNOWN( i_var < numvar  );
    228                         forward_acos_op(p, i_var, arg[0], J, Taylor);
     257                        forward_acos_op(q, p, i_var, arg[0], J, Taylor);
    229258                        break;
    230259                        // -------------------------------------------------
     
    233262                        // sqrt(1 - x * x), asin(x)
    234263                        CPPAD_ASSERT_UNKNOWN( i_var < numvar  );
    235                         forward_asin_op(p, i_var, arg[0], J, Taylor);
     264                        forward_asin_op(q, p, i_var, arg[0], J, Taylor);
    236265                        break;
    237266                        // -------------------------------------------------
     
    240269                        // 1 + x * x, atan(x)
    241270                        CPPAD_ASSERT_UNKNOWN( i_var < numvar  );
    242                         forward_atan_op(p, i_var, arg[0], J, Taylor);
     271                        forward_atan_op(q, p, i_var, arg[0], J, Taylor);
     272                        break;
     273                        // -------------------------------------------------
     274
     275                        case CExpOp:
     276                        forward_cond_op(
     277                                q, p, i_var, arg, num_par, parameter, J, Taylor
     278                        );
     279                        break;
     280                        // ---------------------------------------------------
     281
     282                        case ComOp:
     283                        if( q == 0 ) forward_comp_op_0(
     284                        compareCount, arg, num_par, parameter, J, Taylor
     285                        );
     286                        break;
     287                        // ---------------------------------------------------
     288
     289                        case CosOp:
     290                        // sin(x), cos(x)
     291                        CPPAD_ASSERT_UNKNOWN( i_var < numvar  );
     292                        forward_cos_op(q, p, i_var, arg[0], J, Taylor);
     293                        break;
     294                        // ---------------------------------------------------
     295
     296                        case CoshOp:
     297                        // sinh(x), cosh(x)
     298                        CPPAD_ASSERT_UNKNOWN( i_var < numvar  );
     299                        forward_cosh_op(q, p, i_var, arg[0], J, Taylor);
    243300                        break;
    244301                        // -------------------------------------------------
     
    250307                        Rec->forward_csum(op, arg, i_op, i_var);
    251308                        forward_csum_op(
    252                                 p, i_var, arg, num_par, parameter, J, Taylor
     309                                q, p, i_var, arg, num_par, parameter, J, Taylor
    253310                        );
    254311                        break;
    255312                        // -------------------------------------------------
    256313
    257                         case CExpOp:
    258                         forward_cond_op(
    259                                 p, i_var, arg, num_par, parameter, J, Taylor
    260                         );
    261                         break;
    262                         // ---------------------------------------------------
    263 
    264                         case ComOp:
    265 # if ! USE_FORWARD0SWEEP
    266                         if( p == 0 ) forward_comp_op_0(
    267                         compareCount, arg, num_par, parameter, J, Taylor
    268                         );
    269 # endif
    270                         break;
    271                         // ---------------------------------------------------
    272 
    273                         case CosOp:
    274                         // sin(x), cos(x)
    275                         CPPAD_ASSERT_UNKNOWN( i_var < numvar  );
    276                         forward_cos_op(p, i_var, arg[0], J, Taylor);
    277                         break;
    278                         // ---------------------------------------------------
    279 
    280                         case CoshOp:
    281                         // sinh(x), cosh(x)
    282                         CPPAD_ASSERT_UNKNOWN( i_var < numvar  );
    283                         forward_cosh_op(p, i_var, arg[0], J, Taylor);
    284                         break;
    285                         // -------------------------------------------------
    286 
    287314                        case DisOp:
    288 # if ! CPPAD_USE_FORWARD0SWEEP
    289                         if( p == 0 )
    290                                 forward_dis_op_0(i_var, arg, J, Taylor);
    291                         else
    292 # endif
    293                         {       Taylor[ i_var * J + p] = Base(0);
     315                        i = q;
     316                        if( i == 0 )
     317                        {       forward_dis_op_0(i_var, arg, J, Taylor);
     318                                i++;
     319                        }
     320                        while(i <= p)
     321                        {       Taylor[ i_var * J + i] = Base(0);
     322                                i++;
    294323                        }
    295324                        break;
     
    297326
    298327                        case DivvvOp:
    299                         forward_divvv_op(p, i_var, arg, parameter, J, Taylor);
     328                        forward_divvv_op(q, p, i_var, arg, parameter, J, Taylor);
    300329                        break;
    301330                        // -------------------------------------------------
     
    303332                        case DivpvOp:
    304333                        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_par );
    305                         forward_divpv_op(p, i_var, arg, parameter, J, Taylor);
     334                        forward_divpv_op(q, p, i_var, arg, parameter, J, Taylor);
    306335                        break;
    307336                        // -------------------------------------------------
     
    309338                        case DivvpOp:
    310339                        CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < num_par );
    311                         forward_divvp_op(p, i_var, arg, parameter, J, Taylor);
     340                        forward_divvp_op(q, p, i_var, arg, parameter, J, Taylor);
    312341                        break;
    313342                        // -------------------------------------------------
     
    319348
    320349                        case ExpOp:
    321                         forward_exp_op(p, i_var, arg[0], J, Taylor);
     350                        forward_exp_op(q, p, i_var, arg[0], J, Taylor);
    322351                        break;
    323352                        // -------------------------------------------------
     
    329358
    330359                        case LdpOp:
    331 # if ! CPPAD_USE_FORWARD0SWEEP
    332                         if( p == 0 )
     360                        if( q == 0 )
    333361                        {       non_const_arg = Rec->forward_non_const_arg();
    334362                                forward_load_p_op_0(
     
    343371                                        VectorInd.data()
    344372                                );
     373                                if( q < p )
     374                                        forward_load_op( op, q+1, p, i_var, arg, J, Taylor);
    345375                        }
    346376                        else
    347 # endif
    348                         {       forward_load_op( op, p, i_var, arg, J, Taylor);
     377                        {       forward_load_op( op, q, p, i_var, arg, J, Taylor);
    349378                        }
    350379                        break;
     
    352381
    353382                        case LdvOp:
    354 # if ! CPPAD_USE_FORWARD0SWEEP
    355                         if( p == 0 )
     383                        if( q == 0 )
    356384                        {       non_const_arg = Rec->forward_non_const_arg();
    357385                                forward_load_v_op_0(
     
    366394                                        VectorInd.data()
    367395                                );
     396                                if( q < p )
     397                                        forward_load_op( op, q+1, p, i_var, arg, J, Taylor);
    368398                        }
    369399                        else
    370 # endif
    371                         {       forward_load_op( op, p, i_var, arg, J, Taylor);
     400                        {       forward_load_op( op, q, p, i_var, arg, J, Taylor);
    372401                        }
    373402                        break;
     
    375404
    376405                        case LogOp:
    377                         forward_log_op(p, i_var, arg[0], J, Taylor);
     406                        forward_log_op(q, p, i_var, arg[0], J, Taylor);
    378407                        break;
    379408                        // -------------------------------------------------
    380409
    381410                        case MulvvOp:
    382                         forward_mulvv_op(p, i_var, arg, parameter, J, Taylor);
     411                        forward_mulvv_op(q, p, i_var, arg, parameter, J, Taylor);
    383412                        break;
    384413                        // -------------------------------------------------
     
    386415                        case MulpvOp:
    387416                        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_par );
    388                         forward_mulpv_op(p, i_var, arg, parameter, J, Taylor);
     417                        forward_mulpv_op(q, p, i_var, arg, parameter, J, Taylor);
    389418                        break;
    390419                        // -------------------------------------------------
    391420
    392421                        case ParOp:
    393 # if ! CPPAD_USE_FORWARD0SWEEP
    394                         if( p == 0 ) forward_par_op_0(
    395                                 i_var, arg, num_par, parameter, J, Taylor
    396                         );
    397                         else
    398 # endif
    399                         {       Taylor[ i_var * J + p] = Base(0);
     422                        i = q;
     423                        if( i == 0 )
     424                        {       forward_par_op_0(
     425                                        i_var, arg, num_par, parameter, J, Taylor
     426                                );
     427                                i++;
     428                        }
     429                        while(i <= p)
     430                        {       Taylor[ i_var * J + i] = Base(0);
     431                                i++;
    400432                        }
    401433                        break;
     
    404436                        case PowvpOp:
    405437                        CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < num_par );
    406                         forward_powvp_op(p, i_var, arg, parameter, J, Taylor);
     438                        forward_powvp_op(q, p, i_var, arg, parameter, J, Taylor);
    407439                        break;
    408440                        // -------------------------------------------------
     
    410442                        case PowpvOp:
    411443                        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_par );
    412                         forward_powpv_op(p, i_var, arg, parameter, J, Taylor);
     444                        forward_powpv_op(q, p, i_var, arg, parameter, J, Taylor);
    413445                        break;
    414446                        // -------------------------------------------------
    415447
    416448                        case PowvvOp:
    417                         forward_powvv_op(p, i_var, arg, parameter, J, Taylor);
     449                        forward_powvv_op(q, p, i_var, arg, parameter, J, Taylor);
    418450                        break;
    419451                        // -------------------------------------------------
    420452
    421453                        case PriOp:
    422 # if ! CPPAD_USE_FORWARD0SWEEP
    423                         if( print ) forward_pri_0(s_out,
     454                        if( (q == 0) & print ) forward_pri_0(s_out,
    424455                                i_var, arg, num_text, text, num_par, parameter, J, Taylor
    425456                        );
    426 # endif
    427457                        break;
    428458                        // -------------------------------------------------
    429459
    430460                        case SignOp:
    431                         // cos(x), sin(x)
    432                         CPPAD_ASSERT_UNKNOWN( i_var < numvar  );
    433                         forward_sign_op(p, i_var, arg[0], J, Taylor);
    434                         break;
     461                        // sign(x)
     462                        CPPAD_ASSERT_UNKNOWN( i_var < numvar  );
     463                        forward_sign_op(q, p, i_var, arg[0], J, Taylor);
     464                        break;
     465                        // -------------------------------------------------
    435466
    436467                        case SinOp:
    437468                        // cos(x), sin(x)
    438469                        CPPAD_ASSERT_UNKNOWN( i_var < numvar  );
    439                         forward_sin_op(p, i_var, arg[0], J, Taylor);
     470                        forward_sin_op(q, p, i_var, arg[0], J, Taylor);
    440471                        break;
    441472                        // -------------------------------------------------
     
    444475                        // cosh(x), sinh(x)
    445476                        CPPAD_ASSERT_UNKNOWN( i_var < numvar  );
    446                         forward_sinh_op(p, i_var, arg[0], J, Taylor);
     477                        forward_sinh_op(q, p, i_var, arg[0], J, Taylor);
    447478                        break;
    448479                        // -------------------------------------------------
    449480
    450481                        case SqrtOp:
    451                         forward_sqrt_op(p, i_var, arg[0], J, Taylor);
     482                        forward_sqrt_op(q, p, i_var, arg[0], J, Taylor);
    452483                        break;
    453484                        // -------------------------------------------------
    454485
    455486                        case StppOp:
    456 # if ! CPPAD_USE_FORWARD0SWEEP
    457                         if( p == 0 )
     487                        if( q == 0 )
    458488                        {       forward_store_pp_op_0(
    459489                                        i_var,
     
    467497                                );
    468498                        }
    469 # endif
    470499                        break;
    471500                        // -------------------------------------------------
    472501
    473502                        case StpvOp:
    474 # if ! CPPAD_USE_FORWARD0SWEEP
    475                         if( p == 0 )
     503                        if( q == 0 )
    476504                        {       forward_store_pv_op_0(
    477505                                        i_var,
     
    485513                                );
    486514                        }
    487 # endif
    488515                        break;
    489516                        // -------------------------------------------------
    490517
    491518                        case StvpOp:
    492 # if ! CPPAD_USE_FORWARD0SWEEP
    493                         if( p == 0 )
     519                        if( q == 0 )
    494520                        {       forward_store_vp_op_0(
    495521                                        i_var,
     
    503529                                );
    504530                        }
    505 # endif
    506531                        break;
    507532                        // -------------------------------------------------
    508533
    509534                        case StvvOp:
    510 # if ! CPPAD_USE_FORWARD0SWEEP
    511                         if( p == 0 )
     535                        if( q == 0 )
    512536                        {       forward_store_vv_op_0(
    513537                                        i_var,
     
    521545                                );
    522546                        }
    523 # endif
    524547                        break;
    525548                        // -------------------------------------------------
    526549
    527550                        case SubvvOp:
    528                         forward_subvv_op(p, i_var, arg, parameter, J, Taylor);
     551                        forward_subvv_op(q, p, i_var, arg, parameter, J, Taylor);
    529552                        break;
    530553                        // -------------------------------------------------
     
    532555                        case SubpvOp:
    533556                        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_par );
    534                         forward_subpv_op(p, i_var, arg, parameter, J, Taylor);
     557                        forward_subpv_op(q, p, i_var, arg, parameter, J, Taylor);
    535558                        break;
    536559                        // -------------------------------------------------
     
    538561                        case SubvpOp:
    539562                        CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < num_par );
    540                         forward_subvp_op(p, i_var, arg, parameter, J, Taylor);
     563                        forward_subvp_op(q, p, i_var, arg, parameter, J, Taylor);
    541564                        break;
    542565                        // -------------------------------------------------
     
    545568                        // tan(x)^2, tan(x)
    546569                        CPPAD_ASSERT_UNKNOWN( i_var < numvar  );
    547                         forward_tan_op(p, i_var, arg[0], J, Taylor);
     570                        forward_tan_op(q, p, i_var, arg[0], J, Taylor);
    548571                        break;
    549572                        // -------------------------------------------------
     
    552575                        // tanh(x)^2, tanh(x)
    553576                        CPPAD_ASSERT_UNKNOWN( i_var < numvar  );
    554                         forward_tanh_op(p, i_var, arg[0], J, Taylor);
     577                        forward_tanh_op(q, p, i_var, arg[0], J, Taylor);
    555578                        break;
    556579                        // -------------------------------------------------
     
    565588                                user_n     = arg[2];
    566589                                user_m     = arg[3];
    567                                 if(user_tx.size() < user_n * user_k1)
    568                                         user_tx.resize(user_n * user_k1);
    569                                 if(user_ty.size() < user_m * user_k1)
    570                                         user_ty.resize(user_m * user_k1);
    571                                 if(user_iy.size() < user_m)
     590                                user_atom  = atomic_base<Base>::list(user_index);
     591                                if(user_tx.size() != user_n * user_p1)
     592                                        user_tx.resize(user_n * user_p1);
     593                                if(user_ty.size() != user_m * user_p1)
     594                                        user_ty.resize(user_m * user_p1);
     595                                if(user_iy.size() != user_m)
    572596                                        user_iy.resize(user_m);
    573597                                user_j     = 0;
     
    581605                                CPPAD_ASSERT_UNKNOWN( user_n     == size_t(arg[2]) );
    582606                                CPPAD_ASSERT_UNKNOWN( user_m     == size_t(arg[3]) );
     607
     608                                // call users function for this operation
     609                                user_atom->set_id(user_id);
     610                                CPPAD_ATOMIC_CALL(
     611                                        q, p, user_vx, user_vy, user_tx, user_ty
     612                                );
     613# ifndef NDEBUG
     614                                if( ! user_ok )
     615                                {       std::string msg = user_atom->afun_name()
     616                                        + ": atomic_base.forward: returned false";
     617                                        CPPAD_ASSERT_KNOWN(false, msg.c_str() );
     618                                }
     619# endif
     620                                for(i = 0; i < user_m; i++)
     621                                        if( user_iy[i] > 0 )
     622                                                for(ell = q; ell <= p; ell++)
     623                                                        Taylor[ user_iy[i] * J + ell ] =
     624                                                                user_ty[ i * user_p1 + ell ];
     625
    583626                                user_state = user_start;
    584 
    585                                 // call users function for this operation
    586                                 user_atomic<Base>::forward(user_index, user_id,
    587                                         user_k, user_n, user_m, user_tx, user_ty
    588                                 );
    589                                 for(i = 0; i < user_m; i++) if( user_iy[i] > 0 )
    590                                         Taylor[ user_iy[i] * J + user_k ] =
    591                                                 user_ty[ i * user_k1 + user_k ];
    592627                        }
    593628                        break;
     
    598633                        CPPAD_ASSERT_UNKNOWN( user_j < user_n );
    599634                        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_par );
    600                         user_tx[user_j * user_k1 + 0] = parameter[ arg[0]];
    601                         for(ell = 1; ell < user_k1; ell++)
    602                                 user_tx[user_j * user_k1 + ell] = Base(0);
     635                        user_tx[user_j * user_p1 + 0] = parameter[ arg[0]];
     636                        for(ell = 1; ell < user_p1; ell++)
     637                                user_tx[user_j * user_p1 + ell] = Base(0);
    603638                        ++user_j;
    604639                        if( user_j == user_n )
     
    611646                        CPPAD_ASSERT_UNKNOWN( user_j < user_n );
    612647                        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) <= i_var );
    613                         for(ell = 0; ell < user_k1; ell++)
    614                                 user_tx[user_j * user_k1 + ell] = Taylor[ arg[0] * J + ell];
     648                        for(ell = 0; ell < user_p1; ell++)
     649                                user_tx[user_j * user_p1 + ell] = Taylor[ arg[0] * J + ell];
    615650                        ++user_j;
    616651                        if( user_j == user_n )
     
    623658                        CPPAD_ASSERT_UNKNOWN( user_i < user_m );
    624659                        user_iy[user_i] = 0;
    625                         user_ty[user_i * user_k1 + 0] = parameter[ arg[0]];
    626                         for(ell = 1; ell < user_k; ell++)
    627                                 user_ty[user_i * user_k1 + ell] = Base(0);
     660                        user_ty[user_i * user_p1 + 0] = parameter[ arg[0]];
     661                        for(ell = 1; ell < q; ell++)
     662                                user_ty[user_i * user_p1 + ell] = Base(0);
    628663                        user_i++;
    629664                        if( user_i == user_m )
     
    636671                        CPPAD_ASSERT_UNKNOWN( user_i < user_m );
    637672                        user_iy[user_i] = i_var;
    638                         for(ell = 0; ell < user_k; ell++)
    639                                 user_ty[user_i * user_k1 + ell] = Taylor[ i_var * J + ell];
     673                        for(ell = 0; ell < q; ell++)
     674                                user_ty[user_i * user_p1 + ell] = Taylor[ i_var * J + ell];
    640675                        user_i++;
    641676                        if( user_i == user_m )
     
    672707}
    673708
     709// preprocessor symbols that are local to this file
    674710# undef CPPAD_FORWARD_SWEEP_TRACE
     711# undef CPPAD_ATOMIC_CALL
    675712
    676713/*! \} */
  • trunk/cppad/local/fun_construct.hpp

    r2625 r2859  
    44
    55/* --------------------------------------------------------------------------
    6 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
     6CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-13 Bradley M. Bell
    77
    88CppAD is distributed under multiple licenses. This distribution is under
     
    454454        size_t p = 0;
    455455        compare_change_ = forward_sweep(std::cout, false,
    456                 p, n, total_num_var_, &play_, taylor_col_dim_, taylor_.data()
     456                p, p, n, total_num_var_, &play_, taylor_col_dim_, taylor_.data()
    457457        );
    458458# endif
  • trunk/cppad/local/load_op.hpp

    r2625 r2859  
    44
    55/* --------------------------------------------------------------------------
    6 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
     6CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-13 Bradley M. Bell
    77
    88CppAD is distributed under multiple licenses. This distribution is under
     
    138138(only used for error checking).
    139139
    140 \param d
    141 is the order of the Taylor coefficient that we are computing.
     140\param q
     141is the lowest order of the Taylor coefficient that we are computing.
     142
     143\param p
     144is the highest order of the Taylor coefficient that we are computing.
    142145
    143146\param i_z
     
    146149\param arg
    147150\a arg[2]
    148 If y[x] is a parameter, \a arg[2] is zero
     151If y[x] is a parameter, <code>arg[2]</code> is zero
    149152(which is not a valid variable index).
    150153If y[x] is a variable,
    151 \a arg[2] is the variable index corresponding to y[x].
     154<code>arg[2]</code> is the variable index corresponding to y[x].
    152155
    153156\param nc_taylor
     
    155158
    156159\param taylor
    157 \b Input: if y[x] is a variable, \a taylor[ \a arg[2] * nc_taylor + d ]
    158 is the d-order Taylor coefficient corresponding to y[x].
     160\b Input: if y[x] is a variable,
     161<code>taylor[ arg[2] * nc_taylor + k ]</code>
     162for k = 0 , ... , p,
     163is the k-order Taylor coefficient corresponding to y[x].
    159164\n
    160 \b Output: \a taylor[ i_z * nc_taylor + d ]
    161 is the d-order Taylor coefficient for the variable z.
     165\b Output: <code>taylor[ i_z * nc_taylor + d ]</code>
     166for k = q , ... , p,
     167is the k-order Taylor coefficient for the variable z.
    162168
    163169\par Checked Assertions
    164170\li NumArg(op) == 3
    165171\li NumRes(op) == 1
    166 \li 0 < d < nc_taylor
     172\li p < nc_taylor
     173\li 0 < q <= p
    167174\li size_t(arg[2]) < i_z
    168175*/
     
    170177inline void forward_load_op(
    171178        OpCode         op          ,
    172         size_t         d           ,
     179        size_t         q           ,
     180        size_t         p           ,
    173181        size_t         i_z         ,
    174182        const addr_t*  arg         ,
     
    179187        CPPAD_ASSERT_UNKNOWN( NumArg(op) == 3 );
    180188        CPPAD_ASSERT_UNKNOWN( NumRes(op) == 1 );
    181         CPPAD_ASSERT_UNKNOWN( d > 0 )
    182         CPPAD_ASSERT_UNKNOWN( d < nc_taylor );
     189        CPPAD_ASSERT_UNKNOWN( p < nc_taylor );
     190        CPPAD_ASSERT_UNKNOWN( 0 < q && q <= p );
    183191        CPPAD_ASSERT_UNKNOWN( size_t(arg[2]) < i_z );
    184192
     
    186194        if( arg[2] > 0 )
    187195        {       Base* y_x = taylor + arg[2] * nc_taylor;
    188                 z[d]      = y_x[d];
    189         }
    190         else    z[d]      = Base(0);
     196                for(size_t d = q; d <= p; d++)
     197                        z[d] = y_x[d];
     198        }
     199        else
     200        {       for(size_t d = q; d <= p; d++)
     201                        z[d] = Base(0);
     202        }
    191203}
    192204
  • trunk/cppad/local/log_op.hpp

    r2625 r2859  
    44
    55/* --------------------------------------------------------------------------
    6 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
     6CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-13 Bradley M. Bell
    77
    88CppAD is distributed under multiple licenses. This distribution is under
     
    3434template <class Base>
    3535inline void forward_log_op(
    36         size_t j           ,
     36        size_t q           ,
     37        size_t p           ,
    3738        size_t i_z         ,
    3839        size_t i_x         ,
     
    4647        CPPAD_ASSERT_UNKNOWN( NumRes(LogOp) == 1 );
    4748        CPPAD_ASSERT_UNKNOWN( i_x < i_z );
    48         CPPAD_ASSERT_UNKNOWN( j < nc_taylor );
     49        CPPAD_ASSERT_UNKNOWN( p < nc_taylor );
     50        CPPAD_ASSERT_UNKNOWN( q <= p );
    4951
    5052        // Taylor coefficients corresponding to argument and result
     
    5254        Base* z = taylor + i_z * nc_taylor;
    5355
    54         if( j == 0 )
    55                 z[0] = log( x[0] );
    56         else if ( j == 1 )
    57                 z[1] = x[1] / x[0];
    58         else
     56        if( q == 0 )
     57        {       z[0] = log( x[0] );
     58                q++;
     59                if( p == 0 )
     60                        return;
     61        }
     62        if ( q == 1 )
     63        {       z[1] = x[1] / x[0];
     64                q++;
     65        }
     66        for(size_t j = q; j <= p; j++)
    5967        {
    6068                z[j] = -z[1] * x[j-1];
  • trunk/cppad/local/mul_op.hpp

    r2625 r2859  
    44
    55/* --------------------------------------------------------------------------
    6 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
     6CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-13 Bradley M. Bell
    77
    88CppAD is distributed under multiple licenses. This distribution is under
     
    3939template <class Base>
    4040inline void forward_mulvv_op(
    41         size_t        d           ,
     41        size_t        q           ,
     42        size_t        p           ,
    4243        size_t        i_z         ,
    4344        const addr_t* arg         ,
     
    5152        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < i_z );
    5253        CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < i_z );
    53         CPPAD_ASSERT_UNKNOWN( d < nc_taylor );
     54        CPPAD_ASSERT_UNKNOWN( p < nc_taylor );
     55        CPPAD_ASSERT_UNKNOWN( q <= p );
    5456
    5557        // Taylor coefficients corresponding to arguments and result
     
    5961
    6062        size_t k;
    61         z[d] = Base(0);
    62         for(k = 0; k <= d; k++)
    63                 z[d] += x[d-k] * y[k];
     63        for(size_t d = q; d <= p; d++)
     64        {       z[d] = Base(0);
     65                for(k = 0; k <= d; k++)
     66                        z[d] += x[d-k] * y[k];
     67        }
    6468}
    6569
     
    171175template <class Base>
    172176inline void forward_mulpv_op(
    173         size_t        d           ,
     177        size_t        q           ,
     178        size_t        p           ,
    174179        size_t        i_z         ,
    175180        const addr_t* arg         ,
     
    182187        CPPAD_ASSERT_UNKNOWN( NumRes(MulpvOp) == 1 );
    183188        CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < i_z );
    184         CPPAD_ASSERT_UNKNOWN( d < nc_taylor );
     189        CPPAD_ASSERT_UNKNOWN( p < nc_taylor );
     190        CPPAD_ASSERT_UNKNOWN( q <= p );
    185191
    186192        // Taylor coefficients corresponding to arguments and result
     
    191197        Base x = parameter[ arg[0] ];
    192198
    193         z[d] = x * y[d];
     199        for(size_t d = q; d <= p; d++)
     200                z[d] = x * y[d];
    194201}
    195202/*!
  • trunk/cppad/local/op_code.hpp

    r2625 r2859  
    44
    55/* --------------------------------------------------------------------------
    6 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
     6CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-13 Bradley M. Bell
    77
    88CppAD is distributed under multiple licenses. This distribution is under
     
    495495                CPPAD_ASSERT_UNKNOWN( ind[3+ind[0]+ind[1]] == ind[0]+ind[1] );
    496496                printOpField(os, " pr=", Rec->GetPar(ind[2]), ncol);
    497                 for(i = 0; i < ind[0]; i++)
     497                for(i = 0; i < size_t(ind[0]); i++)
    498498                         printOpField(os, " +v=", ind[3+i], ncol);
    499                 for(i = 0; i < ind[1]; i++)
     499                for(i = 0; i < size_t(ind[1]); i++)
    500500                         printOpField(os, " -v=", ind[3+ind[0]+i], ncol);
    501501                break;
     
    597597                case UserOp:
    598598                CPPAD_ASSERT_UNKNOWN( NumArg(op) == 4 );
    599                 {       const char* name = user_atomic<Base>::name(ind[0]);
    600                         printOpField(os, " f=",   name, ncol);
     599                {       atomic_base<Base>* atom = atomic_base<Base>::list(ind[0]);
     600                        std::string name = atom->name();
     601                        printOpField(os, " f=",   name.c_str(), ncol);
    601602                        printOpField(os, " i=", ind[1], ncol);
    602603                        printOpField(os, " n=", ind[2], ncol);
  • trunk/cppad/local/optimize.hpp

    r2683 r2859  
    44
    55/* --------------------------------------------------------------------------
    6 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
     6CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-13 Bradley M. Bell
    77
    88CppAD is distributed under multiple licenses. This distribution is under
     
    4646$head f$$
    4747The object $icode f$$ has prototype
    48 $icode%
     48$codei%
    4949        ADFun<%Base%> %f%
    5050%$$
     
    10711071                                user_m     = arg[3];
    10721072                                user_q     = 1;
    1073                                 if(user_r.size() < user_n )
    1074                                         user_r.resize(user_n);
    1075                                 if(user_s.size() < user_m )
    1076                                         user_s.resize(user_m);
     1073                                if(user_s.size() != user_n )
     1074                                        user_s.resize(user_n);
     1075                                if(user_r.size() != user_m )
     1076                                        user_r.resize(user_m);
    10771077                                user_j     = user_n;
    10781078                                user_i     = user_m;
     
    11091109                        CPPAD_ASSERT_UNKNOWN( 0 < arg[0] );
    11101110                        --user_j;
    1111                         if( ! user_r[user_j].empty() )
     1111                        if( ! user_s[user_j].empty() )
    11121112                        {       tape[arg[0]].connect = yes_connected;
    11131113                                user_keep.top() = true;
     
    11241124                        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_par );
    11251125                        --user_i;
    1126                         user_s[user_i].clear();
     1126                        user_r[user_i].clear();
    11271127                        if( user_i == 0 )
    11281128                        {       // call users function for this operation
    1129                                 user_atomic<Base>::rev_jac_sparse(user_index, user_id,
    1130                                         user_n, user_m, user_q, user_r, user_s
     1129                                atomic_base<Base>* atom =
     1130                                        atomic_base<Base>::list(user_index);
     1131                                atom->set_id(user_id);
     1132                                atom->rev_sparse_jac(
     1133                                        user_q, user_r, user_s
    11311134                                );
    11321135                                user_state = user_arg;
     
    11391142                        CPPAD_ASSERT_UNKNOWN( 0 < user_i && user_i <= user_m );
    11401143                        --user_i;
    1141                         user_s[user_i].clear();
     1144                        user_r[user_i].clear();
    11421145                        if( tape[i_var].connect != not_connected )
    1143                                 user_s[user_i].insert(0);
     1146                                user_r[user_i].insert(0);
    11441147                        if( user_i == 0 )
    11451148                        {       // call users function for this operation
    1146                                 user_atomic<Base>::rev_jac_sparse(user_index, user_id,
    1147                                         user_n, user_m, user_q, user_r, user_s
     1149                                atomic_base<Base>* atom =
     1150                                        atomic_base<Base>::list(user_index);
     1151                                atom->set_id(user_id);
     1152                                atom->rev_sparse_jac(
     1153                                        user_q, user_r, user_s
    11481154                                );
    11491155                                user_state = user_arg;
  • trunk/cppad/local/pow_op.hpp

    r2625 r2859  
    44
    55/* --------------------------------------------------------------------------
    6 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
     6CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-13 Bradley M. Bell
    77
    88CppAD is distributed under multiple licenses. This distribution is under
     
    3535template <class Base>
    3636inline void forward_powvv_op(
    37         size_t        d           ,
     37        size_t        q           ,
     38        size_t        p           ,
    3839        size_t        i_z         ,
    3940        const addr_t* arg         ,
     
    4344{
    4445        // convert from final result to first result
    45         i_z -= 2; // NumRes(PowvvOp) - 1;
     46        i_z -= 2; // 2 = NumRes(PowvvOp) - 1;
    4647
    4748        // check assumptions
     
    5051        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < i_z );
    5152        CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < i_z );
    52         CPPAD_ASSERT_UNKNOWN( d < nc_taylor );
    53 
    54         // z_0 = log(x)
    55         forward_log_op(d, i_z, arg[0], nc_taylor, taylor);
     53        CPPAD_ASSERT_UNKNOWN( p < nc_taylor );
     54        CPPAD_ASSERT_UNKNOWN( q <= p );
     55
     56        // z_0 = log(x)
     57        forward_log_op(q, p, i_z, arg[0], nc_taylor, taylor);
    5658
    5759        // z_1 = z_0 * y
     
    5961        adr[0] = i_z;
    6062        adr[1] = arg[1];
    61         forward_mulvv_op(d, i_z+1, adr, parameter, nc_taylor, taylor);
    62 
    63         // z_2 = exp(z_1)
    64 # if CPPAD_USE_FORWARD0SWEEP
    65         CPPAD_ASSERT_UNKNOWN( d > 0 );
    66         forward_exp_op(d, i_z+2, i_z+1, nc_taylor, taylor);
    67 # else
     63        forward_mulvv_op(q, p, i_z+1, adr, parameter, nc_taylor, taylor);
     64
     65        // z_2 = exp(z_1)
    6866        // final result for zero order case is exactly the same as for Base
    69         if( d == 0 )
     67        if( q == 0 )
    7068        {       // Taylor coefficients corresponding to arguments and result
    7169                Base* x   = taylor + arg[0]  * nc_taylor;
     
    7472
    7573                z_2[0] = pow(x[0], y[0]);
     74                q++;
    7675        }
    77         else    forward_exp_op(d, i_z+2, i_z+1, nc_taylor, taylor);
    78 # endif
     76        if( q <= p )
     77                forward_exp_op(q, p, i_z+2, i_z+1, nc_taylor, taylor);
    7978}
    8079
     
    194193template <class Base>
    195194inline void forward_powpv_op(
    196         size_t        d           ,
     195        size_t        q           ,
     196        size_t        p           ,
    197197        size_t        i_z         ,
    198198        const addr_t* arg         ,
     
    202202{
    203203        // convert from final result to first result
    204         i_z -= 2; // NumRes(PowpvOp) - 1;
     204        i_z -= 2; // 2 = NumRes(PowpvOp) - 1;
    205205
    206206        // check assumptions
     
    208208        CPPAD_ASSERT_UNKNOWN( NumRes(PowpvOp) == 3 );
    209209        CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < i_z );
    210         CPPAD_ASSERT_UNKNOWN( d < nc_taylor );
     210        CPPAD_ASSERT_UNKNOWN( p < nc_taylor );
     211        CPPAD_ASSERT_UNKNOWN( q <= p );
    211212
    212213        // Taylor coefficients corresponding to arguments and result
     
    214215
    215216        // z_0 = log(x)
    216 # if CPPAD_USE_FORWARD0SWEEP
    217         CPPAD_ASSERT_UNKNOWN( d > 0 );
    218         z_0[d] = Base(0);
    219 # else
    220217        Base x    = parameter[ arg[0] ];
    221         if( d == 0 )
    222                 z_0[0] = log(x);
    223         else    z_0[d] = Base(0);
    224 # endif
     218        size_t d;
     219        for(d = q; d <= p; d++)
     220        {       if( d == 0 )
     221                        z_0[d] = log(x);
     222                else    z_0[d] = Base(0);
     223        }
     224
    225225        // z_1 = z_0 * y
    226226        addr_t adr[2];
    227         adr[0] = i_z * nc_taylor; // offset of z_0[0] in taylor
     227        adr[0] = i_z * nc_taylor; // offset of z_0 in taylor; i.e., log(x)
    228228        adr[1] = arg[1];          // variable index of y in taylor
    229         // use taylor both for parameter and variable values
    230         forward_mulpv_op(d, i_z+1, adr, taylor, nc_taylor, taylor);
    231 
    232         // z_2 = exp(z_1)
    233 # if CPPAD_USE_FORWARD0SWEEP
    234         forward_exp_op(d, i_z+2, i_z+1, nc_taylor, taylor);
    235 # else
     229
     230        // use taylor both for parameter and variable values (trick)
     231        forward_mulpv_op(q, p, i_z+1, adr, taylor, nc_taylor, taylor);
     232
     233        // z_2 = exp(z_1)
    236234        // zero order case exactly same as Base type operation
    237         if( d == 0 )
     235        if( q == 0 )
    238236        {       Base* y   = taylor + arg[1]  * nc_taylor;
    239237                Base* z_2 = taylor + (i_z+2) * nc_taylor;
    240238                z_2[0] = pow(x, y[0]);
     239                q++;
    241240        }
    242         else    forward_exp_op(d, i_z+2, i_z+1, nc_taylor, taylor);
    243 # endif
     241        if( q <= p )
     242                forward_exp_op(q, p, i_z+2, i_z+1, nc_taylor, taylor);
    244243}
    245244/*!
     
    360359template <class Base>
    361360inline void forward_powvp_op(
    362         size_t        d           ,
     361        size_t        q           ,
     362        size_t        p           ,
    363363        size_t        i_z         ,
    364364        const addr_t* arg         ,
     
    368368{
    369369        // convert from final result to first result
    370         i_z -= 2; // NumRes(PowvpOp) - 1;
     370        i_z -= 2; // 2 = NumRes(PowvpOp) - 1
    371371
    372372        // check assumptions
     
    374374        CPPAD_ASSERT_UNKNOWN( NumRes(PowvpOp) == 3 );
    375375        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < i_z );
    376         CPPAD_ASSERT_UNKNOWN( d < nc_taylor );
    377 
    378         // z_0 = log(x)
    379         forward_log_op(d, i_z, arg[0], nc_taylor, taylor);
     376        CPPAD_ASSERT_UNKNOWN( p < nc_taylor );
     377        CPPAD_ASSERT_UNKNOWN( q <= p );
     378
     379        // z_0 = log(x)
     380        forward_log_op(q, p, i_z, arg[0], nc_taylor, taylor);
    380381
    381382        // z_1 = y * z_0
     
    383384        adr[0] = arg[1];
    384385        adr[1] = i_z;
    385         forward_mulpv_op(d, i_z+1, adr, parameter, nc_taylor, taylor);
     386        forward_mulpv_op(q, p, i_z+1, adr, parameter, nc_taylor, taylor);
    386387
    387388        // z_2 = exp(z_1)
    388389        // zero order case exactly same as Base type operation
    389 # if CPPAD_USE_FORWARD0SWEEP
    390         CPPAD_ASSERT_UNKNOWN( d > 0 );
    391         forward_exp_op(d, i_z+2, i_z+1, nc_taylor, taylor);
    392 # else
    393         if( d == 0 )
     390        if( q == 0 )
    394391        {       Base* z_2 = taylor + (i_z+2) * nc_taylor;
    395392                Base* x   = taylor + arg[0] * nc_taylor;
    396393                Base  y   = parameter[ arg[1] ];
    397394                z_2[0]  = pow(x[0], y);
     395                q++;
    398396        }
    399         else    forward_exp_op(d, i_z+2, i_z+1, nc_taylor, taylor);
    400 # endif
    401 
     397        if( q <= p )
     398                forward_exp_op(q, p, i_z+2, i_z+1, nc_taylor, taylor);
    402399}
    403400
  • trunk/cppad/local/print_for.hpp

    r2691 r2859  
    44
    55/* --------------------------------------------------------------------------
    6 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
     6CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-13 Bradley M. Bell
    77
    88CppAD is distributed under multiple licenses. This distribution is under
     
    4444$head Purpose$$
    4545The $cref/zero order forward/ForwardZero/$$ mode command
    46 $icode%
    47         f%.Forward(0, %x%)
     46$codei%
     47        %f%.Forward(0, %x%)
    4848%$$
    4949assigns the
  • trunk/cppad/local/prototype_op.hpp

    r2625 r2859  
    44
    55/* --------------------------------------------------------------------------
    6 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
     6CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-13 Bradley M. Bell
    77
    88CppAD is distributed under multiple licenses. This distribution is under
     
    3434\a Base.
    3535
    36 \param j
    37 order of the Taylor coefficient that we are computing.
     36\param q
     37lowest order of the Taylor coefficient that we are computing.
     38
     39\param p
     40highest order of the Taylor coefficient that we are computing.
    3841
    3942\param i_z
     
    4952
    5053\param taylor
    51 \b Input: \a taylor [ \a i_x * \a nc_taylor + k ]
    52 for k = 0 , ... , \a j
     54\b Input: <code>taylor [ i_x * nc_taylor + k ]</code>,
     55for k = 0 , ... , p,
    5356is the k-th order Taylor coefficient corresponding to x.
    5457\n
    55 \b Input: \a taylor [ \a i_z * \a nc_taylor + k ]
    56 for k = 0 , ... , \a j - 1
     58\b Input: <code>taylor [ i_z * nc_taylor + k ]</code>,
     59for k = 0 , ... , q-1,
    5760is the k-th order Taylor coefficient corresponding to z.
    5861\n
    59 \b Output: \a taylor [ \a i_z * \a nc_taylor + j ]
    60 is the j-th order Taylor coefficient corresponding to z.
     62\b Output: <code>taylor [ i_z * nc_taylor + k ]</code>,
     63for k = q , ... , p,
     64is the k-th order Taylor coefficient corresponding to z.
    6165
    6266\par Checked Assertions
    6367\li NumArg(op) == 1
    6468\li NumRes(op) == 1
    65 \li \a i_x < \a i_z
    66 \li \a j < \a nc_taylor
     69\li i_x < i_z
     70\li p < nc_taylor
     71\li q <= p
    6772*/
    6873template <class Base>
    6974inline void forward_unary1_op(
    70         size_t j           ,
     75        size_t q           ,
     76        size_t p           ,
    7177        size_t i_z         ,
    7278        size_t i_x         ,
     
    214220\a Base.
    215221
    216 \param j
    217 order of the Taylor coefficients that we are computing.
     222\param q
     223lowest order of the Taylor coefficients that we are computing.
     224
     225\param p
     226highest order of the Taylor coefficients that we are computing.
    218227
    219228\param i_z
     
    230239
    231240\param taylor
    232 \b Input: \a taylor [ \a i_x * \a nc_taylor + k ]
    233 for k = 0 , ... , \a j
     241\b Input: <code>taylor [ i_x * nc_taylor + k ]</code>
     242for k = 0 , ... , p,
    234243is the k-th order Taylor coefficient corresponding to x.
    235244\n
    236 \b Input: \a taylor [ \a i_z * \a nc_taylor + k ]
    237 for k = 0 , ... , \a j - 1
     245\b Input: <code>taylor [ i_z * nc_taylor + k ]</code>
     246for k = 0 , ... , q - 1,
    238247is the k-th order Taylor coefficient corresponding to z.
    239248\n
    240 \b Input: \a taylor [ ( \a i_z - 1) * \a nc_taylor + k ]
    241 for k = 0 , ... , \a j - 1
     249\b Input: <code>taylor [ ( i_z - 1) * nc_taylor + k ]</code>
     250for k = 0 , ... , q-1,
    242251is the k-th order Taylor coefficient corresponding to the auxillary result y.
    243252\n
    244 \b Output: \a taylor [ \a i_z * \a nc_taylor + j ]
    245 is the j-th order Taylor coefficient corresponding to z.
    246 \n
    247 \b Output: \a taylor [ ( \a i_z - 1 ) * \a nc_taylor + j ]
    248 is the j-th order Taylor coefficient corresponding to
     253\b Output: <code>taylor [ i_z * nc_taylor + k ]</code>,
     254for k = q , ... , p,
     255is the k-th order Taylor coefficient corresponding to z.
     256\n
     257\b Output: <code>taylor [ ( i_z - 1 ) * nc_taylor + k ]</code>,
     258for k = q , ... , p,
     259is the k-th order Taylor coefficient corresponding to
    249260the autillary result y.
    250261
     
    252263\li NumArg(op) == 1
    253264\li NumRes(op) == 2
    254 \li \a i_x + 1 < \a i_z
    255 \li \a j < \a nc_taylor
     265\li i_x + 1 < i_z
     266\li p < nc_taylor
     267\li q <= p
    256268*/
    257269template <class Base>
    258270inline void forward_unary2_op(
    259         size_t j           ,
     271        size_t q           ,
     272        size_t p           ,
    260273        size_t i_z         ,
    261274        size_t i_x         ,
     
    419432\a Base.
    420433
    421 \param d
    422 order of the Taylor coefficient that we are computing.
     434\param q
     435lowest order of the Taylor coefficient that we are computing.
     436
     437\param p
     438highest order of the Taylor coefficient that we are computing.
    423439
    424440\param i_z
     
    446462
    447463\param taylor
    448 \b Input: If x is a variable, \a taylor [ \a arg[0] * \a nc_taylor + k ]
    449 for k = 0 , ... , \a d
     464\b Input: If x is a variable,
     465<code>taylor [ arg[0] * nc_taylor + k ]</code>,
     466for k = 0 , ... , p,
    450467is the k-th order Taylor coefficient corresponding to x.
    451468\n
    452 \b Input: If y is a variable, \a taylor [ \a arg[1] * \a nc_taylor + k ]
    453 for k = 0 , ... , \a d
     469\b Input: If y is a variable,
     470<code>taylor [ arg[1] * nc_taylor + k ]</code>,
     471for k = 0 , ... , p,
    454472is the k-th order Taylor coefficient corresponding to y.
    455473\n
    456 \b Input: \a taylor [ \a i_z * \a nc_taylor + k ]
    457 for k = 0 , ... , \a d - 1
     474\b Input: <code>taylor [ i_z * nc_taylor + k ]</code>,
     475for k = 0 , ... , q-1,
    458476is the k-th order Taylor coefficient corresponding to z.
    459477\n
    460 \b Output: \a taylor [ \a i_z * \a nc_taylor + \a d ]
    461 is the \a d-th order Taylor coefficient corresponding to z.
     478\b Output: <code>taylor [ i_z * nc_taylor + k ]</code>,
     479for k = q, ... , p,
     480is the k-th order Taylor coefficient corresponding to z.
    462481
    463482\par Checked Assertions
    464483\li NumArg(op) == 2
    465484\li NumRes(op) == 1
    466 \li If x is a variable, \a arg[0] < \a i_z
    467 \li If y is a variable, \a arg[1] < \a i_z
    468 \li \a d < \a nc_taylor
     485\li If x is a variable, arg[0] < i_z
     486\li If y is a variable, arg[1] < i_z
     487\li p <  nc_taylor
     488\li q <=  p
    469489*/
    470490template <class Base>
    471491inline void forward_binary_op(
    472         size_t         d         ,
     492        size_t        q          ,
     493        size_t        p          ,
    473494        size_t        i_z        ,
    474495        const addr_t* arg        ,
     
    661682\a Base.
    662683
    663 \param d
    664 order of the Taylor coefficient that we are computing.
     684\param q
     685lowest order of the Taylor coefficient that we are computing.
     686
     687\param p
     688highest order of the Taylor coefficient that we are computing.
    665689
    666690\param i_z
     
    696720
    697721\param taylor
    698 \b Input: If x is a variable, \a taylor [ \a arg[0] * \a nc_taylor + k ]
    699 for k = 0 , ... , \a d
     722\b Input: If x is a variable,
     723<code>taylor [ arg[0] * nc_taylor + k ]</code>
     724for k = 0 , ... , p,
    700725is the k-th order Taylor coefficient corresponding to x.
    701726\n
    702 \b Input: If y is a variable, \a taylor [ \a arg[1] * \a nc_taylor + k ]
    703 for k = 0 , ... , \a d
     727\b Input: If y is a variable,
     728<code>taylor [ arg[1] * nc_taylor + k ]</code>
     729for k = 0 , ... , p
    704730is the k-th order Taylor coefficient corresponding to y.
    705731\n
    706 \b Input: \a taylor [ \a (i_z - 2 + j) * \a nc_taylor + k ]
    707 for j = 0, 1, 2 , and for k = 0 , ... , \a d - 1
     732\b Input: <code>taylor [ (i_z-2+j) * nc_taylor + k ]</code>,
     733for j = 0, 1, 2 , for k = 0 , ... , q-1,
    708734is the k-th order Taylor coefficient corresponding to z_j.
    709735\n
    710 \b Output: \a taylor [ \a (i_z - 2 + j) * \a nc_taylor + \a d ]
    711 is the \a d-th order Taylor coefficient corresponding to z_j.
     736\b Output: <code>taylor [ (i_z-2+j) * nc_taylor + k ]</code>,
     737is the k-th order Taylor coefficient corresponding to z_j.
    712738
    713739\par Checked Assertions
    714740\li NumArg(op) == 2
    715741\li NumRes(op) == 3
    716 \li If x is a variable, \a arg[0] < \a i_z - 2
    717 \li If y is a variable, \a arg[1] < \a i_z - 2
    718 \li \a d < \a nc_taylor
     742\li If x is a variable, arg[0] < i_z - 2
     743\li If y is a variable, arg[1] < i_z - 2
     744\li p < nc_taylor
     745\li q <= p
    719746*/
    720747template <class Base>
    721748inline void forward_pow_op(
    722         size_t        d          ,
     749        size_t        q          ,
     750        size_t        p          ,
    723751        size_t        i_z        ,
    724752        const addr_t* arg        ,
  • trunk/cppad/local/rev_hes_sweep.hpp

    r2794 r2859  
    168168
    169169        // work space used by UserOp.
     170        vector<size_t>     user_ix;  // variable indices for argument vector x
    170171        typedef std::set<size_t> size_set;
    171         const size_t user_q = limit; // maximum element plus one
    172172        size_set::iterator set_itr;  // iterator for a standard set
    173173        size_set::iterator set_end;  // end of iterator sequence
    174         vector<size_t>     user_ix;  // variable indices for argument vector x
    175         vector< size_set > user_r;   // forward Jacobian sparsity pattern for x
     174        vector< size_set > set_r;    // forward Jacobian sparsity for x
     175        vector< size_set > set_u;    // reverse Hessian sparsity for y
     176        vector< size_set > set_v;    // reverse Hessian sparsity for x
     177        //
     178        vector<bool>       bool_r;   // bool forward Jacobian sparsity for x
     179        vector<bool>       bool_u;   // bool reverse Hessian sparsity for y
     180        vector<bool>       bool_v;   // bool reverse Hessian sparsity for x
     181        //
     182        vector<bool>       user_vx;  // which components of x are variables
    176183        vector<bool>       user_s;   // reverse Jacobian sparsity for y
    177184        vector<bool>       user_t;   // reverse Jacobian sparsity for x
    178         vector< size_set > user_u;   // reverse Hessian sparsity for y
    179         vector< size_set > user_v;   // reverse Hessian sparsity for x
    180         size_t user_index = 0;       // indentifier for this user_atomic operation
     185        const size_t user_q = limit; // maximum element plus one
     186        size_t user_index = 0;       // indentifier for this atomic operation
    181187        size_t user_id    = 0;       // user identifier for this call to operator
    182188        size_t user_i     = 0;       // index in result vector
     
    184190        size_t user_m     = 0;       // size of result vector
    185191        size_t user_n     = 0;       // size of arugment vector
     192        //
     193        atomic_base<Base>* user_atom = CPPAD_NULL; // user's atomic op calculator
     194        bool               user_bool = false;      // use bool or set sparsity ?
     195# ifndef NDEBUG
     196        bool               user_ok   = false;      // atomic op return value
     197# endif
    186198        // next expected operator in a UserOp sequence
    187199        enum { user_start, user_arg, user_ret, user_end } user_state = user_end;
     
    564576                                user_n     = arg[2];
    565577                                user_m     = arg[3];
    566                                 if(user_ix.size() < user_n)
    567                                 {       user_ix.resize(user_n);
    568                                         user_r.resize(user_n);
    569                                         user_t.resize(user_n);
    570                                         user_v.resize(user_n);
     578                                user_atom  = atomic_base<Base>::list(user_index);
     579                                user_bool  = user_atom->sparsity() ==
     580                                                        atomic_base<Base>::bool_sparsity_enum;
     581                                user_ix.resize(user_n);
     582                                user_vx.resize(user_n);
     583                                user_s.resize(user_m);
     584                                user_t.resize(user_n);
     585
     586                                // simpler to initialize all sparsity patterns as empty
     587                                for(i = 0; i < user_m; i++)
     588                                        user_s[i] = false;
     589                                for(i = 0; i < user_n; i++)
     590                                        user_t[i] = false;
     591                                if( user_bool )
     592                                {       bool_r.resize(user_n * user_q);
     593                                        bool_u.resize(user_m * user_q);
     594                                        bool_v.resize(user_n * user_q);
     595                                        // simpler to initialize all patterns as empty
     596                                        for(i = 0; i < user_m; i++)
     597                                        {
     598                                                for(j = 0; j < user_q; j++)
     599                                                        bool_u[ i * user_q + j] = false;
     600                                        }
     601                                        for(i = 0; i < user_n; i++)
     602                                        {
     603                                                for(j = 0; j < user_q; j++)
     604                                                {       bool_r[ i * user_q + j] = false;
     605                                                        bool_v[ i * user_q + j] = false;
     606                                                }
     607                                        }
    571608                                }
    572                                 if(user_s.size() < user_m)
    573                                 {       user_s.resize(user_m);
    574                                         user_u.resize(user_m);
     609                                else
     610                                {       set_r.resize(user_n);
     611                                        set_u.resize(user_m);
     612                                        set_v.resize(user_n);
     613                                        for(i = 0; i < user_m; i++)
     614                                                set_u[i].clear();
     615                                        for(i = 0; i < user_n; i++)
     616                                        {       set_r[i].clear();
     617                                                set_v[i].clear();
     618                                        }
    575619                                }
    576620                                user_j     = user_n;
     
    587631
    588632                                // call users function for this operation
    589                                 user_atomic<Base>::rev_hes_sparse(user_index, user_id,
    590                                         user_n, user_m,
    591                                         user_q, user_r, user_s, user_t, user_u, user_v
     633                                user_atom->set_id(user_id);
     634# ifdef NDEBUG
     635                                if( user_bool )
     636                                        user_atom->rev_sparse_hes(user_vx,
     637                                                user_s, user_t, user_q, bool_r, bool_u, bool_v
    592638                                );
    593                                 for(j = 0; j < user_n; j++) if( user_ix[j] > 0 )
    594                                 {       size_t i_x = user_ix[j];
    595                                         RevJac[i_x] = user_t[j];
    596                                         set_itr = user_v[j].begin();
    597                                         set_end = user_v[j].end();
    598                                         while( set_itr != set_end )
    599                                                 rev_hes_sparse.add_element(i_x, *set_itr++);
     639                                else
     640                                        user_atom->rev_sparse_hes(user_vx,
     641                                                user_s, user_t, user_q, set_r, set_u, set_v
     642                                );
     643# else
     644                                if( user_bool )
     645                                        user_ok = user_atom->rev_sparse_hes(user_vx,
     646                                                user_s, user_t, user_q, bool_r, bool_u, bool_v
     647                                );
     648                                else
     649                                        user_ok = user_atom->rev_sparse_hes(user_vx,
     650                                                user_s, user_t, user_q, set_r, set_u, set_v
     651                                );
     652                                if( ! user_ok )
     653                                {       std::string msg = user_atom->afun_name()
     654                                        + ": atomic_base.rev_sparse_hes: returned false";
     655                                        CPPAD_ASSERT_KNOWN(false, msg.c_str() );
     656                                }
     657# endif
     658                                for(i = 0; i < user_n; i++) if( user_ix[i] > 0 )
     659                                {
     660                                        size_t  i_x = user_ix[i];
     661                                        if( user_t[i] )
     662                                                RevJac[i_x] = true;
     663                                        if( user_bool )
     664                                        {
     665                                                for(j = 0; j < user_q; j++)
     666                                                        if( bool_v[ i * user_q + j ] )
     667                                                                rev_hes_sparse.add_element(i_x, j);
     668                                        }
     669                                        else
     670                                        {
     671                                                set_itr = set_v[i].begin();
     672                                                set_end = set_v[i].end();
     673                                                while( set_itr != set_end )
     674                                                        rev_hes_sparse.add_element(i_x, *set_itr++);
     675                                        }
    600676                                }
    601677               }
     
    610686                        --user_j;
    611687                        user_ix[user_j] = 0;
    612                         user_r[user_j].clear();
     688                        user_vx[user_j] = false;
    613689                        if( user_j == 0 )
    614690                                user_state = user_start;
     
    624700                        --user_j;
    625701                        user_ix[user_j] = arg[0];
    626                         user_r[user_j].clear();
     702                        user_vx[user_j] = true;
    627703                        for_jac_sparse.begin(arg[0]);
    628704                        i = for_jac_sparse.next_element();
    629705                        while( i < user_q )
    630                         {       user_r[user_j].insert(i);
     706                        {       if( user_bool )
     707                                        bool_r[ user_j * user_q + i ] = true;
     708                                else
     709                                        set_r[user_j].insert(i);
    631710                                i = for_jac_sparse.next_element();
    632711                        }
     
    642721                        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_par );
    643722                        --user_i;
    644                         user_s[user_i] = false;
    645                         user_u[user_i].clear();
    646723                        if( user_i == 0 )
    647724                                user_state = user_arg;
     
    653730                        CPPAD_ASSERT_UNKNOWN( 0 < user_i && user_i <= user_m );
    654731                        --user_i;
    655                         user_s[user_i] = RevJac[i_var];
    656                         user_u[user_i].clear();
     732                        if( RevJac[i_var] )
     733                        {
     734                                user_s[user_i] = true;
     735                        }
    657736                        rev_hes_sparse.begin(i_var);
    658                         i = rev_hes_sparse.next_element();
    659                         while( i < user_q )
    660                         {       user_u[user_i].insert(i);
    661                                 i = rev_hes_sparse.next_element();
     737                        j = rev_hes_sparse.next_element();
     738                        while( j < user_q )
     739                        {       if( user_bool )
     740                                        bool_u[user_i * user_q + j] = true;
     741                                else
     742                                        set_u[user_i].insert(j);
     743                                j = rev_hes_sparse.next_element();
    662744                        }
    663745                        if( user_i == 0 )
  • trunk/cppad/local/rev_jac_sweep.hpp

    r2794 r2859  
    2929*/
    3030# define CPPAD_REV_JAC_SWEEP_TRACE 0
     31
     32/*
     33\def CPPAD_ATOMIC_CALL
     34This avoids warnings when NDEBUG is defined and user_ok is not used.
     35If \c NDEBUG is defined, this resolves to
     36\code
     37        user_atom->rev_sparse_jac
     38\endcode
     39otherwise, it respolves to
     40\code
     41        user_ok = user_atom->rev_sparse_jac
     42\endcode
     43This maco is undefined at the end of this file to facillitate is
     44use with a different definition in other files.
     45*/
     46# ifdef NDEBUG
     47# define CPPAD_ATOMIC_CALL user_atom->rev_sparse_jac
     48# else
     49# define CPPAD_ATOMIC_CALL user_ok = user_atom->rev_sparse_jac
     50# endif
    3151
    3252/*!
     
    137157        // work space used by UserOp.
    138158        typedef std::set<size_t> size_set;
    139         const size_t user_q = limit; // maximum element plus one
    140159        size_set::iterator set_itr;  // iterator for a standard set
    141160        size_set::iterator set_end;  // end of iterator sequence
    142         vector< size_set > user_r;   // sparsity pattern for the argument x
    143         vector< size_set > user_s;   // sparisty pattern for the result y
    144         size_t user_index = 0;       // indentifier for this user_atomic operation
     161        vector< size_set > set_r;   // set sparsity pattern for the argument x
     162        vector< size_set > set_s;   // set sparisty pattern for the result y
     163        //
     164        vector<bool>       bool_r;   // bool sparsity pattern for the argument x
     165        vector<bool>       bool_s;   // bool sparisty pattern for the result y
     166        //
     167        const size_t user_q = limit; // maximum element plus one
     168        size_t user_index = 0;       // indentifier for this atomic operation
    145169        size_t user_id    = 0;       // user identifier for this call to operator
    146170        size_t user_i     = 0;       // index in result vector
     
    148172        size_t user_m     = 0;       // size of result vector
    149173        size_t user_n     = 0;       // size of arugment vector
     174        //
     175        atomic_base<Base>* user_atom = CPPAD_NULL; // user's atomic op calculator
     176        bool               user_bool = false;      // use bool or set sparsity ?
     177# ifndef NDEBUG
     178        bool               user_ok   = false;      // atomic op return value
     179# endif
    150180        // next expected operator in a UserOp sequence
    151181        enum { user_start, user_arg, user_ret, user_end } user_state = user_end;
     
    515545                                user_n     = arg[2];
    516546                                user_m     = arg[3];
    517                                 if(user_r.size() < user_n )
    518                                         user_r.resize(user_n);
    519                                 if(user_s.size() < user_m )
    520                                         user_s.resize(user_m);
     547                                user_atom  = atomic_base<Base>::list(user_index);
     548                                user_bool  = user_atom->sparsity() ==
     549                                                        atomic_base<Base>::bool_sparsity_enum;
     550                                if( user_bool )
     551                                {       if( bool_r.size() != user_m * user_q )
     552                                                bool_r.resize( user_m * user_q );
     553                                        if( bool_s.size() != user_n * user_q )
     554                                                bool_s.resize( user_n * user_q );
     555                                        for(i = 0; i < user_m; i++)
     556                                                for(j = 0; j < user_q; j++)
     557                                                        bool_r[ i * user_q + j] = false;
     558                                }
     559                                else
     560                                {       if(set_r.size() != user_m )
     561                                                set_r.resize(user_m);
     562                                        if(set_s.size() != user_n )
     563                                                set_s.resize(user_n);
     564                                        for(i = 0; i < user_m; i++)
     565                                                set_r[i].clear();
     566                                }
    521567                                user_j     = user_n;
    522568                                user_i     = user_m;
     
    529575                                CPPAD_ASSERT_UNKNOWN( user_n     == size_t(arg[2]) );
    530576                                CPPAD_ASSERT_UNKNOWN( user_m     == size_t(arg[3]) );
     577# ifndef NDEBUG
     578                                if( ! user_ok )
     579                                {       std::string msg = user_atom->afun_name()
     580                                        + ": atomic_base.rev_sparse_jac: returned false";
     581                                        CPPAD_ASSERT_KNOWN(false, msg.c_str() );
     582                                }
     583# endif
    531584                                user_state = user_end;
    532585               }
     
    554607                        // It might be faster if we add set union to var_sparsity
    555608                        // where one of the sets is not in var_sparsity.
    556                         set_itr = user_r[user_j].begin();
    557                         set_end = user_r[user_j].end();
    558                         while( set_itr != set_end )
    559                                 var_sparsity.add_element(arg[0], *set_itr++);   
     609                        if( user_bool )
     610                        {       for(j = 0; j < user_q; j++)
     611                                        if( bool_s[ user_j * user_q + j ] )     
     612                                                var_sparsity.add_element(arg[0], j);   
     613                        }
     614                        else
     615                        {       set_itr = set_s[user_j].begin();
     616                                set_end = set_s[user_j].end();
     617                                while( set_itr != set_end )
     618                                        var_sparsity.add_element(arg[0], *set_itr++);   
     619                        }
    560620                        if( user_j == 0 )
    561621                                user_state = user_start;
     
    569629                        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_par );
    570630                        --user_i;
    571                         user_s[user_i].clear();
    572631                        if( user_i == 0 )
    573632                        {       // call users function for this operation
    574                                 user_atomic<Base>::rev_jac_sparse(user_index, user_id,
    575                                         user_n, user_m, user_q, user_r, user_s
     633                                user_atom->set_id(user_id);
     634                                if( user_bool)
     635                                        CPPAD_ATOMIC_CALL(
     636                                                user_q, bool_r, bool_s
     637                                );
     638                                else
     639                                        CPPAD_ATOMIC_CALL(
     640                                                user_q, set_r, set_s
    576641                                );
    577642                                user_state = user_arg;
     
    584649                        CPPAD_ASSERT_UNKNOWN( 0 < user_i && user_i <= user_m );
    585650                        --user_i;
    586                         user_s[user_i].clear();
    587651                        var_sparsity.begin(i_var);
    588652                        i = var_sparsity.next_element();
    589653                        while( i < user_q )
    590                         {       user_s[user_i].insert(i);
     654                        {       if( user_bool )
     655                                        bool_r[ user_i * user_q + i ] = true;
     656                                        else
     657                                                set_r[user_i].insert(i);
    591658                                i = var_sparsity.next_element();
    592659                        }
    593660                        if( user_i == 0 )
    594661                        {       // call users function for this operation
    595                                 user_atomic<Base>::rev_jac_sparse(user_index, user_id,
    596                                         user_n, user_m, user_q, user_r, user_s
     662                                user_atom->set_id(user_id);
     663                                if( user_bool)
     664                                        CPPAD_ATOMIC_CALL(
     665                                                user_q, bool_r, bool_s
     666                                );
     667                                else
     668                                        CPPAD_ATOMIC_CALL(
     669                                                user_q, set_r, set_s
    597670                                );
    598671                                user_state = user_arg;
     
    637710// preprocessor symbols that are local to this file
    638711# undef CPPAD_REV_JAC_SWEEP_TRACE
     712# undef CPPAD_ATOMIC_CALL
    639713
    640714# endif
  • trunk/cppad/local/rev_sparse_hes.hpp

    r2625 r2859  
    44
    55/* --------------------------------------------------------------------------
    6 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
     6CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-13 Bradley M. Bell
    77
    88CppAD is distributed under multiple licenses. This distribution is under
     
    3838
    3939$head Syntax$$
    40 $icode%h% = %f%.RevSparseHes(%q%, %s%)%$$
     40$icode%h% = %f%.RevSparseHes(%q%, %s%)
     41%$$
     42$icode%h% = %f%.RevSparseHes(%q%, %s%, %transpose%)%$$
    4143
    4244
     
    4446We use $latex F : \B{R}^n \rightarrow \B{R}^m$$ to denote the
    4547$cref/AD function/glossary/AD Function/$$ corresponding to $icode f$$.
    46 We define the matrix $latex H(x) \in \B{R}^{q \times n}$$
    47 as the partial with respect to $latex x$$,
    48 of the partial with respect to $latex u$$ (at $latex u = 0$$),
    49 of $latex S * F[ x + R * u ]$$ where
    50 $latex R \in \B{R}^{n \times q}$$ and $latex S \in \B{R}^{1 \times m}$$; i.e.,
     48For a fixed matrix $latex R \in \B{R}^{n \times q}$$
     49and a fixed vector $latex S \in \B{R}^{1 \times m}$$,
     50we define
    5151$latex \[
    52         H(x)  =  R^\R{T} (S * F)^{(2)} ( x )
     52\begin{array}{rcl}
     53H(x)
     54& = & \partial_x \left[ \partial_u S * F[ x + R * u ] \right]_{u=0}
     55\\
     56& = & R^\R{T} * (S * F)^{(2)} ( x )
     57\\
     58H(x)^\R{T}
     59& = & (S * F)^{(2)} ( x ) * R
     60\end{array}
    5361\] $$
    5462Given a
     
    7482%$$
    7583It specifies the number of columns in $latex R \in \B{R}^{n \times q}$$
    76 and the number of rows in
    77 $latex H(x) \in \B{R}^{q \times n}$$.
     84and the number of rows in $latex H(x) \in \B{R}^{q \times n}$$.
    7885It must be the same value as in the previous $cref ForSparseJac$$ call
    7986$codei%
    80         %f%.ForSparseJac(%q%, %r%)
    81 %$$
     87        %f%.ForSparseJac(%q%, %r%, %r_transpose%)
     88%$$
     89Note that if $icode r_transpose$$ is true, $icode r$$ in the call above
     90corresponding to $latex R^\R{T} \in \B{R}^{q \times n}$$
     91
     92$head transpose$$
     93The argument $icode transpose$$ has prototype
     94$codei%
     95        bool %transpose%
     96%$$
     97The default value $code false$$ is used when $icode transpose$$ is not present.
     98
    8299
    83100$head r$$
    84 The argument $icode r$$ in the previous call
     101The matrix $latex R$$ is specified by the previous call
    85102$codei%
    86         %f%.ForSparseJac(%q%, %r%)
    87 %$$
    88 is a $cref/sparsity pattern/glossary/Sparsity Pattern/$$
    89 for the matrix $latex R$$ above.
     103        %f%.ForSparseJac(%q%, %r%, %transpose%)
     104%$$
     105see $cref/r/ForSparseJac/r/$$.
    90106The type of the elements of
    91107$cref/VectorSet/RevSparseHes/VectorSet/$$ must be the
     
    113129%$$
    114130(see $cref/VectorSet/RevSparseHes/VectorSet/$$ below).
    115 If it has elements of type $code bool$$,
     131
     132$subhead transpose false$$
     133If $icode h$$ has elements of type $code bool$$,
    116134its size is $latex q * n$$.
    117135If it has elements of type $code std::set<size_t>$$,
    118 its size is $latex q$$.
     136its size is $latex q$$ and all the set elements are between
     137zero and $icode%n%-1%$$ inclusive.
    119138It specifies a
    120139$cref/sparsity pattern/glossary/Sparsity Pattern/$$
    121140for the matrix $latex H(x)$$.
     141
     142$subhead transpose true$$
     143If $icode h$$ has elements of type $code bool$$,
     144its size is $latex n * q$$.
     145If it has elements of type $code std::set<size_t>$$,
     146its size is $latex n$$ and all the set elements are between
     147zero and $icode%q%-1%$$ inclusive.
     148It specifies a
     149$cref/sparsity pattern/glossary/Sparsity Pattern/$$
     150for the matrix $latex H(x)^\R{T}$$.
    122151
    123152$head VectorSet$$
     
    189218is either \c sparse_pack, \c sparse_set, or \c sparse_list.
    190219
     220\param transpose
     221is true (false) if \c is is equal to \f$ H(x) \f$ (\f$ H(x)^T \f$)
     222where
     223\f[
     224        H(x) = R^T (S * F)^{(2)} (x)
     225\f]
     226where \f$ F \f$ is the function corresponding to the operation sequence
     227and \a x is any argument value.
     228
    191229\param q
    192230is the value of \a q in the
     
    206244the input value of \a h must be a vector with size \c q*n.
    207245The input value of its elements does not matter.
    208 On output, \a h is the sparsity pattern for the matrix
    209 \f[
    210         H(x) = R^T (S * F)^{(2)} (x)
    211 \f]
    212 where \f$ F \f$ is the function corresponding to the operation sequence
    213 and \a x is any argument value.
     246On output, \a h is the sparsity pattern for the matrix \f$ H(x) \f$
     247or \f$ H(x)^T \f$ depending on \c transpose.
    214248
    215249\param total_num_var
     
    236270template <class Base, class VectorSet, class Sparsity>
    237271void RevSparseHesBool(
     272        bool                      transpose         ,
    238273        size_t                    q                 ,
    239274        const VectorSet&          s                 ,
     
    257292        CPPAD_ASSERT_KNOWN(
    258293                q == for_jac_sparsity.end(),
    259                 "RevSparseHes: q (first argument) is not equal to its value\n"
     294                "RevSparseHes: q is not equal to its value\n"
    260295                "in the previous call to ForSparseJac with this ADFun object."
    261296        );
    262297        CPPAD_ASSERT_KNOWN(
    263298                size_t(s.size()) == m,
    264                 "RevSparseHes: s (second argument) length is not equal to\n"
     299                "RevSparseHes: size of s is not equal to\n"
    265300                "range dimension for ADFun object."
    266301        );
     
    277312        }
    278313
    279 
    280314        // vector of sets that will hold reverse Hessain values
    281315        Sparsity       rev_hes_sparsity;
     
    294328        // return values corresponding to independent variables
    295329        CPPAD_ASSERT_UNKNOWN( size_t(h.size()) == n * q );
     330        for(j = 0; j < n; j++)
     331        {       for(i = 0; i < q; i++)
     332                {       if( transpose )
     333                                h[ j * q + i ] = false;
     334                        else    h[ i * n + j ] = false;
     335                }
     336        }
    296337
    297338        // j is index corresponding to reverse mode partial
     
    304345
    305346                // extract the result from rev_hes_sparsity
    306                 for(i = 0; i < q; i++)
    307                         h[ i * n + j ] = false;
    308347                CPPAD_ASSERT_UNKNOWN( rev_hes_sparsity.end() == q );
    309348                rev_hes_sparsity.begin(j + 1);
    310349                i = rev_hes_sparsity.next_element();
    311350                while( i < q )
    312                 {       h[ i * n + j ] = true;
     351                {       if( transpose )
     352                                h[ j * q + i ] = true;
     353                        else    h[ i * n + j ] = true;
    313354                        i = rev_hes_sparsity.next_element();
    314355                }
     
    334375\tparam Sparsity
    335376is either \c sparse_pack, \c sparse_set, or \c sparse_list.
     377
     378\param transpose
     379is true (false) if \c is is equal to \f$ H(x) \f$ (\f$ H(x)^T \f$)
     380where
     381\f[
     382        H(x) = R^T (S * F)^{(2)} (x)
     383\f]
     384where \f$ F \f$ is the function corresponding to the operation sequence
     385and \a x is any argument value.
    336386
    337387\param q
     
    350400
    351401\param h
    352 the input value of \a h must be a vector with size \a q.
     402If \c transpose, the input value of \a h must be a vector with size \a q.
     403Otherwise, its input value must have size \c n;
    353404On input, each element of \a h must be an empty set.
    354 The input value of its elements does not matter.
    355 On output, \a h is the sparsity pattern for the matrix
    356 \f[
    357         H(x) = R^T (S * F)^{(2)} (x)
    358 \f]
    359 where \f$ F \f$ is the function corresponding to the operation sequence
    360 and \a x is any argument value.
     405On output, \a h is the sparsity pattern for the matrix \f$ H(x) \f$
     406or \f$ H(x)^T \f$ depending on \c transpose.
    361407
    362408\param total_num_var
     
    383429template <class Base, class VectorSet, class Sparsity>
    384430void RevSparseHesSet(
     431        bool                      transpose         ,
    385432        size_t                    q                 ,
    386433        const VectorSet&          s                 ,
     
    409456        CPPAD_ASSERT_KNOWN(
    410457                q == for_jac_sparsity.end(),
    411                 "RevSparseHes: q (first argument) is not equal to its value\n"
     458                "RevSparseHes: q is not equal to its value\n"
    412459                "in the previous call to ForSparseJac with this ADFun object."
    413460        );
    414461        CPPAD_ASSERT_KNOWN(
    415462                s.size() == 1,
    416                 "RevSparseHes: s (second argument) length is not equal to one."
     463                "RevSparseHes: size of s is not equal to one."
    417464        );
    418465
     
    452499        // return values corresponding to independent variables
    453500        // j is index corresponding to reverse mode partial
    454         CPPAD_ASSERT_UNKNOWN( size_t(h.size()) == q );
     501        CPPAD_ASSERT_UNKNOWN( size_t(h.size()) == q || transpose );
     502        CPPAD_ASSERT_UNKNOWN( size_t(h.size()) == n || ! transpose );
    455503        for(j = 0; j < n; j++)
    456504        {       CPPAD_ASSERT_UNKNOWN( ind_taddr[j] < total_num_var );
     
    464512                i = rev_hes_sparsity.next_element();
    465513                while( i < q )
    466                 {       h[i].insert(j);
     514                {       if( transpose )
     515                                h[j].insert(i);
     516                        else    h[i].insert(j);
    467517                        i = rev_hes_sparsity.next_element();
    468518                }
     
    488538is a simple vector with elements of type \c bool
    489539or \c std::set<size_t>.
     540
     541\param transpose
     542is true (false) if \c is is equal to \f$ H(x) \f$ (\f$ H(x)^T \f$)
     543where
     544\f[
     545        H(x) = R^T (S * F)^{(2)} (x)
     546\f]
     547where \f$ F \f$ is the function corresponding to the operation sequence
     548and \a x is any argument value.
    490549
    491550\param q
     
    506565
    507566\return
    508 is a vector with size \c q*n.
    509 containing a sparsity pattern for the matrix
     567If \c transpose is false (true),
     568the return vector is a sparsity pattern for \f$ H(x) \f$ (\f$ H(x)^T \f$).
    510569\f[
    511570        H(x) = R^T ( S * F)^{(2)} (x)
     
    517576template <class Base>
    518577template <class VectorSet>
    519 VectorSet ADFun<Base>::RevSparseHes(size_t q,  const VectorSet& s)
     578VectorSet ADFun<Base>::RevSparseHes(
     579        size_t q,  const VectorSet& s, bool transpose
     580)
    520581{       VectorSet h;
    521582        typedef typename VectorSet::value_type Set_type;
     
    525586        RevSparseHesCase(
    526587                Set_type()    ,
     588                transpose     ,
    527589                q             ,
    528590                s             ,
     
    542604is a \c bool value. This argument is used to dispatch to the proper source
    543605code depending on the vlaue of \c VectorSet::value_type.
     606
     607\param transpose
     608See \c RevSparseHes(q, s).
    544609
    545610\param q
     
    556621void ADFun<Base>::RevSparseHesCase(
    557622        bool              set_type         ,
     623        bool              transpose        , 
    558624        size_t            q                , 
    559625        const VectorSet&  s                ,
     
    572638        // use sparse_pack for the calculation
    573639        CppAD::RevSparseHesBool(
     640                transpose                ,
    574641                q                        ,
    575642                s                        ,
     
    594661code depending on the vlaue of \c VectorSet::value_type.
    595662
     663\param transpose
     664See \c RevSparseHes(q, s).
     665
    596666\param q
    597667See \c RevSparseHes(q, s).
     
    607677void ADFun<Base>::RevSparseHesCase(
    608678        const std::set<size_t>&   set_type         ,
     679        bool                      transpose        , 
    609680        size_t                    q                , 
    610681        const VectorSet&          s                ,
    611682        VectorSet&                h                )
    612 {       h.resize(q);
     683{       size_t n = Domain();
     684        if( transpose )
     685                h.resize(n);
     686        else    h.resize(q);
    613687
    614688        CPPAD_ASSERT_KNOWN(
     
    622696        // use sparse_pack for the calculation
    623697        CppAD::RevSparseHesSet(
     698                transpose                ,
    624699                q                        ,
    625700                s                        ,
  • trunk/cppad/local/rev_sparse_jac.hpp

    r2811 r2859  
    3636
    3737$head Syntax$$
    38 $icode%r% = %F%.RevSparseJac(%p%, %s%)%$$
     38$icode%s% = %f%.RevSparseJac(%q%, %r%)
     39%$$
     40$icode%s% = %f%.RevSparseJac(%q%, %r%, %transpose%)%$$
    3941
    4042$head Purpose$$
    41 We use $latex F : \B{R}^n \rightarrow R^m$$ to denote the
     43We use $latex F : B^n \rightarrow B^m$$ to denote the
    4244$cref/AD function/glossary/AD Function/$$ corresponding to $icode f$$.
    43 For a fixed $latex p \times m$$ matrix $latex S$$,
    44 the Jacobian of $latex S * F( x )$$
     45For a fixed matrix $latex R \in B^{q \times m}$$,
     46the Jacobian of $latex R * F( x )$$
    4547with respect to $latex x$$ is
    4648$latex \[
    47         J(x) = S * F^{(1)} ( x )
     49        S(x) = R * F^{(1)} ( x )
    4850\] $$
    4951Given a
    5052$cref/sparsity pattern/glossary/Sparsity Pattern/$$
    51 for $latex S$$,
    52 $code RevSparseJac$$ returns a sparsity pattern for the $latex J(x)$$.
     53for $latex R$$,
     54$code RevSparseJac$$ returns a sparsity pattern for the $latex S(x)$$.
    5355
    5456$head f$$
     
    6062$head x$$
    6163the sparsity pattern is valid for all values of the independent
    62 variables in $latex x \in \B{R}^n$$
     64variables in $latex x \in B^n$$
    6365(even if it has $cref CondExp$$ or $cref VecAD$$ operations).
    6466
    65 $head p$$
    66 The argument $icode p$$ has prototype
     67$head q$$
     68The argument $icode q$$ has prototype
    6769$codei%
    68         size_t %p%
     70        size_t %q%
    6971%$$
    7072It specifies the number of rows in
    71 $latex S \in \B{R}^{p \times m}$$ and the
    72 Jacobian $latex J(x) \in \B{R}^{p \times n}$$.
    73 
    74 $head s$$
     73$latex R \in B^{q \times m}$$ and the
     74Jacobian $latex S(x) \in B^{q \times n}$$.
     75
     76$head transpose$$
     77The argument $icode transpose$$ has prototype
     78$codei%
     79        bool %transpose%
     80%$$
     81The default value $code false$$ is used when $icode transpose$$ is not present.
     82
     83$head r$$
    7584The argument $icode s$$ has prototype
    7685$codei%
    77         const %VectorSet%& %s%
     86        const %VectorSet%& %r%
    7887%$$
    79 (see $cref/VectorSet/RevSparseJac/VectorSet/$$ below).
    80 If it has elements of type $code bool$$,
    81 its size is $latex p * m$$.
     88see $cref/VectorSet/RevSparseJac/VectorSet/$$ below.
     89
     90$subhead transpose false$$
     91If $icode r$$ has elements of type $code bool$$,
     92its size is $latex q * m$$.
    8293If it has elements of type $code std::set<size_t>$$,
    83 its size is $icode p$$ and all its set elements are between
     94its size is $icode q$$ and all its set elements are between
    8495zero and $latex m - 1$$.
    8596It specifies a
    8697$cref/sparsity pattern/glossary/Sparsity Pattern/$$
    87 for the matrix $icode S$$.
    88 
    89 $head r$$
    90 The return value $icode r$$ has prototype
    91 $codei%
    92         %VectorSet% %r%
    93 %$$
    94 (see $cref/VectorSet/RevSparseJac/VectorSet/$$ below).
    95 If it has elements of type $code bool$$,
    96 its size is $latex p * n$$.
     98for the matrix $latex R \in B^{q \times m}$$.
     99
     100$subhead transpose true$$
     101If $icode r$$ has elements of type $code bool$$,
     102its size is $latex m * q$$.
    97103If it has elements of type $code std::set<size_t>$$,
    98 its size is $icode p$$.
     104its size is $icode m$$ and all its set elements are between
     105zero and $latex q - 1$$.
    99106It specifies a
    100107$cref/sparsity pattern/glossary/Sparsity Pattern/$$
    101 for the matrix $latex J(x)$$.
     108for the matrix $latex R^\R{T} \in B^{m \times q}$$.
     109
     110$head s$$
     111The return value $icode s$$ has prototype
     112$codei%
     113        %VectorSet% %s%
     114%$$
     115see $cref/VectorSet/RevSparseJac/VectorSet/$$ below.
     116
     117$subhead transpose false$$
     118If it has elements of type $code bool$$,
     119its size is $latex q * n$$.
     120If it has elements of type $code std::set<size_t>$$,
     121its size is $icode q$$ and all its set elements are between
     122zero and $latex n - 1$$.
     123It specifies a
     124$cref/sparsity pattern/glossary/Sparsity Pattern/$$
     125for the matrix $latex S(x) \in {q \times n}$$.
     126
     127$subhead transpose true$$
     128If it has elements of type $code bool$$,
     129its size is $latex n * q$$.
     130If it has elements of type $code std::set<size_t>$$,
     131its size is $icode n$$ and all its set elements are between
     132zero and $latex q - 1$$.
     133It specifies a
     134$cref/sparsity pattern/glossary/Sparsity Pattern/$$
     135for the matrix $latex S(x)^\R{T} \in {n \times q}$$.
    102136
    103137$head VectorSet$$
     
    109143
    110144$head Entire Sparsity Pattern$$
    111 Suppose that $latex p = m$$ and
    112 $latex S$$ is the $latex m \times m$$ identity matrix.
     145Suppose that $latex q = m$$ and
     146$latex R$$ is the $latex m \times m$$ identity matrix.
    113147In this case,
    114 the corresponding value for $icode r$$ is a
    115 sparsity pattern for the Jacobian $latex J(x) = F^{(1)} ( x )$$.
     148the corresponding value for $icode s$$ is a
     149sparsity pattern for the Jacobian $latex S(x) = F^{(1)} ( x )$$.
    116150
    117151$head Example$$
     
    153187is a simple vector class with elements of type \c bool.
    154188
    155 \param p
    156 is the number of rows in the matrix \f$ S \f$.
     189\param transpose
     190are the sparsity patterns transposed.
     191
     192\param q
     193is the number of rows in the matrix \f$ R \f$.
     194
     195\param r
     196is a sparsity pattern for the matrix \f$ R \f$.
    157197
    158198\param s
    159 is a sparsity pattern for the matrix \f$ S \f$.
    160 
    161 \param r
    162 the input value of \a r must be a vector with size \c p*n
     199the input value of \a s must be a vector with size <tt>p*n</tt>
    163200where \c n is the number of independent variables
    164201corresponding to the operation sequence stored in \a play.
    165 The input value of the components of \c r does not matter.
    166 On output, \a r is the sparsity pattern for the matrix
     202The input value of the components of \c s does not matter.
     203On output, \a s is the sparsity pattern for the matrix
    167204\f[
    168         J(x) = S * F^{(1)} (x)
     205        S(x) = R * F^{(1)} (x)
    169206\f]
    170207where \f$ F \f$ is the function corresponding to the operation sequence
     
    189226template <class Base, class VectorSet>
    190227void RevSparseJacBool(
    191         size_t                 p                ,
    192         const VectorSet&       s                ,
    193         VectorSet&             r                ,
     228        bool                   transpose        ,
     229        size_t                 q                ,
     230        const VectorSet&       r                ,
     231        VectorSet&             s                ,
    194232        size_t                 total_num_var    ,
    195233        CppAD::vector<size_t>& dep_taddr        ,
     
    208246
    209247        CPPAD_ASSERT_KNOWN(
    210                 p > 0,
    211                 "RevSparseJac: p (first argument) is not greater than zero"
    212         );
    213 
     248                q > 0,
     249                "RevSparseJac: q is not greater than zero"
     250        );
    214251        CPPAD_ASSERT_KNOWN(
    215                 size_t(s.size()) == p * m,
    216                 "RevSparseJac: s (second argument) length is not equal to\n"
    217                 "p (first argument) times range dimension for ADFun object."
     252                size_t(r.size()) == q * m,
     253                "RevSparseJac: size of r is not equal to\n"
     254                "q times range dimension for ADFun object."
    218255        );
    219256
    220257        // vector of sets that will hold the results
    221258        sparse_pack    var_sparsity;
    222         var_sparsity.resize(total_num_var, p);
     259        var_sparsity.resize(total_num_var, q);
    223260
    224261        // The sparsity pattern corresponding to the dependent variables
    225262        for(i = 0; i < m; i++)
    226263        {       CPPAD_ASSERT_UNKNOWN( dep_taddr[i] < total_num_var );
    227 
    228                 for(j = 0; j < p; j++) if( s[ i * p + j ] )
    229                         var_sparsity.add_element( dep_taddr[i], j );
     264                if( transpose )
     265                {       for(j = 0; j < q; j++) if( r[ j * m + i ] )
     266                                var_sparsity.add_element( dep_taddr[i], j );
     267                }
     268                else
     269                {       for(j = 0; j < q; j++) if( r[ i * q + j ] )
     270                                var_sparsity.add_element( dep_taddr[i], j );
     271                }
    230272        }
    231273
     
    239281
    240282        // return values corresponding to dependent variables
    241         CPPAD_ASSERT_UNKNOWN( size_t(r.size()) == p * n );
     283        CPPAD_ASSERT_UNKNOWN( size_t(s.size()) == q * n );
    242284        for(j = 0; j < n; j++)
    243285        {       CPPAD_ASSERT_UNKNOWN( ind_taddr[j] == (j+1) );
     
    247289
    248290                // extract the result from var_sparsity
    249                 for(i = 0; i < p; i++)
    250                         r[ i * n + j ] = false;
    251                 CPPAD_ASSERT_UNKNOWN( var_sparsity.end() == p );
     291                if( transpose )
     292                {       for(i = 0; i < q; i++)
     293                                s[ j * q + i ] = false;
     294                }
     295                else
     296                {       for(i = 0; i < q; i++)
     297                                s[ i * n + j ] = false;
     298                }
     299                CPPAD_ASSERT_UNKNOWN( var_sparsity.end() == q );
    252300                var_sparsity.begin(j+1);
    253301                i = var_sparsity.next_element();
    254                 while( i < p )
    255                 {       r[ i * n + j ] = true;
    256                         i              = var_sparsity.next_element();
     302                while( i < q )
     303                {       if( transpose )
     304                                s[ j * q + i ] = true;
     305                        else    s[ i * n + j ] = true;
     306                        i  = var_sparsity.next_element();
    257307                }
    258308        }
     
    272322is a simple vector class with elements of type \c std::set<size_t>.
    273323
    274 \param p
     324\param transpose
     325see \c RevSparseJacBool.
     326
     327\param q
     328see \c RevSparseJacBool.
     329
     330\param r
    275331see \c RevSparseJacBool.
    276332
    277333\param s
    278 see \c RevSparseJacBool.
    279 
    280 \param r
    281334see \c RevSparseJacBool.
    282335
     
    295348template <class Base, class VectorSet>
    296349void RevSparseJacSet(
    297         size_t                 p                ,
    298         const VectorSet&       s                ,
    299         VectorSet&             r                ,
     350        bool                   transpose        ,
     351        size_t                 q                ,
     352        const VectorSet&       r                ,
     353        VectorSet&             s                ,
    300354        size_t                 total_num_var    ,
    301355        CppAD::vector<size_t>& dep_taddr        ,
     
    314368        // domain dimensions for F
    315369        size_t n = ind_taddr.size();
     370        size_t m = dep_taddr.size();
    316371
    317372        CPPAD_ASSERT_KNOWN(
    318                 p > 0,
    319                 "RevSparseJac: p (first argument) is not greater than zero"
    320         );
    321 
     373                q > 0,
     374                "RevSparseJac: q is not greater than zero"
     375        );
    322376        CPPAD_ASSERT_KNOWN(
    323                 size_t(s.size()) == p,
    324                 "RevSparseJac: s (second argument) length is not equal to "
    325                 "p (first argument)."
     377                size_t(r.size()) == q || transpose,
     378                "RevSparseJac: size of r is not equal to q and transpose is false."
     379        );
     380        CPPAD_ASSERT_KNOWN(
     381                size_t(r.size()) == m || ! transpose,
     382                "RevSparseJac: size of r is not equal to m and transpose is true."
    326383        );
    327384
    328385        // vector of lists that will hold the results
    329386        CPPAD_INTERNAL_SPARSE_SET    var_sparsity;
    330         var_sparsity.resize(total_num_var, p);
     387        var_sparsity.resize(total_num_var, q);
    331388
    332389        // The sparsity pattern corresponding to the dependent variables
    333         for(i = 0; i < p; i++)
    334         {       itr = s[i].begin();
    335                 while(itr != s[i].end())
    336                 {       j = *itr++;
    337                         CPPAD_ASSERT_KNOWN(
    338                                 j < dep_taddr.size(),
    339                                 "RevSparseJac: an element of the set s[i] "
    340                                 "has value greater than or equal Range dimension."
    341                         );
    342                         CPPAD_ASSERT_UNKNOWN( dep_taddr[j] < total_num_var );
    343                         var_sparsity.add_element( dep_taddr[j], i );
     390        if( transpose )
     391        {       for(i = 0; i < m; i++)
     392                {       itr = r[i].begin();
     393                        while(itr != r[i].end())
     394                        {       j = *itr++;
     395                                CPPAD_ASSERT_KNOWN(
     396                                j < q,
     397                                "RevSparseJac: transpose is true and element of the set\n"
     398                                "r[i] has value greater than or equal q."
     399                                );
     400                                CPPAD_ASSERT_UNKNOWN( dep_taddr[i] < total_num_var );
     401                                var_sparsity.add_element( dep_taddr[i], j );
     402                        }
    344403                }
    345404        }
    346 
     405        else
     406        {       for(i = 0; i < q; i++)
     407                {       itr = r[i].begin();
     408                        while(itr != r[i].end())
     409                        {       j = *itr++;
     410                                CPPAD_ASSERT_KNOWN(
     411                                j < m,
     412                                "RevSparseJac: transpose is false and element of the set\n"
     413                                "r[i] has value greater than or equal range dimension."
     414                                );
     415                                CPPAD_ASSERT_UNKNOWN( dep_taddr[j] < total_num_var );
     416                                var_sparsity.add_element( dep_taddr[j], i );
     417                        }
     418                }
     419        }
    347420        // evaluate the sparsity patterns
    348421        RevJacSweep(
     
    354427
    355428        // return values corresponding to dependent variables
    356         CPPAD_ASSERT_UNKNOWN( size_t(r.size()) == p );
     429        CPPAD_ASSERT_UNKNOWN( size_t(s.size()) == q || transpose );
     430        CPPAD_ASSERT_UNKNOWN( size_t(s.size()) == n || ! transpose );
    357431        for(j = 0; j < n; j++)
    358432        {       CPPAD_ASSERT_UNKNOWN( ind_taddr[j] == (j+1) );
     
    361435                CPPAD_ASSERT_UNKNOWN( play.GetOp( ind_taddr[j] ) == InvOp );
    362436
    363                 // extract result from rev_hes_sparsity
    364                 // and add corresponding elements to sets in r
    365                 CPPAD_ASSERT_UNKNOWN( var_sparsity.end() == p );
     437                CPPAD_ASSERT_UNKNOWN( var_sparsity.end() == q );
    366438                var_sparsity.begin(j+1);
    367439                i = var_sparsity.next_element();
    368                 while( i < p )
    369                 {       r[i].insert(j);
     440                while( i < q )
     441                {       if( transpose )
     442                                s[j].insert(i);
     443                        else    s[i].insert(j);
    370444                        i = var_sparsity.next_element();
    371445                }
     
    375449
    376450/*!
    377 Private helper function for \c RevSparseJac(p, s).
    378 
    379 All of the description in the public member function \c RevSparseJac(p, s)
    380 applies.
     451Private helper function for \c RevSparseJac(q, r, transpose).
     452
     453All of the description in the public member function
     454\c RevSparseJac(q, r, transpose) apply.
    381455
    382456\param set_type
    383457is a \c bool value.
    384 This arugment is used to dispatch to the proper source code
     458This argument is used to dispatch to the proper source code
    385459depending on the value of \c VectorSet::value_type.
    386460
    387 \param p
    388 See \c RevSparseJac(p, s).
     461\param transpose
     462See \c RevSparseJac(q, r, transpose)
     463
     464\param q
     465See \c RevSparseJac(q, r, transpose)
     466
     467\param r
     468See \c RevSparseJac(q, r, transpose)
    389469
    390470\param s
    391 See \c RevSparseJac(p, s).
    392 
    393 \param r
    394 is the return value for the corresponding call to RevSparseJac(p, s);
     471is the return value for the corresponding call to
     472RevSparseJac(q, r, transpose).
    395473*/
    396474
     
    399477void ADFun<Base>::RevSparseJacCase(
    400478        bool                set_type          ,
    401         size_t              p                 ,
    402         const VectorSet&    s                 ,
    403         VectorSet&          r                 )
     479        bool                transpose         ,
     480        size_t              q                 ,
     481        const VectorSet&    r                 ,
     482        VectorSet&          s                 )
    404483{       size_t n = Domain();
    405484
    406485        // dimension of the result vector
    407         r.resize( p * n );
    408 
    409         // store results in r
     486        s.resize( q * n );
     487
     488        // store results in s
    410489        RevSparseJacBool(
    411                 p              ,
     490                transpose      ,
     491                q              ,
     492                r              ,
    412493                s              ,
    413                 r              ,
    414494                total_num_var_ ,
    415495                dep_taddr_     ,
     
    420500
    421501/*!
    422 Private helper function for \c RevSparseJac(p, s).
    423 
    424 All of the description in the public member function \c RevSparseJac(p, s)
    425 applies.
     502Private helper function for \c RevSparseJac(q, r, transpose).
     503
     504All of the description in the public member function
     505\c RevSparseJac(q, r, transpose) apply.
    426506
    427507\param set_type
    428508is a \c std::set<size_t> object.
    429 This arugment is used to dispatch to the proper source code
     509This argument is used to dispatch to the proper source code
    430510depending on the value of \c VectorSet::value_type.
    431511
    432 \param p
    433 See \c RevSparseJac(p, s).
     512\param transpose
     513See \c RevSparseJac(q, r, transpose)
     514
     515\param q
     516See \c RevSparseJac(q, r, transpose)
     517
     518\param r
     519See \c RevSparseJac(q, r, transpose)
    434520
    435521\param s
    436 See \c RevSparseJac(p, s).
    437 
    438 \param r
    439 is the return value for the corresponding call to RevSparseJac(p, s);
     522is the return value for the corresponding call to RevSparseJac(q, r, transpose)
    440523*/
    441524
     
    444527void ADFun<Base>::RevSparseJacCase(
    445528        const std::set<size_t>&      set_type          ,
    446         size_t                       p                 ,
    447         const VectorSet&             s                 ,
    448         VectorSet&                   r                 )
     529        bool                         transpose         ,
     530        size_t                       q                 ,
     531        const VectorSet&             r                 ,
     532        VectorSet&                   s                 )
    449533{       // dimension of the result vector
    450         r.resize( p );
     534        if( transpose )
     535                s.resize( Domain() );
     536        else    s.resize( q );
    451537
    452538        // store results in r
    453539        RevSparseJacSet(
    454                 p              ,
     540                transpose      ,
     541                q              ,
     542                r              ,
    455543                s              ,
    456                 r              ,
    457544                total_num_var_ ,
    458545                dep_taddr_     ,
     
    467554The C++ source code corresponding to this operation is
    468555\verbatim
    469         s = f.RevSparseJac(q, r)
     556        s = f.RevSparseJac(q, r, transpose)
    470557\endverbatim
    471558
     
    477564or \c std::set<size_t>.
    478565
    479 \param p
    480 is the number of rows in the matrix \f$ S \f$.
    481 
    482 \param s
    483 is a sparsity pattern for the matrix \f$ S \f$.
     566\param q
     567is the number of rows in the matrix \f$ R \f$.
     568
     569\param r
     570is a sparsity pattern for the matrix \f$ R \f$.
     571
     572\param transpose
     573are the sparsity patterns for \f$ R \f$ and \f$ S(x) \f$ transposed.
    484574
    485575\return
     576If \c transpose is false (true), the return value is a sparsity pattern
     577for \f$ S(x) \f$ (\f$ S(x)^T \f$) where
     578\f[
     579        S(x) = R * F^{(1)} (x)
     580\f]
     581and \f$ F \f$ is the function corresponding to the operation sequence
     582and \a x is any argument value.
    486583If \c VectorSet::value_type is \c bool,
    487 the return value \c r is a vector with size \c p*n
    488 where \c n is the number of independent variables
    489 corresponding to the operation sequence stored in \c f.
     584the return value has size \f$ q * n \f$ ( \f$ n * q \f$).
    490585If \c VectorSet::value_type is \c std::set<size_t>,
    491 the