Changeset 2401


Ignore:
Timestamp:
May 24, 2012 12:30:25 PM (8 years ago)
Author:
bradbell
Message:

Merge in changes from branches/sparse.

  1. Update version number to today.
  2. Add sparse return user API interfaces to SparseJacobian? and SparseHessian?.
  3. Include new interface and imporve sparse Hessian and Jacobian examples.
  4. More testing of sparse Hessian and Jacobians.
  5. Add is_element to sparse_set and sparse_pack.

.: merge information.
search.sh: search cppad/speed directory, drop makefile.in from results.
svn_merge.sh: parameters for this merge.
parallel_ad.hpp: initilize is_element statics.
sparse_evaluate.hpp: fix typo in documentation.
glossary.omh: convert to modern omhelp, improve AD Function entry.
link_sparse_hessian.cpp: fix typo in documentation.

Location:
trunk
Files:
23 edited

Legend:

Unmodified
Added
Removed
  • trunk

  • trunk/AUTHORS

    r2355 r2401  
    22             ===========================================
    33
    4 To date, 2012-04-23, Bradley M. Bell is the sole author of CppAD.
     4To date, 2012-05-24, 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/bin/commit.sh

    r2239 r2401  
    1313# replacement text for this commit
    1414cat << EOF > bin/commit.user.$$
    15 This is a template file for making commits to the CppAD repository.
    16 Lines with no 'at' characters, are general comments not connected to
    17 a specific file. Lines containing an 'at' character are "file name"
    18 followed by comment. Lines before the first 'at' character are preserved
    19 during
    20         bin/commit.sh edit
    21 for example this entire paragraph is preserved.
    22 
    23 dir/file.ext@ optional comment about this file.
     15Merge in changes from branches/sparse.
     161. Update version number to today.
     172. Add sparse return user API interfaces to SparseJacobian and SparseHessian.
     183. Include new interface and imporve sparse Hessian and Jacobian examples.
     194. More testing of sparse Hessian and Jacobians.
     205. Add is_element to sparse_set and sparse_pack.
     21
     22.@ merge information.
     23AUTHORS@
     24bin/search.sh@ search cppad/speed directory, drop makefile.in from results.
     25bin/svn_merge.sh@ parameters for this merge.
     26configure@
     27configure.ac@
     28cppad/configure.hpp@
     29cppad/local/ad_fun.hpp@
     30cppad/local/declare_ad.hpp@
     31cppad/local/parallel_ad.hpp@ initilize is_element statics.
     32cppad/local/sparse_hessian.hpp@
     33cppad/local/sparse_jacobian.hpp@
     34cppad/local/sparse_pack.hpp@
     35cppad/local/sparse_set.hpp@
     36cppad/speed/sparse_evaluate.hpp@ fix typo in documentation.
     37example/sparse_hessian.cpp@
     38example/sparse_jacobian.cpp@
     39omh/glossary.omh@ convert to modern omhelp, improve AD Function entry.
     40omh/whats_new_12.omh@
     41speed/src/link_sparse_hessian.cpp@ fix typo in documentation.
     42test_more/sparse_hessian.cpp@
     43test_more/sparse_jacobian.cpp@
    2444EOF
    2545# -----------------------------------------------------------------------------
     
    106126# -------------------------------------------------
    107127# list of files that changed
    108 svn status | sed -n -e '/^[ADMRC] /p' | \
    109         sed -e 's/^[ADMRC] [+ ]*//' \
     128svn status | sed -n -e '/^[ADMRC][ADMRC]* /p' | \
     129        sed -e 's/^[ADMRC][ADMRC]* [+ ]*//' \
    110130                -e '/^bin\/commit.sh$/d' -e '/^bin\/commit.sed$/d' | \
    111131        sort -u > bin/commit.list.$$
  • trunk/bin/search.sh

    • Property svn:keywords set to Id
    r2184 r2401  
    11#! /bin/bash -e
    2 # $Id: search.sh 2082 2011-08-31 17:50:58Z bradbell $
     2# $Id$
    33# -----------------------------------------------------------------------------
    4 # CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-11 Bradley M. Bell
     4# CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
    55#
    66# CppAD is distributed under multiple licenses. This distribution is under
     
    2727        cppad
    2828        cppad/local
     29        cppad/speed
    2930        cppad_ipopt/example
    3031        cppad_ipopt/speed
     
    5354'
    5455#
    55 find_files.sh "$pattern" "$extensions" "$directories"
     56find_files.sh "$pattern" "$extensions" "$directories" | \
     57         sed -e '/\/makefile.in/d'
  • trunk/bin/svn_merge.sh

    r2178 r2401  
    22# $Id$
    33# -----------------------------------------------------------------------------
    4 # CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-11 Bradley M. Bell
     4# CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
    55#
    66# CppAD is distributed under multiple licenses. This distribution is under
     
    3030#
    3131# Name of the directory where the changes have been committed
    32 from_branch=branches/thread
     32from_branch=branches/sparse
    3333#
    3434# Version of the repository corresponding to from_branch just before changes
    35 Start=2124
     35Start=2357
    3636#
    3737# Version of the repository corresponding to from_branch after the changes
    38 End=2177
     38End=2400
    3939#
    4040# the svn merge command
  • trunk/configure

    r2355 r2401  
    11#! /bin/sh
    22# Guess values for system-dependent variables and create Makefiles.
    3 # Generated by GNU Autoconf 2.68 for CppAD 20120423.
     3# Generated by GNU Autoconf 2.68 for CppAD 20120524.
    44#
    55# Report bugs to <cppad@list.coin-or.org>.
     
    561561PACKAGE_NAME='CppAD'
    562562PACKAGE_TARNAME='cppad'
    563 PACKAGE_VERSION='20120423'
    564 PACKAGE_STRING='CppAD 20120423'
     563PACKAGE_VERSION='20120524'
     564PACKAGE_STRING='CppAD 20120524'
    565565PACKAGE_BUGREPORT='cppad@list.coin-or.org'
    566566PACKAGE_URL=''
     
    13381338  # This message is too long to be a string in the A/UX 3.1 sh.
    13391339  cat <<_ACEOF
    1340 \`configure' configures CppAD 20120423 to adapt to many kinds of systems.
     1340\`configure' configures CppAD 20120524 to adapt to many kinds of systems.
    13411341
    13421342Usage: $0 [OPTION]... [VAR=VALUE]...
     
    14041404if test -n "$ac_init_help"; then
    14051405  case $ac_init_help in
    1406      short | recursive ) echo "Configuration of CppAD 20120423:";;
     1406     short | recursive ) echo "Configuration of CppAD 20120524:";;
    14071407   esac
    14081408  cat <<\_ACEOF
     
    15241524if $ac_init_version; then
    15251525  cat <<\_ACEOF
    1526 CppAD configure 20120423
     1526CppAD configure 20120524
    15271527generated by GNU Autoconf 2.68
    15281528
     
    21482148running configure, to aid debugging if configure makes a mistake.
    21492149
    2150 It was created by CppAD $as_me 20120423, which was
     2150It was created by CppAD $as_me 20120524, which was
    21512151generated by GNU Autoconf 2.68.  Invocation command line was
    21522152
     
    48304830# Define the identity of the package.
    48314831 PACKAGE='cppad'
    4832  VERSION='20120423'
     4832 VERSION='20120524'
    48334833
    48344834
     
    72887288# values after options handling.
    72897289ac_log="
    7290 This file was extended by CppAD $as_me 20120423, which was
     7290This file was extended by CppAD $as_me 20120524, which was
    72917291generated by GNU Autoconf 2.68.  Invocation command line was
    72927292
     
    73457345ac_cs_config="`$as_echo "$ac_configure_args" | sed 's/^ //; s/[\\""\`\$]/\\\\&/g'`"
    73467346ac_cs_version="\\
    7347 CppAD config.status 20120423
     7347CppAD config.status 20120524
    73487348configured by $0, generated by GNU Autoconf 2.68,
    73497349  with options \\"\$ac_cs_config\\"
  • trunk/configure.ac

    r2355 r2401  
    1313dnl Process this file with autoconf to produce a configure script.
    1414dnl   package   version              bug-report
    15 AC_INIT(CppAD, 20120423, cppad@list.coin-or.org)
     15AC_INIT(CppAD, 20120524, cppad@list.coin-or.org)
    1616AC_SUBST(PACKAGE_URL, "http://www.coin-or.org/CppAD")
    1717AC_SUBST(PACKAGE_DESCRIPTION, "Differentiation of C++ Algorithms")
  • trunk/cppad/configure.hpp

    r2355 r2401  
    4343cppad-yyyymmdd as a C string where yyyy is year, mm is month, and dd is day.
    4444*/
    45 # define CPPAD_PACKAGE_STRING "cppad-20120423"
     45# define CPPAD_PACKAGE_STRING "cppad-20120524"
    4646
    4747/*!
  • trunk/cppad/local/ad_fun.hpp

    r2178 r2401  
    44
    55/* --------------------------------------------------------------------------
    6 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-11 Bradley M. Bell
     6CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
    77
    88CppAD is distributed under multiple licenses. This distribution is under
     
    178178        );
    179179        // ------------------------------------------------------------
     180        // Forward mode version of SparseJacobian
     181        // (see doxygen in sparse_jacobian.hpp)
     182        template <class VectorBase, class VectorSet>
     183        size_t SparseJacobianFor(
     184                const VectorBase&     x               ,
     185                VectorSet&            p_transpose     ,
     186                VectorBase&           jac             ,
     187                sparse_jacobian_work& work
     188        );
     189        // Reverse mode version of SparseJacobian
     190        // (see doxygen in sparse_jacobian.hpp)
     191        template <class VectorBase, class VectorSet>
     192        size_t SparseJacobianRev(
     193                const VectorBase&     x               ,
     194                VectorSet&            p               ,
     195                VectorBase&           jac             ,
     196                sparse_jacobian_work& work
     197        );
     198        // ------------------------------------------------------------
    180199        // vector of bool version of SparseJacobian
    181200        // (see doxygen in sparse_jacobian.hpp)
    182201        template <class VectorBase, class VectorSet>
     202        size_t SparseJacobianCase(
     203                bool                     set_type    ,
     204                const VectorBase&        x           ,
     205                const VectorSet&         p           ,
     206                VectorBase&              jac         ,
     207                sparse_jacobian_work&    work
     208        );
     209        // vector of std::set<size_t> version of SparseJacobian
     210        // (see doxygen in sparse_jacobian.hpp)
     211        template <class VectorBase, class VectorSet>
     212        size_t SparseJacobianCase(
     213                const std::set<size_t>&  set_type    ,
     214                const VectorBase&        x           ,
     215                const VectorSet&         p           ,
     216                VectorBase&              jac         ,
     217                sparse_jacobian_work&    work
     218        );
     219        // vector of bool version of SparseJacobian
     220        // (see doxygen in sparse_jacobian.hpp)
     221        template <class VectorBase, class VectorSet>
    183222        void SparseJacobianCase(
    184223                bool                     set_type    ,
     
    197236        );
    198237        // ------------------------------------------------------------
     238        // combined sparse_set and sparse_pack version of SparseHessian
     239        // (see doxygen in sparse_hessian.hpp)
     240        template <class VectorBase, class VectorSet>
     241        size_t SparseHessianCompute(
     242                const VectorBase&        x           ,
     243                const VectorBase&        w           ,
     244                VectorSet&               sparsity    ,
     245                VectorBase&              hes         ,
     246                sparse_hessian_work&     work
     247        );
     248        // vector of bool version of SparseHessian
     249        // (see doxygen in sparse_hessian.hpp)
     250        template <class VectorBase, class VectorSet>
     251        size_t SparseHessianCase(
     252                bool                     set_type    ,
     253                const VectorBase&        x           ,
     254                const VectorBase&        w           ,
     255                const VectorSet&         p           ,
     256                VectorBase&              hes         ,
     257                sparse_hessian_work&     work
     258        );
     259        // vector of std::set<size_t> version of SparseHessian
     260        // (see doxygen in sparse_hessian.hpp)
     261        template <class VectorBase, class VectorSet>
     262        size_t SparseHessianCase(
     263                const std::set<size_t>&  set_type    ,
     264                const VectorBase&        x           ,
     265                const VectorBase&        w           ,
     266                const VectorSet&         p           ,
     267                VectorBase&              hes         ,
     268                sparse_hessian_work&     work
     269        );
    199270        // vector of bool version of SparseHessian
    200271        // (see doxygen in sparse_hessian.hpp)
     
    406477        /// calculate sparse Jacobians
    407478        template <typename VectorBase>
    408         VectorBase SparseJacobian(const VectorBase &x);
    409 
    410         /// calculate sparse Jacobians
     479        VectorBase SparseJacobian(
     480                const VectorBase &x
     481        );
    411482        template <typename VectorBase, typename VectorSet>
    412         VectorBase SparseJacobian(const VectorBase &x, const VectorSet &p);
     483        VectorBase SparseJacobian(
     484                const VectorBase &x ,
     485                const VectorSet  &p
     486        );
     487        template <class VectorBase, class VectorSet, class VectorSize>
     488        size_t SparseJacobianForward(
     489                const VectorBase&     x     ,
     490                const VectorSet&      p     ,
     491                const VectorSize&     r     ,
     492                const VectorSize&     c     ,
     493                VectorBase&           jac   ,
     494                sparse_jacobian_work& work
     495        );
     496        template <class VectorBase, class VectorSet, class VectorSize>
     497        size_t SparseJacobianReverse(
     498                const VectorBase&     x    ,
     499                const VectorSet&      p    ,
     500                const VectorSize&     r    ,
     501                const VectorSize&     c    ,
     502                VectorBase&           jac  ,
     503                sparse_jacobian_work& work
     504        );
    413505
    414506        /// calculate sparse Hessians
    415507        template <typename VectorBase>
    416         VectorBase SparseHessian(const VectorBase &x, const VectorBase &w);
    417 
    418         /// calculate sparse Hessians
     508        VectorBase SparseHessian(
     509                const VectorBase&    x  ,
     510                const VectorBase&    w
     511        );
    419512        template <typename VectorBase, typename VectorBool>
    420513        VectorBase SparseHessian(
    421                 const VectorBase &x, const VectorBase &w, const VectorBool &p
     514                const VectorBase&    x  ,
     515                const VectorBase&    w  ,
     516                const VectorBool&    p
     517        );
     518        template <class VectorBase, class VectorSet, class VectorSize>
     519        size_t SparseHessian(
     520                const VectorBase&    x   ,
     521                const VectorBase&    w   ,
     522                const VectorSet&     p   ,
     523                const VectorSize&    r   ,
     524                const VectorSize&    c   ,
     525                VectorBase&          hes ,
     526                sparse_hessian_work& work
    422527        );
    423528
  • trunk/cppad/local/declare_ad.hpp

    r2038 r2401  
    44
    55/* --------------------------------------------------------------------------
    6 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-11 Bradley M. Bell
     6CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
    77
    88CppAD is distributed under multiple licenses. This distribution is under
     
    1919
    2020        // classes
     21        class sparse_jacobian_work;
     22        class sparse_hessian_work;
    2123        template <class Base> class AD;
    2224        template <class Base> class ADFun;
  • trunk/cppad/local/parallel_ad.hpp

    r2344 r2401  
    111111        sp.next_element();     // has static data
    112112        sp.clear(0);           // has static data
    113 
     113        sp.is_element(0, 0);   // has static data
    114114
    115115        // statics that depend on the value of Base
  • trunk/cppad/local/sparse_hessian.hpp

    r2006 r2401  
    44
    55/* --------------------------------------------------------------------------
    6 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-11 Bradley M. Bell
     6CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
    77
    88CppAD is distributed under multiple licenses. This distribution is under
     
    1717$begin sparse_hessian$$
    1818$spell
     19        recomputed
    1920        CppAD
    2021        valarray
     
    3132
    3233$head Syntax$$
    33 $codei%%hes% = %f%.SparseHessian(%x%, %w%)
     34$icode%hes% = %f%.SparseHessian(%x%, %w%)
     35%hes% = %f%.SparseHessian(%x%, %w%, %p%)
     36%n_sweep% = %f%.SparseHessian(%x%, %w%, %p%, %row%, %col%, %hes%, %work%)
    3437%$$
    35 $codei%%hes% = %f%.SparseHessian(%x%, %w%, %p%)%$$
    3638
    3739$head Purpose$$
     40We use $latex n$$ for the $cref/domain/seq_property/Domain/$$ size,
     41and $latex m$$ for the $cref/range/seq_property/Range/$$ size of $icode f$$.
    3842We use $latex F : \R^n \rightarrow \R^m$$ do denote the
    3943$cref/AD function/glossary/AD Function/$$
     
    4145The syntax above sets $icode hes$$ to the Hessian
    4246$latex \[
    43         hes = \dpow{2}{x} \sum_{i=1}^m w_i F_i (x)
     47        H(x) = \dpow{2}{x} \sum_{i=1}^m w_i F_i (x)
    4448\] $$
    45 This routine assumes  that the matrix $icode hes$$ is sparse
    46 and uses this assumption sparse to reduce the amount of computation necessary.
    47 One should use speed tests to verify that results are computed faster
    48 than when using the routine $cref/Hessian/$$.
     49This routine takes advantage of the sparsity of the Hessian
     50in order to reduce the amount of computation necessary.
     51If $icode row$$ and $icode col$$ are present, it also takes
     52advantage of the reduced set of elements of the Hessian that
     53need to be computed.
     54One can use speed tests (e.g. $cref speed_test$$)
     55to verify that results are computed faster
     56than when using the routine $cref Hessian$$.
    4957
    5058$head f$$
     
    5361        ADFun<%Base%> %f%
    5462%$$
    55 Note that the $cref/ADFun/$$ object $icode f$$ is not $code const$$
     63Note that the $cref ADFun$$ object $icode f$$ is not $code const$$
    5664(see $cref/Uses Forward/sparse_hessian/Uses Forward/$$ below).
    5765
     
    93101It specifies a
    94102$cref/sparsity pattern/glossary/Sparsity Pattern/$$
    95 for the Hessian $icode hes$$.
     103for the Hessian $latex H(x)$$.
    96104$pre
    97105
     
    108116for the internal calculations is unspecified.
    109117
     118$head row, col$$
     119The arguments $icode row$$ and $icode col$$ are optional and have prototype
     120$codei%
     121        const %VectorSize%& %row%
     122        const %VectorSize%& %col%
     123%$$
     124(see $cref/VectorSize/sparse_hessian/VectorSize/$$ below).
     125They specify which rows and columns of $latex H (x)$$ are
     126returned and in what order.
     127We use $latex K$$ to denote the value $icode%hes%.size()%$$
     128which must also equal the size of $icode row$$ and $icode col$$.
     129Furthermore,
     130for $latex k = 0 , \ldots , K-1$$, it must hold that
     131$latex row[k] < n$$ and $latex col[k] < n$$.
     132In addition,
     133all of the $latex (row[k], col[k])$$ pairs must correspond to a true value
     134in the sparsity pattern $icode p$$.
     135
    110136$head hes$$
    111137The result $icode hes$$ has prototype
     
    113139        %VectorBase% %hes%
    114140%$$
    115 and its size is $latex n * n$$.
    116 For $latex j = 0 , \ldots , n - 1 $$
    117 and $latex \ell = 0 , \ldots , n - 1$$
     141In the case where $icode row$$ and $icode col$$ are not present,
     142the size of $icode hes$$ is $latex n * n$$ and
     143its size is $latex n * n$$.
     144In this case, for $latex i = 0 , \ldots , n - 1 $$
     145and $latex ell = 0 , \ldots , n - 1$$
    118146$latex \[
    119147        hes [ j * n + \ell ] = \DD{ w^{\rm T} F }{ x_j }{ x_\ell } ( x )
    120148\] $$
     149$pre
     150
     151$$
     152In the case where the arguments $icode row$$ and $icode col$$ are present,
     153we use $latex K$$ to denote the size of $icode hes$$.
     154The input value of its elements does not matter.
     155Upon return, for $latex k = 0 , \ldots , K - 1$$,
     156$latex \[
     157        hes [ k ] = \DD{ w^{\rm T} F }{ x_j }{ x_\ell } (x)
     158        \; , \;
     159        \; {\rm where} \;
     160        j = row[k]
     161        \; {\rm and } \;
     162        \ell = col[k]
     163\] $$
     164
     165$head work$$
     166If this argument is present, it has prototype
     167$icode%
     168        sparse_hessian_work& %work%
     169%$$
     170This object can only be used with the routines $code SparseHessian$$.
     171During its the first use, information is stored in $icode work$$.
     172This is used to reduce the work done by future calls to $code SparseHessian$$
     173with the same $icode f$$, $icode p$$, $icode row$$, and $icode col$$.
     174If a future call is make where any of these values have changed,
     175you must first call $icode%work%.clear()%$$
     176to inform CppAD that this information needs to be recomputed.
     177
     178$head n_sweep$$
     179The return value $icode n_sweep$$ has prototype
     180$codei%
     181        size_t %n_sweep%
     182%$$
     183It is the number of first order forward sweeps
     184used to compute the requested Hessian values.
     185Each first forward sweep is followed by a second order reverse sweep
     186so it is also the number of reverse sweeps.
     187This is proportional to the total work that $code SparseHessian$$ does,
     188not counting the zero order forward sweep,
     189or the work to combine multiple columns into a single
     190forward-reverse sweep pair.
    121191
    122192$head VectorBase$$
    123 The type $icode VectorBase$$ must be a $cref/SimpleVector/$$ class with
     193The type $icode VectorBase$$ must be a $cref SimpleVector$$ class with
    124194$cref/elements of type/SimpleVector/Elements of Specified Type/$$
    125195$icode Base$$.
    126 The routine $cref/CheckSimpleVector/$$ will generate an error message
     196The routine $cref CheckSimpleVector$$ will generate an error message
    127197if this is not the case.
    128198
    129199$head VectorSet$$
    130 The type $icode VectorSet$$ must be a $xref/SimpleVector/$$ class with
    131 $xref/SimpleVector/Elements of Specified Type/elements of type/$$
     200The type $icode VectorSet$$ must be a $cref SimpleVector$$ class with
     201$cref/elements of type/SimpleVector/Elements of Specified Type/$$
    132202$code bool$$ or $code std::set<size_t>$$;
    133203see $cref/sparsity pattern/glossary/Sparsity Pattern/$$ for a discussion
    134204of the difference.
    135 The routine $xref/CheckSimpleVector/$$ will generate an error message
     205The routine $cref CheckSimpleVector$$ will generate an error message
    136206if this is not the case.
    137207
     
    144214this condition.
    145215
     216$head VectorSize$$
     217The type $icode VectorSize$$ must be a $cref SimpleVector$$ class with
     218$cref/elements of type/SimpleVector/Elements of Specified Type/$$
     219$code size_t$$.
     220The routine $cref CheckSimpleVector$$ will generate an error message
     221if this is not the case.
     222
    146223$head Uses Forward$$
    147 After each call to $cref/Forward/$$,
     224After each call to $cref Forward$$,
    148225the object $icode f$$ contains the corresponding
    149226$cref/Taylor coefficients/glossary/Taylor Coefficient/$$.
    150227After $code SparseHessian$$,
    151 the previous calls to $xref/Forward/$$ are undefined.
     228the previous calls to $cref Forward$$ are undefined.
    152229
    153230$head Example$$
     
    156233%$$
    157234The routine
    158 $cref/sparse_hessian.cpp/$$
     235$cref sparse_hessian.cpp$$
    159236is examples and tests of $code sparse_hessian$$.
    160237It return $code true$$, if it succeeds and $code false$$ otherwise.
     
    171248Sparse Hessian driver routine and helper functions.
    172249*/
    173 
     250// ===========================================================================
     251/*!
     252class used by SparseHessian to hold information
     253so it does not need to be recomputed.
     254*/
     255class sparse_hessian_work {
     256        public:
     257                /// version of user r array sorted by row or column
     258                CppAD::vector<size_t> r_sort;
     259                /// version of user c array sorted by row or column
     260                CppAD::vector<size_t> c_sort;
     261                /// mapping from sorted array indices to user array indices
     262                CppAD::vector<size_t> k_sort;
     263                /// results of the coloring algorithm
     264                CppAD::vector<size_t> color;
     265                /// inform CppAD that this information needs to be recomputed
     266                void clear(void)
     267                {       r_sort.resize(0);
     268                        c_sort.resize(0);
     269                        k_sort.resize(0);
     270                        color.resize(0);
     271                }
     272};
     273// ===========================================================================
     274/*!
     275Private helper function that does computation for all Sparse Hessian cases.
     276
     277\tparam Base
     278See \c SparseHessian(x, w, p, row, col, hes, work).
     279
     280\tparam VectorBase
     281See \c SparseHessian(x, w, p, row, col, hes, work).
     282
     283\tparam VectorSet
     284is either \c sparse_pack or \c sparse_set.
     285
     286\param x
     287See \c SparseHessian(x, w, p, row, col, hes, work).
     288
     289\param w
     290See \c SparseHessian(x, w, p, row, col, hes, work).
     291
     292\param sparsity
     293If <code>work.color.size() != 0</code>, then \c sparsity is not used.
     294Otherwise, it is a
     295sparsity pattern for the Hessian of this ADFun<Base> object.
     296
     297\param hes
     298See \c SparseHessian(x, w, p, row, col, hes, work).
     299
     300\param work
     301See \c SparseHessian(x, w, p, row, col, hes, work).
     302
     303\return
     304See \c SparseHessian(x, w, p, row, col, hes, work).
     305*/
     306template<class Base>
     307template <class VectorBase, class VectorSet>
     308size_t ADFun<Base>::SparseHessianCompute(
     309        const VectorBase&     x           ,
     310        const VectorBase&     w           ,
     311        VectorSet&            sparsity    ,
     312        VectorBase&           hes         ,
     313        sparse_hessian_work& work         )
     314{
     315        using   CppAD::vectorBool;
     316        size_t i, j, k, ell;
     317
     318        CppAD::vector<size_t>& row(work.r_sort);
     319        CppAD::vector<size_t>& col(work.c_sort);
     320        CppAD::vector<size_t>& user_k(work.k_sort);
     321        CppAD::vector<size_t>& color(work.color);
     322
     323        size_t n = Domain();
     324
     325        // some values
     326        const Base zero(0);
     327        const Base one(1);
     328
     329        // check VectorBase is Simple Vector class with Base type elements
     330        CheckSimpleVector<Base, VectorBase>();
     331
     332        CPPAD_ASSERT_UNKNOWN( x.size() == n );
     333        CPPAD_ASSERT_UNKNOWN( color.size() == 0 || color.size() == n );
     334
     335        // number of components of Hessian that are required
     336        size_t K = hes.size();
     337        CPPAD_ASSERT_UNKNOWN( row.size() == K+1 );
     338        CPPAD_ASSERT_UNKNOWN( col.size() == K );
     339        CPPAD_ASSERT_UNKNOWN( row[K] == n );
     340
     341        // Point at which we are evaluating the Hessian
     342        Forward(0, x);
     343
     344        // Rows of the Hessian (i below) correspond to the forward mode index
     345        // and columns (j below) correspond to the reverse mode index.
     346        if( color.size() == 0 )
     347        {
     348                CPPAD_ASSERT_UNKNOWN( sparsity.n_set() ==  n );
     349                CPPAD_ASSERT_UNKNOWN( sparsity.end() ==  n );
     350
     351                // rows and columns that are in the returned hessian
     352                VectorSet r_used, c_used;
     353                r_used.resize(n, n);
     354                c_used.resize(n, n);
     355                k = 0;
     356                while( k < K )
     357                {       CPPAD_ASSERT_UNKNOWN( row[k] < n && col[k] < n );
     358                        CPPAD_ASSERT_UNKNOWN( k == 0 || row[k-1] <= row[k] );
     359                        CPPAD_ASSERT_KNOWN(
     360                                sparsity.is_element(row[k], col[k]) ,
     361                                "SparseHessian: an (row, col) pair is not in sparsity pattern."
     362                        );
     363                        r_used.add_element(col[k], row[k]);
     364                        c_used.add_element(row[k], col[k]);
     365                        k++;
     366                }
     367       
     368                // given a column index, which rows are non-zero and not used
     369                VectorSet not_used;
     370                not_used.resize(n, n);
     371                for(i = 0; i < n; i++)
     372                {       sparsity.begin(i);
     373                        j = sparsity.next_element();
     374                        while( j != sparsity.end() )
     375                        {       if( ! r_used.is_element(j , i) )
     376                                        not_used.add_element(j, i);
     377                                j = sparsity.next_element();
     378                        }
     379                }
     380
     381                // initial coloring
     382                color.resize(n);
     383                for(i = 0; i < n; i++)
     384                        color[i] = i;
     385       
     386                // See GreedyPartialD2Coloring Algorithm Section 3.6.2 of
     387                // Graph Coloring in Optimization Revisited by
     388                // Assefaw Gebremedhin, Fredrik Maane, Alex Pothen
     389                vectorBool forbidden(n);
     390                for(i = 1; i < n; i++)
     391                {
     392                        // initial all colors as ok for this row
     393                        // (value of forbidden for ell > i does not matter)
     394                        for(ell = 0; ell <= i; ell++)
     395                                forbidden[ell] = false;
     396       
     397                        // -----------------------------------------------------
     398                        // Forbid colors that this row would destroy results for.
     399                        // for each column that is non-zero for this row
     400                        sparsity.begin(i);
     401                        j = sparsity.next_element();
     402                        while( j != sparsity.end() )
     403                        {       // for each row that this column uses
     404                                r_used.begin(j);
     405                                ell = r_used.next_element();
     406                                while( ell != r_used.end() )
     407                                {       // if this is not the same row, forbid its color
     408                                        if( ell < i )
     409                                                forbidden[ color[ell] ] = true;
     410                                        ell = r_used.next_element();
     411                                }
     412                                j = sparsity.next_element();
     413                        }
     414       
     415                        // -------------------------------------------------------
     416                        // Forbid colors that would destroy the results for this row.
     417                        // for each column that this row used
     418                        c_used.begin(i);
     419                        j = c_used.next_element();
     420                        while( j != c_used.end() )
     421                        {       // For each row that is non-zero for this column
     422                                // (the used rows have already been checked above).
     423                                not_used.begin(j);
     424                                ell = not_used.next_element();
     425                                while( ell != not_used.end() )
     426                                {       // if this is not the same row, forbid its color
     427                                        if( ell < i )
     428                                                forbidden[ color[ell] ] = true;
     429                                        ell = not_used.next_element();
     430                                }
     431                                j = c_used.next_element();
     432                        }
     433
     434                        // pick the color with the smallest index
     435                        ell = 0;
     436                        while( forbidden[ell] )
     437                        {       ell++;
     438                                CPPAD_ASSERT_UNKNOWN( ell <= i );
     439                        }
     440                        color[i] = ell;
     441                }
     442        }
     443        size_t n_color = 1;
     444        for(ell = 0; ell < n; ell++)
     445                n_color = std::max(n_color, color[ell] + 1);
     446
     447        // direction vector for calls to forward (rows of the Hessian)
     448        VectorBase u(n);
     449
     450        // location for return values from reverse (columns of the Hessian)
     451        VectorBase ddw(2 * n);
     452
     453        // initialize the return value
     454        for(k = 0; k < K; k++)
     455                hes[k] = zero;
     456
     457        // loop over colors
     458        size_t n_sweep = 0;
     459        for(ell = 0; ell < n_color; ell++)
     460        {       bool any = false;
     461                k = 0;
     462                for(i = 0; i < n; i++) if( color[i] == ell )
     463                {       // find first k such that row[k] has color ell
     464                        if( ! any )
     465                        {       while( row[k] < i )
     466                                        k++;
     467                                any = row[k] == i;
     468                        }
     469                }
     470                if( any )
     471                {       n_sweep++;
     472                        // combine all rows with this color
     473                        for(i = 0; i < n; i++)
     474                        {       u[i] = zero;
     475                                if( color[i] == ell )
     476                                        u[i] = one;
     477                        }
     478                        // call forward mode for all these rows at once
     479                        Forward(1, u);
     480
     481                        // evaluate derivative of w^T * F'(x) * u
     482                        ddw = Reverse(2, w);
     483
     484                        // set the corresponding components of the result
     485                        for(i = 0; i < n; i++) if( color[i] == ell )
     486                        {       // find first index in c for this column
     487                                while( row[k] < i )
     488                                        k++;
     489                                // extract the results for this row
     490                                while( row[k] == i )
     491                                {       hes[ user_k[k] ] = ddw[ col[k] * 2 + 1 ];
     492                                        k++;
     493                                }
     494                        }
     495                }
     496        }
     497        return n_sweep;
     498}
     499// ===========================================================================
     500/*!
     501Private helper function for vector of \c bool sparsity pattern cases.
     502
     503\tparam Base
     504See \c SparseHessian(x, w, p, row, col, hes, work).
     505
     506\tparam VectorBase
     507See \c SparseHessian(x, w, p, row, col, hes, work).
     508
     509\tparam VectorSet
     510is a simple vector with elements of type \c bool.
     511
     512\param set_type
     513has element type for vector representing the sparsity sets.
     514
     515\param x
     516See \c SparseHessian(x, w, p, row, col, hes, work).
     517
     518\param w
     519See \c SparseHessian(x, w, p, row, col, hes, work).
     520
     521\param p
     522Sparsity pattern for the Hessian of this ADFun<Base> object.
     523
     524\param hes
     525See \c SparseHessian(x, w, p, row, col, hes, work).
     526
     527\param work
     528See \c SparseHessian(x, w, p, row, col, hes, work).
     529
     530\return
     531See \c SparseHessian(x, w, p, row, col, hes, work).
     532*/
     533template <class Base>
     534template <class VectorBase, class VectorSet>
     535size_t ADFun<Base>::SparseHessianCase(
     536        bool                  set_type  ,
     537        const VectorBase&     x         ,
     538        const VectorBase&     w         ,
     539        const VectorSet&      p         ,
     540        VectorBase&           hes       ,
     541        sparse_hessian_work&  work      )
     542{
     543        size_t n = Domain();
     544
     545        // check VectorSet is Simple Vector class with bool elements
     546        CheckSimpleVector<bool, VectorSet>();
     547
     548        // check VectorBase is Simple Vector class with Base type elements
     549        CheckSimpleVector<Base, VectorBase>();
     550
     551        CPPAD_ASSERT_UNKNOWN( x.size() == n );
     552        CPPAD_ASSERT_UNKNOWN( w.size() == Range() );
     553        CPPAD_ASSERT_KNOWN(
     554                p.size() == n * n,
     555                "SparseHessian: using bool values for sparsity and p.size() "
     556                "not equal square of domain dimension for f"
     557        );
     558 
     559        sparse_pack sparsity;
     560        if( work.color.size() == 0 )
     561        {       bool transpose = false;
     562                vec_bool_to_sparse_pack(sparsity, p, n, n, transpose);
     563        }
     564       
     565        // compute the Hessian
     566        size_t n_sweep = SparseHessianCompute(x, w, sparsity, hes, work);
     567
     568        return n_sweep;
     569}
     570/*!
     571Private helper function for vector of std::set<size_t> sparsity pattern cases.
     572
     573\tparam Base
     574See \c SparseHessian(x, w, p, row, col, hes, work).
     575
     576\tparam VectorBase
     577See \c SparseHessian(x, w, p, row, col, hes, work).
     578
     579\tparam VectorSet
     580is a simple vector with elements of type <code>std::set<size_t></code>.
     581
     582\param set_type
     583has element type for vector representing the sparsity sets.
     584
     585\param x
     586See \c SparseHessian(x, w, p, row, col, hes, work).
     587
     588\param w
     589See \c SparseHessian(x, w, p, row, col, hes, work).
     590
     591\param p
     592Sparsity pattern for the Hessian of this ADFun<Base> object.
     593
     594\param hes
     595See \c SparseHessian(x, w, p, row, col, hes, work).
     596
     597\param work
     598See \c SparseHessian(x, w, p, row, col, hes, work).
     599
     600\return
     601See \c SparseHessian(x, w, p, row, col, hes, work).
     602*/
     603template <class Base>
     604template <class VectorBase, class VectorSet>
     605size_t ADFun<Base>::SparseHessianCase(
     606        const std::set<size_t>&     set_type  ,
     607        const VectorBase&           x         ,
     608        const VectorBase&           w         ,
     609        const VectorSet&            p         ,
     610        VectorBase&                 hes       ,
     611        sparse_hessian_work&        work      )
     612{
     613        size_t n = Domain();
     614
     615        // check VectorSet is Simple Vector class with std::set<size_t> elements
     616        CheckSimpleVector<std::set<size_t>, VectorSet>(
     617                one_element_std_set<size_t>(), two_element_std_set<size_t>()
     618        );
     619
     620        // check VectorBase is Simple Vector class with Base type elements
     621        CheckSimpleVector<Base, VectorBase>();
     622
     623        CPPAD_ASSERT_UNKNOWN( x.size() == n );
     624        CPPAD_ASSERT_UNKNOWN( w.size() == Range() );
     625        CPPAD_ASSERT_KNOWN(
     626                p.size() == n,
     627                "SparseHessian: using std::set<size_t> for sparsity and p.size() "
     628                "not equal domain dimension for f"
     629        );
     630 
     631        sparse_set sparsity;
     632        if( work.color.size() == 0 )
     633        {       bool transpose = false;
     634                vec_set_to_sparse_set(sparsity, p, n, n, transpose);
     635        }
     636       
     637        // compute the Hessian
     638        size_t n_sweep = SparseHessianCompute(x, w, sparsity, hes, work);
     639
     640        return n_sweep;
     641}
     642// ===========================================================================
    174643/*!
    175644Private helper function for SparseHessian(x, w, p).
     
    207676        VectorBase&        hes              )
    208677{
    209         typedef CppAD::vector<size_t> SizeVector;
    210         typedef CppAD::vectorBool     VectorBool;
    211         size_t j, k, l;
    212 
    213678        size_t n = Domain();
     679        size_t i, j, k;
    214680
    215681        // check Vector is Simple Vector class with bool elements
     
    232698                "not equal square of domain dimension for f"
    233699        );
    234        
    235         // initial coloring
    236         SizeVector color(n);
    237         for(j = 0; j < n; j++)
    238                 color[j] = j;
    239 
    240         // See GreedyPartialD2Coloring Algorithm Section 3.6.2 of
    241         // Graph Coloring in Optimization Revisited by
    242         // Assefaw Gebremedhin, Fredrik Maane, Alex Pothen
    243         VectorBool forbidden(n);
    244         for(j = 0; j < n; j++)
    245         {       // initial all colors as ok for this column
    246                 for(k = 0; k < n; k++)
    247                         forbidden[k] = false;
    248                 // for each row that is connected to column j
    249                 for(l = 0; l < n; l++) if( p[l * n + j] )
    250                 {       // for each column that is connected to row l
    251                         for(k = 0; k < n; k++) if( p[l * n + k] & (j != k) )   
    252                                 forbidden[ color[k] ] = true;
    253                 }
    254                 k = 0;
    255                 while( forbidden[k] && k < n )
    256                 {       k++;
    257                         CPPAD_ASSERT_UNKNOWN( k < n );
    258                 }
    259                 color[j] = k;
    260         }
    261         size_t n_color = 1;
    262         for(k = 0; k < n; k++) n_color = std::max(n_color, color[k] + 1);
    263 
    264         // some values
    265         const Base zero(0);
    266         const Base one(1);
    267 
    268         // point at which we are evaluating the Hessian
    269         Forward(0, x);
    270 
    271         // initialize the return value
    272         for(j = 0; j < n; j++)
    273                 for(k = 0; k < n; k++)
    274                         hes[j * n + k] = zero;
    275 
    276         // direction vector for calls to forward
    277         VectorBase u(n);
    278 
    279         // location for return values from Reverse
    280         VectorBase ddw(n * 2);
    281 
    282         // loop over colors
    283         size_t c;
    284         for(c = 0; c < n_color; c++)
    285         {       // determine all the colums with this color
    286                 for(j = 0; j < n; j++)
    287                 {       if( color[j] == c )
    288                                 u[j] = one;
    289                         else    u[j] = zero;
    290                 }
    291                 // call forward mode for all these columns at once
    292                 Forward(1, u);
    293 
    294                 // evaluate derivative of w^T * F'(x) * u
    295                 ddw = Reverse(2, w);
    296 
    297                 // set the corresponding components of the result
    298                 for(j = 0; j < n; j++) if( color[j] == c )
    299                 {       for(l = 0; l < n; l++)
    300                                 if( p[ l * n + j ] )
    301                                         hes[l * n + j] = ddw[l * 2 + 1];
    302                 }
    303         }
     700        sparse_hessian_work work;
     701        CppAD::vector<size_t>& row(work.r_sort);
     702        CppAD::vector<size_t>& col(work.c_sort);
     703        CppAD::vector<size_t>& user_k(work.k_sort);
     704
     705        k = 0;
     706        for(i = 0; i < n; i++)
     707        {       for(j = 0; j < n; j++)
     708                        if( p[i * n + j] )
     709                                k++;
     710        }
     711        size_t K = k;
     712        VectorBase H(K);
     713        row.resize(K+1);
     714        col.resize(K);
     715        user_k.resize(K);
     716        k = 0;
     717        for(i = 0; i < n; i++)
     718        {       for(j = 0; j < n; j++)
     719                        if( p[i * n + j] )
     720                        {       row[k]      = i;
     721                                col[k]      = j;
     722                                user_k[k] = k; 
     723                                k++;
     724                        }
     725        }
     726        row[K] = n;
     727 
     728        SparseHessianCase(set_type, x, w, p, H, work);
     729
     730        Base zero(0);
     731        for(i = 0; i < n; i++)
     732        {       for(j = 0; j < n; j++)
     733                        hes[i * n + j] = zero;
     734        }
     735        for(k = 0; k < K; k++)
     736                hes[row[k] * n + col[k]] = H[k];
    304737}
    305738/*!
     
    339772        VectorBase&              hes              )
    340773{
    341         typedef CppAD::vector<size_t> SizeVector;
    342         typedef CppAD::vectorBool     VectorBool;
    343         size_t j, k, l;
    344 
    345774        size_t n = Domain();
     775        size_t i, j, k;
    346776
    347777        // check VectorSet is Simple Vector class with sets for elements
     
    363793        CPPAD_ASSERT_KNOWN(
    364794                p.size() == n,
    365                 "SparseHessian: using set values and size of p "
     795                "SparseHessian: using size_t sets for sparsity pattern and p.size() "
    366796                "not equal domain dimension for f"
    367797        );
    368        
    369         // initial coloring
    370         SizeVector color(n);
    371         for(j = 0; j < n; j++)
    372                 color[j] = j;
    373 
    374         // See GreedyPartialD2Coloring Algorithm Section 3.6.2 of
    375         // Graph Coloring in Optimization Revisited by
    376         // Assefaw Gebremedhin, Fredrik Maane, Alex Pothen
    377         VectorBool forbidden(n);
    378         std::set<size_t>::const_iterator itr_k, itr_l;
    379         for(j = 0; j < n; j++)
    380         {       // initial all colors as ok for this column
    381                 for(k = 0; k < n; k++)
    382                         forbidden[k] = false;
    383 
    384                 // for each row that is connected to column j
    385                 itr_l = p[j].begin();
    386                 while( itr_l != p[j].end() )
    387                 {       l = *itr_l++;
    388                         // for each column that is connected to row l
    389                         // (same as each row that is connect to column l)
    390                         itr_k = p[l].begin();
    391                         while( itr_k != p[l].end() )
    392                         {       k = *itr_k++;
    393                                 forbidden[ color[k] ] = (j != k);
    394                         }
     798        sparse_hessian_work work;
     799        CppAD::vector<size_t>& row(work.r_sort);
     800        CppAD::vector<size_t>& col(work.c_sort);
     801        CppAD::vector<size_t>& user_k(work.k_sort);
     802
     803        k = 0;
     804        std::set<size_t>::const_iterator itr;
     805        for(i = 0; i < n; i++)
     806        {       itr = p[i].begin();
     807                while( itr != p[i].end() )
     808                {       itr++;
     809                        k++;
    395810                }
    396                 k = 0;
    397                 while( forbidden[k] && k < n )
    398                 {       k++;
    399                         CPPAD_ASSERT_UNKNOWN( k < n );
     811        }
     812        size_t K = k;
     813        VectorBase H(K);
     814        row.resize(K+1);
     815        col.resize(K);
     816        user_k.resize(K);
     817        k = 0;
     818        for(i = 0; i < n; i++)
     819        {       itr = p[i].begin();
     820                while( itr != p[i].end() )
     821                {       row[k]      = i;
     822                        col[k]      = *itr;
     823                        user_k[k] = k; 
     824                        itr++;
     825                        k++;
    400826                }
    401                 color[j] = k;
    402         }
    403         size_t n_color = 1;
    404         for(k = 0; k < n; k++) n_color = std::max(n_color, color[k] + 1);
    405 
    406         // some values
    407         const Base zero(0);
    408         const Base one(1);
    409 
    410         // point at which we are evaluating the Hessian
    411         Forward(0, x);
    412 
    413         // initialize the return value
    414         for(j = 0; j < n; j++)
    415                 for(k = 0; k < n; k++)
    416                         hes[j * n + k] = zero;
    417 
    418         // direction vector for calls to forward
    419         VectorBase u(n);
    420 
    421         // location for return values from Reverse
    422         VectorBase ddw(n * 2);
    423 
    424         // loop over colors
    425         size_t c;
    426         for(c = 0; c < n_color; c++)
    427         {       // determine all the colums with this color
    428                 for(j = 0; j < n; j++)
    429                 {       if( color[j] == c )
    430                                 u[j] = one;
    431                         else    u[j] = zero;
     827        }
     828        row[K] = n;
     829 
     830        SparseHessianCase(set_type, x, w, p, H, work);
     831
     832        Base zero(0);
     833        for(i = 0; i < n; i++)
     834        {       for(j = 0; j < n; j++)
     835                        hes[i * n + j] = zero;
     836        }
     837        for(k = 0; k < K; k++)
     838                hes[row[k] * n + col[k]] = H[k];
     839}
     840// ===========================================================================
     841// Public Member Functions
     842// ===========================================================================
     843/*!
     844Compute user specified subset of a sparse Hessian.
     845
     846The C++ source code corresponding to this operation is
     847\verbatim
     848        SparceHessian(x, w, p, row, col, hes, work)
     849\endverbatim
     850
     851\tparam Base
     852is the base type for the recording that is stored in this ADFun<Base object.
     853
     854\tparam VectorBase
     855is a simple vector class with elements of type \a Base.
     856
     857\tparam VectorSet
     858is a simple vector class with elements of type
     859\c bool or \c std::set<size_t>.
     860
     861\tparam VectorSize
     862is a simple vector class with elements of type \c size_t.
     863
     864\param x
     865is a vector specifing the point at which to compute the Hessian.
     866
     867\param w
     868is the weighting vector that defines a scalar valued function by
     869a weighted sum of the components of the vector valued function
     870$latex F(x)$$.
     871
     872\param p
     873is the sparsity pattern for the Hessian that we are calculating.
     874
     875\param row
     876is the vector of row indices for the returned Hessian values.
     877
     878\param col
     879is the vector of columns indices for the returned Hessian values.
     880It must have the same size are r.
     881
     882\param hes
     883is the vector of Hessian values.
     884It must have the same size are r.
     885The return value <code>hes[k]</code> is the second partial of
     886\f$ w^{\rm T} F(x)\f$ with respect to the
     887<code>row[k]</code> and <code>col[k]</code> component of \f$ x\f$.
     888
     889\param work
     890contains information that depends on the function object, sparsity pattern,
     891\c row, and \c col vector.
     892If these values are the same, \c work does not need to be recomputed.
     893To be more specific,
     894\c r_sort is sorted copy of \c row ,
     895\c c_sort is sorted copy of \c col ,
     896<code>k_sort[k]</code> is the original index corresponding to the
     897values <code>r_sort[k]</code> and <code>c_sort[k]</code>.
     898The order for the sort is by columns.
     899Let \c n the domain dimension,
     900and \c K the size of \c row , \c col , and \c hes.
     901There is one extra entry
     902in the sorted row array and it has value <code>r_sort[K]=n</code>.
     903The \c color vector is set and used by \c SparseHessianCompute.
     904
     905\return
     906Is the number of first order forward sweeps used to compute the
     907requested Hessian values.
     908(This is also equal to the number of second order reverse sweeps.)
     909The total work, not counting the zero order
     910forward sweep, or the time to combine computations, is proportional to this
     911return value.
     912*/
     913template<class Base>
     914template <class VectorBase, class VectorSet, class VectorSize>
     915size_t ADFun<Base>::SparseHessian(
     916        const VectorBase&     x    ,
     917        const VectorBase&     w    ,
     918        const VectorSet&      p    ,
     919        const VectorSize&     row  ,
     920        const VectorSize&     col  ,
     921        VectorBase&           hes  ,
     922        sparse_hessian_work&  work )
     923{
     924        size_t n = Domain();
     925        size_t k, K = hes.size();
     926        if( work.r_sort.size() == 0 )
     927        {       // create version of (row, col,k) sorted by rows
     928                size_t min_bytes = 3 * K * sizeof(size_t);
     929                size_t cap_bytes;
     930                void*  v_ptr  = thread_alloc::get_memory(min_bytes, cap_bytes);
     931                size_t* rck  = reinterpret_cast<size_t*>(v_ptr);
     932                for(k = 0; k < K; k++)
     933                {       // must put column first so it is used for sorting
     934                        rck[3 * k + 0] = row[k];
     935                        rck[3 * k + 1] = col[k];
     936                        rck[3 * k + 2] = k;
    432937                }
    433                 // call forward mode for all these columns at once
    434                 Forward(1, u);
    435 
    436                 // evaluate derivative of w^T * F'(x) * u
    437                 ddw = Reverse(2, w);
    438 
    439                 // set the corresponding components of the result
    440                 for(j = 0; j < n; j++) if( color[j] == c )
    441                 {       itr_l = p[j].begin();
    442                         while( itr_l != p[j].end() )
    443                         {       l = *itr_l++;
    444                                 hes[l * n + j] = ddw[l * 2 + 1];
    445                         }
     938                std::qsort(rck, K, 3 * sizeof(size_t), std_qsort_compare<size_t>);
     939                work.c_sort.resize(K);
     940                work.r_sort.resize(K+1);
     941                work.k_sort.resize(K);
     942                for(k = 0; k < K; k++)
     943                {       work.r_sort[k] = rck[3 * k + 0]; 
     944                        work.c_sort[k] = rck[3 * k + 1]; 
     945                        work.k_sort[k] = rck[3 * k + 2]; 
    446946                }
    447         }
     947                work.r_sort[K] = n;
     948                thread_alloc::return_memory(v_ptr);
     949        }
     950# ifndef NDEBUG
     951        CPPAD_ASSERT_KNOWN(
     952                x.size() == n ,
     953                "SparseHessian: size of x not equal domain dimension for f."
     954        );
     955        CPPAD_ASSERT_KNOWN(
     956                row.size() == K && col.size() == K ,
     957                "SparseHessian: either r or c does not have the same size as ehs."
     958        );
     959        CPPAD_ASSERT_KNOWN(
     960                work.r_sort.size() == K+1 &&
     961                work.c_sort.size() == K   &&
     962                work.k_sort.size() == K   ,
     963                "SparseHessian: invalid value in work."
     964        );
     965        CPPAD_ASSERT_KNOWN(
     966                work.color.size() == 0 || work.color.size() == n,
     967                "SparseHessian: invalid value in work."
     968        );
     969        for(k = 0; k < K; k++)
     970        {       CPPAD_ASSERT_KNOWN(
     971                        row[k] < n,
     972                        "SparseHessian: invalid value in r."
     973                );
     974                CPPAD_ASSERT_KNOWN(
     975                        col[k] < n,
     976                        "SparseHessian: invalid value in c."
     977                );
     978                CPPAD_ASSERT_KNOWN(
     979                        work.k_sort[k] < K,
     980                        "SparseHessian: invalid value in work."
     981                );
     982                CPPAD_ASSERT_KNOWN(
     983                        work.r_sort[k] == row[ work.k_sort[k] ] ,
     984                        "SparseHessian: invalid value in work."
     985                );
     986                CPPAD_ASSERT_KNOWN(
     987                        work.c_sort[k] == col[ work.k_sort[k] ],
     988                        "SparseHessian: invalid value in work."
     989                );
     990        }
     991        if( work.color.size() != 0 )
     992                for(size_t j = 0; j < n; j++) CPPAD_ASSERT_KNOWN(
     993                        work.color[j] < n,
     994                        "SparseHessian: invalid value in work."
     995        );
     996# endif
     997 
     998        typedef typename VectorSet::value_type Set_type;
     999        size_t n_sweep = SparseHessianCase(Set_type(), x, w, p, hes, work);
     1000        return n_sweep;
    4481001}
    449 
    4501002/*!
    4511003Compute a sparse Hessian.
     
    4531005The C++ source code coresponding to this operation is
    4541006\verbatim
    455         hes = SparseHessian(x, w)
     1007        hes = SparseHessian(x, w, p)
    4561008\endverbatim
    4571009
     
    4981050        return hes;
    4991051}
    500 
    5011052/*!
    5021053Compute a sparse Hessian
  • trunk/cppad/local/sparse_jacobian.hpp

    r2006 r2401  
    44
    55/* --------------------------------------------------------------------------
    6 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-11 Bradley M. Bell
     6CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
    77
    88CppAD is distributed under multiple licenses. This distribution is under
     
    1313Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
    1414-------------------------------------------------------------------------- */
    15 
    1615/*
    1716$begin sparse_jacobian$$
    1817$spell
     18        recomputed
    1919        valarray
    2020        std
     
    3232
    3333$head Syntax$$
    34 $codei%%jac% = %f%.SparseJacobian(%x%)
     34$icode%jac% = %f%.SparseJacobian(%x%)
     35%jac% = %f%.SparseJacobian(%x%, %p%)
     36%n_sweep% = %f%.SparseJacobianForward(%x%, %p%, %row%, %col%, %jac%, %work%)
     37%n_sweep% = %f%.SparseJacobianReverse(%x%, %p%, %row%, %col%, %jac%, %work%)
    3538%$$
    36 $codei%%jac% = %f%.SparseJacobian(%x%, %p%)%$$
    3739
    3840$head Purpose$$
     41We use $latex n$$ for the $cref/domain/seq_property/Domain/$$ size,
     42and $latex m$$ for the $cref/range/seq_property/Range/$$ size of $icode f$$.
    3943We use $latex F : \R^n \rightarrow \R^m$$ do denote the
    4044$cref/AD function/glossary/AD Function/$$
     
    4448        jac = F^{(1)} (x)
    4549\] $$
    46 This routine assumes
    47 that the matrix $latex F^{(1)} (x) \in \R^{m \times n}$$ is sparse
    48 and uses this assumption to reduce the amount of computation necessary.
    49 One should use speed tests (e.g. $cref/speed_test/$$)
     50This routine takes advantage of the sparsity of the Jacobian
     51in order to reduce the amount of computation necessary.
     52If $icode row$$ and $icode col$$ are present, it also takes
     53advantage of the reduced set of elements of the Jacobian that
     54need to be computed.
     55One can use speed tests (e.g. $cref speed_test$$)
    5056to verify that results are computed faster
    51 than when using the routine $cref/Jacobian/$$.
     57than when using the routine $cref Jacobian$$.
    5258
    5359$head f$$
     
    5662        ADFun<%Base%> %f%
    5763%$$
    58 Note that the $cref/ADFun/$$ object $icode f$$ is not $code const$$
     64Note that the $cref ADFun$$ object $icode f$$ is not $code const$$
    5965(see $cref/Uses Forward/sparse_jacobian/Uses Forward/$$ below).
    6066
     
    7379$head p$$
    7480The argument $icode p$$ is optional and has prototype
    75 $syntax%
     81$codei%
    7682        const %VectorSet%& %p%
    7783%$$
     
    9096If this sparsity pattern does not change between calls to
    9197$codei SparseJacobian$$, it should be faster to calculate $icode p$$ once
    92 (using $cref/ForSparseJac/$$ or $cref/RevSparseJac/$$)
     98(using $cref ForSparseJac$$ or $cref/RevSparseJac/$$)
    9399and then pass $icode p$$ to $codei SparseJacobian$$.
    94100In addition,
     
    100106for the internal calculations is unspecified.
    101107
     108$head row, col$$
     109The arguments $icode row$$ and $icode col$$ are optional and have prototype
     110$codei%
     111        const %VectorSize%& %row%
     112        const %VectorSize%& %col%
     113%$$
     114(see $cref/VectorSize/sparse_jacobian/VectorSize/$$ below).
     115They specify which rows and columns of $latex F^{(1)} (x)$$ are
     116returned and in what order.
     117We use $latex K$$ to denote the value $icode%jac%.size()%$$
     118which must also equal the size of $icode row$$ and $icode col$$.
     119Furthermore,
     120for $latex k = 0 , \ldots , K-1$$, it must hold that
     121$latex row[k] < m$$ and $latex col[k] < n$$.
     122In addition,
     123all of the $latex (row[k], col[k])$$ pairs must correspond to a true value
     124in the sparsity pattern $icode p$$.
     125
    102126$head jac$$
    103127The result $icode jac$$ has prototype
    104128$codei%
    105         %VectorBase% %jac%
     129        %VectorBase%& %jac%
    106130%$$
    107 and its size is $latex m * n$$.
    108 For $latex i = 0 , \ldots , m - 1$$,
    109 and $latex j = 0 , \ldots , n - 1 $$
     131In the case where the arguments $icode row$$ and $icode col$$ are not present,
     132the size of $icode jac$$ is $latex m * n$$ and
     133for $latex i = 0 , \ldots , m-1$$,
     134$latex j = 0 , \ldots , n-1$$,
    110135$latex \[
    111136        jac [ i * n + j ] = \D{ F_i }{ x_j } (x)
    112137\] $$
     138$pre
     139
     140$$
     141In the case where the arguments $icode row$$ and $icode col$$ are present,
     142we use $latex K$$ to denote the size of $icode jac$$.
     143The input value of its elements does not matter.
     144Upon return, for $latex k = 0 , \ldots , K - 1$$,
     145$latex \[
     146        jac [ k ] = \D{ F_i }{ x_j } (x)
     147        \; , \;
     148        \; {\rm where} \;
     149        i = row[k]
     150        \; {\rm and } \;
     151        j = col[k]
     152\] $$
     153
     154$head work$$
     155If this argument is present, it has prototype
     156$icode%
     157        sparse_jacobian_work& %work%
     158%$$
     159This object can only be used with the routines
     160$code SparseJacobianForward$$ and $code SparseJacobianReverse$$.
     161During its the first use, information is stored in $icode work$$.
     162This is used to reduce the work done by future calls to the same mode
     163(forward or reverse),
     164the same $icode f$$, $icode p$$, $icode row$$, and $icode col$$.
     165If a future call is for a different mode,
     166or any of these values have changed,
     167you must first call $icode%work%.clear()%$$
     168to inform CppAD that this information needs to be recomputed.
     169
     170$head n_sweep$$
     171The return value $icode n_sweep$$ has prototype
     172$codei%
     173        size_t %n_sweep%
     174%$$
     175If $code SparseJacobianForward$$ ($code SparseJacobianReverse$$) is used,
     176$icode n_sweep$$ is the number of first order forward (reverse) sweeps
     177used to compute the requested Jacobian values.
     178This is proportional to the total work that $code SparseJacobian$$ does,
     179not counting the zero order forward sweep,
     180or the work to combine multiple columns (rows) into a single sweep.
    113181
    114182$head VectorBase$$
    115 The type $icode VectorBase$$ must be a $cref/SimpleVector/$$ class with
     183The type $icode VectorBase$$ must be a $cref SimpleVector$$ class with
    116184$cref/elements of type/SimpleVector/Elements of Specified Type/$$
    117185$icode Base$$.
    118 The routine $cref/CheckSimpleVector/$$ will generate an error message
     186The routine $cref CheckSimpleVector$$ will generate an error message
    119187if this is not the case.
    120188
    121189$head VectorSet$$
    122 The type $icode VectorSet$$ must be a $xref/SimpleVector/$$ class with
    123 $xref/SimpleVector/Elements of Specified Type/elements of type/$$
     190The type $icode VectorSet$$ must be a $cref SimpleVector$$ class with
     191$cref/elements of type/SimpleVector/Elements of Specified Type/$$
    124192$code bool$$ or $code std::set<size_t>$$;
    125193see $cref/sparsity pattern/glossary/Sparsity Pattern/$$ for a discussion
    126194of the difference.
    127 The routine $cref/CheckSimpleVector/$$ will generate an error message
     195The routine $cref CheckSimpleVector$$ will generate an error message
    128196if this is not the case.
    129197
     
    136204this condition.
    137205
     206$head VectorSize$$
     207The type $icode VectorSize$$ must be a $cref SimpleVector$$ class with
     208$cref/elements of type/SimpleVector/Elements of Specified Type/$$
     209$code size_t$$.
     210The routine $cref CheckSimpleVector$$ will generate an error message
     211if this is not the case.
     212
    138213$head Uses Forward$$
    139 After each call to $cref/Forward/$$,
     214After each call to $cref Forward$$,
    140215the object $icode f$$ contains the corresponding
    141216$cref/Taylor coefficients/glossary/Taylor Coefficient/$$.
    142217After $code SparseJacobian$$,
    143 the previous calls to $xref/Forward/$$ are undefined.
     218the previous calls to $cref Forward$$ are undefined.
    144219
    145220$head Example$$
     
    148223%$$
    149224The routine
    150 $cref/sparse_jacobian.cpp/$$
     225$cref sparse_jacobian.cpp$$
    151226is examples and tests of $code sparse_jacobian$$.
    152227It return $code true$$, if it succeeds and $code false$$ otherwise.
    153228
    154229$end
    155 -----------------------------------------------------------------------------
     230==============================================================================
    156231*/
    157232# include <cppad/local/std_set.hpp>
     
    162237Sparse Jacobian driver routine and helper functions.
    163238*/
    164 
     239// ===========================================================================
    165240/*!
    166 Private helper function for SparseJacobian(x, p).
    167 
    168 All of the description in the public member function SparseJacobian(x, p)
    169 applies.
     241class used by SparseJacobian to hold information
     242so it does not need to be recomputed.
     243*/
     244class sparse_jacobian_work {
     245        public:
     246                /// version of user r array sorted by row or column
     247                CppAD::vector<size_t> r_sort;
     248                /// version of user c array sorted by row or column
     249                CppAD::vector<size_t> c_sort;
     250                /// mapping from sorted array indices to user array indices
     251                CppAD::vector<size_t> k_sort;
     252                /// results of the coloring algorithm
     253                CppAD::vector<size_t> color;
     254                /// inform CppAD that this information needs to be recomputed
     255                void clear(void)
     256                {       r_sort.resize(0);
     257                        c_sort.resize(0);
     258                        k_sort.resize(0);
     259                        color.resize(0);
     260                }
     261};
     262// ===========================================================================
     263/*!
     264std::qsort comparision function
     265
     266\tparam Key
     267Is the type used for the key for the comparision operation.
     268There must be a conversion operation from the type
     269\c Key to the type \c int.
     270
     271\param a_vptr
     272is a pointer to first \c Key value to be compared
     273(this must be the beginning of the corresponding record being sorted).
     274
     275\param b_vptr
     276is a pointer to the second \c Key value to be compared
     277(this must be the beginning of the corresponding record being sorted).
     278
     279\return
     280is negative, zero, position if the first value is less than,
     281equal, or greater than the second (respectively).
     282*/
     283template <class Key>
     284int std_qsort_compare(const void* a_vptr, const void* b_vptr)
     285{       const Key* a_ptr = reinterpret_cast<const Key*>(a_vptr);
     286        const Key* b_ptr = reinterpret_cast<const Key*>(b_vptr);
     287        int a = static_cast<int>(*a_ptr);
     288        int b = static_cast<int>(*b_ptr);
     289        return a - b;
     290}
     291// ===========================================================================
     292/*!
     293Private helper function forward mode cases
     294
     295\tparam Base
     296See \c SparseJacobianForward(x, p, row, col, jac, work).
     297
     298\tparam VectorBase
     299See \c SparseJacobianForward(x, p, row, col, jac, work).
     300
     301\tparam VectorSet
     302is either \c sparse_pack or \c sparse_set.
     303
     304\param x
     305See \c SparseJacobianForward(x, p, row, col, jac, work).
     306
     307\param p_transpose
     308If <code>work.color.size() != 0</code>, then \c p_transpose is not used.
     309Otherwise, it is a
     310sparsity pattern for the transpose of the Jacobian of this ADFun<Base> object.
     311Note that we do not change the values in \c p_transpose,
     312but is not \c const because we use its iterator facility.
     313
     314
     315\param jac
     316See \c SparseJacobianForward(x, p, row, col, jac, work).
     317
     318\param work
     319See \c SparseJacobianForward(x, p, row, col, jac, work).
     320
     321\return
     322See \c SparseJacobianForward(x, p, row, col, jac, work).
     323*/
     324template<class Base>
     325template <class VectorBase, class VectorSet>
     326size_t ADFun<Base>::SparseJacobianFor(
     327        const VectorBase&     x           ,
     328        VectorSet&            p_transpose ,
     329        VectorBase&           jac         ,
     330        sparse_jacobian_work& work        )
     331{
     332        using   CppAD::vectorBool;
     333        size_t i, j, k, ell;
     334
     335        CppAD::vector<size_t>& row(work.r_sort);
     336        CppAD::vector<size_t>& col(work.c_sort);
     337        CppAD::vector<size_t>& user_k(work.k_sort);
     338        CppAD::vector<size_t>& color(work.color);
     339
     340        size_t m = Range();
     341        size_t n = Domain();
     342
     343        // some values
     344        const Base zero(0);
     345        const Base one(1);
     346
     347        // check VectorBase is Simple Vector class with Base type elements
     348        CheckSimpleVector<Base, VectorBase>();
     349
     350        CPPAD_ASSERT_UNKNOWN( x.size() == n );
     351        CPPAD_ASSERT_UNKNOWN( color.size() == 0 || color.size() == n );
     352
     353        // number of components of Jacobian that are required
     354        size_t K = jac.size();
     355        CPPAD_ASSERT_UNKNOWN( row.size() == K );
     356        CPPAD_ASSERT_UNKNOWN( col.size() == K+1 );
     357        CPPAD_ASSERT_UNKNOWN( col[K] == n );
     358
     359        // Point at which we are evaluating the Jacobian
     360        Forward(0, x);
     361
     362        if( color.size() == 0 )
     363        {
     364                CPPAD_ASSERT_UNKNOWN( p_transpose.n_set() ==  n );
     365                CPPAD_ASSERT_UNKNOWN( p_transpose.end() ==  m );
     366
     367                // rows and columns that are in the returned jacobian
     368                VectorSet r_used, c_used;
     369                r_used.resize(n, m);
     370                c_used.resize(m, n);
     371                k = 0;
     372                while( k < K )
     373                {       CPPAD_ASSERT_UNKNOWN( row[k] < m && col[k] < n );
     374                        CPPAD_ASSERT_UNKNOWN( k == 0 || col[k-1] <= col[k] );
     375                        CPPAD_ASSERT_KNOWN(
     376                                p_transpose.is_element(col[k], row[k]) ,
     377                                "SparseJacobianForward: "
     378                                "an (row, col) pair is not in sparsity pattern."
     379                        );
     380                        r_used.add_element(col[k], row[k]);
     381                        c_used.add_element(row[k], col[k]);
     382                        k++;
     383                }
     384       
     385                // given a row index, which columns are non-zero and not used
     386                VectorSet not_used;
     387                not_used.resize(m, n);
     388                for(j = 0; j < n; j++)
     389                {       p_transpose.begin(j);
     390                        i = p_transpose.next_element();
     391                        while( i != p_transpose.end() )
     392                        {       if( ! c_used.is_element(i, j) )
     393                                        not_used.add_element(i, j);
     394                                i = p_transpose.next_element();
     395                        }
     396                }
     397
     398                // initial coloring
     399                color.resize(n);
     400                for(j = 0; j < n; j++)
     401                        color[j] = j;
     402       
     403                // See GreedyPartialD2Coloring Algorithm Section 3.6.2 of
     404                // Graph Coloring in Optimization Revisited by
     405                // Assefaw Gebremedhin, Fredrik Maane, Alex Pothen
     406                vectorBool forbidden(n);
     407                for(j = 1; j < n; j++)
     408                {
     409                        // initial all colors as ok for this column
     410                        // (value of forbidden for ell > j does not matter)
     411                        for(ell = 0; ell <= j; ell++)
     412                                forbidden[ell] = false;
     413       
     414                        // for each row that is non-zero for this column
     415                        p_transpose.begin(j);
     416                        i = p_transpose.next_element();
     417                        while( i != p_transpose.end() )
     418                        {       // for each column that this row uses
     419                                c_used.begin(i);
     420                                ell = c_used.next_element();
     421                                while( ell != c_used.end() )
     422                                {       // if this is not the same column, forbid its color
     423                                        if( ell < j )
     424                                                forbidden[ color[ell] ] = true;
     425                                        ell = c_used.next_element();
     426                                }
     427                                i = p_transpose.next_element();
     428                        }
     429       
     430                        // for each row that this column uses
     431                        r_used.begin(j);
     432                        i = r_used.next_element();
     433                        while( i != r_used.end() )
     434                        {       // For each column that is non-zero for this row
     435                                // (the used columns have already been checked above).
     436                                not_used.begin(i);
     437                                ell = not_used.next_element();
     438                                while( ell != not_used.end() )
     439                                {       // if this is not the same column, forbid its color
     440                                        if( ell < j )
     441                                                forbidden[ color[ell] ] = true;
     442                                        ell = not_used.next_element();
     443                                }
     444                                i = r_used.next_element();
     445                        }
     446
     447                        // pick the color with smallest index
     448                        ell = 0;
     449                        while( forbidden[ell] )
     450                        {       ell++;
     451                                CPPAD_ASSERT_UNKNOWN( ell <= j );
     452                        }
     453                        color[j] = ell;
     454                }
     455        }
     456        size_t n_color = 1;
     457        for(ell = 0; ell < n; ell++)
     458                n_color = std::max(n_color, color[ell] + 1);
     459
     460        // direction vector for calls to forward
     461        VectorBase dx(n);
     462
     463        // location for return values from forward
     464        VectorBase dy(m);
     465
     466        // initialize the return value
     467        for(k = 0; k < K; k++)
     468                jac[k] = zero;
     469
     470        // loop over colors
     471        size_t n_sweep = 0;
     472        for(ell = 0; ell < n_color; ell++)
     473        {       bool any = false;
     474                k = 0;
     475                for(j = 0; j < n; j++) if( color[j] == ell )
     476                {       // find first k such that col[k] has color ell
     477                        if( ! any )
     478                        {       while( col[k] < j )
     479                                        k++;
     480                                any = col[k] == j;
     481                        }
     482                }
     483                if( any )
     484                {       n_sweep++;
     485                        // combine all columns with this color
     486                        for(j = 0; j < n; j++)
     487                        {       dx[j] = zero;
     488                                if( color[j] == ell )
     489                                        dx[j] = one;
     490                        }
     491                        // call forward mode for all these columns at once
     492                        dy = Forward(1, dx);
     493
     494                        // set the corresponding components of the result
     495                        for(j = 0; j < n; j++) if( color[j] == ell )
     496                        {       // find first index in c for this column
     497                                while( col[k] < j )
     498                                        k++;
     499                                // extract the row results for this column
     500                                while( col[k] == j )
     501                                {       jac[ user_k[k] ] = dy[row[k]];
     502                                        k++;
     503                                }
     504                        }
     505                }
     506        }
     507        return n_sweep;
     508}
     509/*!
     510Private helper function for reverse mode cases.
     511
     512\tparam Base
     513See \c SparseJacobianForward(x, p, row, col, jac, work).
     514
     515\tparam VectorBase
     516See \c SparseJacobianForward(x, p, row, col, jac, work).
     517
     518\tparam VectorSet
     519is either \c sparse_pack or \c sparse_set.
     520
     521\param x
     522See \c SparseJacobianForward(x, p, row, col, jac, work).
     523
     524\param p
     525If <code>work.color.size() != 0</code>, then \c p is not used.
     526Otherwise, it is a
     527sparsity pattern for the Jacobian of this ADFun<Base> object.
     528Note that we do not change the values in \c p_transpose,
     529but is not \c const because we use its iterator facility.
     530
     531\param jac
     532See \c SparseJacobianForward(x, p, row, col, jac, work).
     533
     534\param work
     535See \c SparseJacobianForward(x, p, row, col, jac, work).
     536
     537\return
     538See \c SparseJacobianForward(x, p, row, col, jac, work).
     539*/
     540template<class Base>
     541template <class VectorBase, class VectorSet>
     542size_t ADFun<Base>::SparseJacobianRev(
     543        const VectorBase&     x           ,
     544        VectorSet&            p           ,
     545        VectorBase&           jac         ,
     546        sparse_jacobian_work& work        )
     547{
     548        using   CppAD::vectorBool;
     549        size_t i, j, k, ell;
     550
     551        CppAD::vector<size_t>& row(work.r_sort);
     552        CppAD::vector<size_t>& col(work.c_sort);
     553        CppAD::vector<size_t>& user_k(work.k_sort);
     554        CppAD::vector<size_t>& color(work.color);
     555
     556        size_t m = Range();
     557        size_t n = Domain();
     558
     559        // some values
     560        const Base zero(0);
     561        const Base one(1);
     562
     563        // check VectorBase is Simple Vector class with Base type elements
     564        CheckSimpleVector<Base, VectorBase>();
     565
     566        CPPAD_ASSERT_UNKNOWN( x.size() == n );
     567        CPPAD_ASSERT_UNKNOWN (color.size() == m || color.size() == 0 );
     568
     569        // number of components of Jacobian that are required
     570        size_t K = jac.size();
     571        CPPAD_ASSERT_UNKNOWN( row.size() == K+1 );
     572        CPPAD_ASSERT_UNKNOWN( col.size() == K );
     573        CPPAD_ASSERT_UNKNOWN( row[K] == m );
     574
     575        // Point at which we are evaluating the Jacobian
     576        Forward(0, x);
     577
     578        if( color.size() == 0 )
     579        {       CPPAD_ASSERT_UNKNOWN( p.n_set() ==  m );
     580                CPPAD_ASSERT_UNKNOWN( p.end() ==  n );
     581
     582                // rows and columns that are in the returned jacobian
     583                VectorSet r_used, c_used;
     584                r_used.resize(n, m);
     585                c_used.resize(m, n);
     586                k = 0;
     587                while( k < K )
     588                {       CPPAD_ASSERT_UNKNOWN( row[k] < m && col[k] < n );
     589                        CPPAD_ASSERT_UNKNOWN( k == 0 || row[k-1] <= row[k] );
     590                        CPPAD_ASSERT_KNOWN(
     591                                p.is_element(row[k], col[k]) ,
     592                                "SparseJacobianReverse: "
     593                                "an (row, col) pair is not in sparsity pattern."
     594                        );
     595                        r_used.add_element(col[k], row[k]);
     596                        c_used.add_element(row[k], col[k]);
     597                        k++;
     598                }
     599       
     600                // given a column index, which rows are non-zero and not used
     601                VectorSet not_used;
     602                not_used.resize(n, m);
     603                for(i = 0; i < m; i++)
     604                {       p.begin(i);
     605                        j = p.next_element();
     606                        while( j != p.end() )
     607                        {       if( ! r_used.is_element(j , i) )
     608                                        not_used.add_element(j, i);
     609                                j = p.next_element();
     610                        }
     611                }
     612       
     613                // initial coloring
     614                color.resize(m);
     615                for(i = 0; i < m; i++)
     616                        color[i] = i;
     617       
     618                // See GreedyPartialD2Coloring Algorithm Section 3.6.2 of
     619                // Graph Coloring in Optimization Revisited by
     620                // Assefaw Gebremedhin, Fredrik Maane, Alex Pothen
     621                vectorBool forbidden(m);
     622                for(i = 1; i < m; i++)
     623                {
     624                        // initial all colors as ok for this row
     625                        // (value of forbidden for ell > i does not matter)
     626                        for(ell = 0; ell <= i; ell++)
     627                                forbidden[ell] = false;
     628       
     629                        // -----------------------------------------------------
     630                        // Forbid colors for which this row would destroy results
     631                        // for each column that is non-zero for this row
     632                        p.begin(i);
     633                        j = p.next_element();
     634                        while( j != p.end() )
     635                        {       // for each row that this column uses
     636                                r_used.begin(j);
     637                                ell = r_used.next_element();
     638                                while( ell != r_used.end() )
     639                                {       // if this is not the same row, forbid its color
     640                                        if( ell < i )
     641                                                forbidden[ color[ell] ] = true;
     642                                        ell = r_used.next_element();
     643                                }
     644                                j = p.next_element();
     645                        }
     646
     647       
     648                        // -----------------------------------------------------
     649                        // Forbid colors that would destroy results for this row.
     650                        // for each column that this row uses
     651                        c_used.begin(i);
     652                        j = c_used.next_element();
     653                        while( j != c_used.end() )
     654                        {       // For each row that is non-zero for this column
     655                                // (the used rows have already been checked above).
     656                                not_used.begin(j);
     657                                ell = not_used.next_element();
     658                                while( ell != not_used.end() )
     659                                {       // if this is not the same row, forbid its color
     660                                        if( ell < i )
     661                                                forbidden[ color[ell] ] = true;
     662                                        ell = not_used.next_element();
     663                                }
     664                                j = c_used.next_element();
     665                        }
     666
     667                        // pick the color with smallest index
     668                        ell = 0;
     669                        while( forbidden[ell] )
     670                        {       ell++;
     671                                CPPAD_ASSERT_UNKNOWN( ell <= i );
     672                        }
     673                        color[i] = ell;
     674                }
     675        }
     676        size_t n_color = 1;
     677        for(ell = 0; ell < m; ell++)
     678                n_color = std::max(n_color, color[ell] + 1);
     679
     680        // weighting vector for calls to reverse
     681        VectorBase w(m);
     682
     683        // location for return values from Reverse
     684        VectorBase dw(n);
     685
     686        // initialize the return value
     687        for(k = 0; k < K; k++)
     688                jac[k] = zero;
     689
     690        // loop over colors
     691        size_t n_sweep = 0;
     692        for(ell = 0; ell < n_color; ell++)
     693        {       bool any = false;
     694                k = 0;
     695                for(i = 0; i < m; i++) if( color[i] == ell )
     696                {       // find first k such that row[k] has color ell
     697                        if( ! any )
     698                        {       while( row[k] < i )
     699                                        k++;
     700                                any = row[k] == i;
     701                        }
     702                }
     703                if( any )
     704                {       n_sweep++;
     705                        // combine all the rows with this color
     706                        for(i = 0; i < m; i++)
     707                        {       w[i] = zero;
     708                                if( color[i] == ell )
     709                                        w[i] = one;
     710                        }
     711                        // call reverse mode for all these rows at once
     712                        dw = Reverse(1, w);
     713
     714                        // set the corresponding components of the result
     715                        for(i = 0; i < m; i++) if( color[i] == ell )
     716                        {       // find first index in r for this row
     717                                while( row[k] < i )
     718                                        k++;
     719                                // extract the row results for this row
     720                                while( row[k] == i )
     721                                {       jac[ user_k[k] ] = dw[col[k]];
     722                                        k++;
     723                                }
     724                        }
     725                }
     726        }
     727        return n_sweep;
     728}
     729// ===========================================================================
     730/*!
     731Private helper function for vector of bool sparsity pattern cases.
     732
     733\tparam Base
     734See \c SparseJacobianForward(x, p, row, col, jac, work).
     735
     736\tparam VectorBase
     737See \c SparseJacobianForward(x, p, row, col, jac, work).
     738
     739\tparam VectorSet
     740is a simple vector with elements of type \c bool.
    170741
    171742\param set_type
    172 is a \c bool value. This argument is used to dispatch to the proper souce
     743has element type for vector representing the sparsity sets.
     744
     745\param x
     746See \c SparseJacobianForward(x, p, row, col, jac, work).
     747
     748\param p
     749Sparsity pattern for the Jacobian of this ADFun<Base> object.
     750
     751\param jac
     752See \c SparseJacobianForward(x, p, row, col, jac, work).
     753
     754\param work
     755See \c SparseJacobianForward(x, p, row, col, jac, work).
     756
     757\return
     758See \c SparseJacobianForward(x, p, row, col, jac, work).
     759*/
     760template <class Base>
     761template <class VectorBase, class VectorSet>
     762size_t ADFun<Base>::SparseJacobianCase(
     763        bool                  set_type    ,
     764        const VectorBase&     x           ,
     765        const VectorSet&      p           ,
     766        VectorBase&           jac         ,
     767        sparse_jacobian_work& work        )
     768{       size_t n = Domain();
     769        size_t m = Range();
     770        size_t n_sweep;
     771
     772        // check VectorSet is Simple Vector class with bool elements
     773        CheckSimpleVector<bool, VectorSet>();
     774
     775        // check VectorBase is Simple Vector class with Base type elements
     776        CheckSimpleVector<Base, VectorBase>();
     777
     778        CPPAD_ASSERT_KNOWN(
     779                p.size() == m * n ,
     780                "SparseJacobian: using bool values and size of p "
     781                " not equal range dimension times domain dimension for f"
     782        );
     783        CPPAD_ASSERT_UNKNOWN( x.size() == n );
     784        size_t K = jac.size();
     785
     786        if( work.c_sort.size() == K+1 )
     787        {       // row major, use forward mode ----------------------------------
     788                CPPAD_ASSERT_UNKNOWN( work.r_sort.size() == K );
     789                CPPAD_ASSERT_UNKNOWN( work.c_sort[K] == n );
     790
     791                // convert the sparsity pattern to a sparse_pack object
     792                // so can fold vector of bools and vector of sets into same function
     793                sparse_pack s_transpose;
     794                if( work.color.size() == 0 )
     795                {       bool transpose = true;
     796                        vec_bool_to_sparse_pack(s_transpose, p, m, n, transpose);
     797                }
     798       
     799                // now we have folded this into the following case
     800                n_sweep = SparseJacobianFor(x, s_transpose, jac, work);
     801        }
     802        else
     803        {       // column major, use reverse mode --------------------------------
     804                CPPAD_ASSERT_UNKNOWN( work.r_sort.size() == K+1 );
     805                CPPAD_ASSERT_UNKNOWN( work.c_sort.size() == K   );
     806                CPPAD_ASSERT_UNKNOWN( work.r_sort[K] == m );
     807
     808                // convert the sparsity pattern to a sparse_pack object
     809                // so can fold vector of bools and vector of sets into same function
     810                sparse_pack s;
     811                if( work.color.size() == 0 )
     812                {       bool transpose = false;
     813                        vec_bool_to_sparse_pack(s, p, m, n, transpose);
     814                }
     815       
     816                // now we have folded this into the following case
     817                n_sweep = SparseJacobianRev(x, s, jac, work);
     818        }
     819        return n_sweep;
     820}
     821/*!
     822Private helper function for vector of std::set<size_t> sparsity pattern cases.
     823
     824\tparam Base
     825See \c SparseJacobianForward(x, p, row, col, jac, work).
     826
     827\tparam VectorBase
     828See \c SparseJacobianForward(x, p, row, col, jac, work).
     829
     830\tparam VectorSet
     831is a simple vector with elements of type <code>std::set<size_t></code>.
     832
     833\param set_type
     834has element type for vector representing the sparsity sets.
     835
     836\param x
     837See \c SparseJacobianForward(x, p, row, col, jac, work).
     838
     839\param p
     840Sparsity pattern for the Jacobian of this ADFun<Base> object.
     841
     842\param jac
     843See \c SparseJacobianForward(x, p, row, col, jac, work).
     844
     845\param work
     846See \c SparseJacobianForward(x, p, row, col, jac, work).
     847
     848\return
     849See \c SparseJacobianForward(x, p, row, col, jac, work).
     850*/
     851template <class Base>
     852template <class VectorBase, class VectorSet>
     853size_t ADFun<Base>::SparseJacobianCase(
     854        const std::set<size_t>&  set_type        ,
     855        const VectorBase&        x               ,
     856        const VectorSet&         p               ,
     857        VectorBase&              jac             ,
     858        sparse_jacobian_work&    work            )
     859{       size_t n = Domain();
     860        size_t m = Range();
     861        size_t n_sweep;
     862
     863        // check VectorSet is Simple Vector class with sets for elements
     864        CheckSimpleVector<std::set<size_t>, VectorSet>(
     865                one_element_std_set<size_t>(), two_element_std_set<size_t>()
     866        );
     867
     868        // check VectorBase is Simple Vector class with Base type elements
     869        CheckSimpleVector<Base, VectorBase>();
     870
     871        CPPAD_ASSERT_KNOWN(
     872                p.size() == Range() ,
     873                "SparseJacobian: using size_t sets for sparsity pattern and p.size() "
     874                "not equal range dimension for f"
     875        );
     876        CPPAD_ASSERT_UNKNOWN( x.size() == n );
     877        size_t K = jac.size();
     878
     879        if( work.c_sort.size() == K+1 )
     880        {       // row major, use forward mode ---------------------------------
     881                CPPAD_ASSERT_UNKNOWN( work.r_sort.size() == K );
     882                CPPAD_ASSERT_UNKNOWN( work.c_sort[K] == n );
     883
     884                // convert the sparsity pattern to a sparse_pack object
     885                // so can fold vector of bools and vector of sets into same function
     886                sparse_set s_transpose;
     887                if( work.color.size() == 0 )
     888                {       bool transpose = true;
     889                        vec_set_to_sparse_set(s_transpose, p, m, n, transpose);
     890                }
     891       
     892                // now we have folded this into the following case
     893                n_sweep = SparseJacobianFor(x, s_transpose, jac, work);
     894        }
     895        else
     896        {       // column major, use reverse mode -----------------------------------
     897                CPPAD_ASSERT_UNKNOWN( work.c_sort.size() == K );
     898                CPPAD_ASSERT_UNKNOWN( work.r_sort.size() == K+1 );
     899                CPPAD_ASSERT_UNKNOWN( work.r_sort[K] == m );
     900
     901                // convert the sparsity pattern to a sparse_pack object
     902                // so can fold vector of bools and vector of sets into same function
     903                sparse_set s;
     904                if( work.color.size() == 0 )
     905                {       bool transpose = false;
     906                        vec_set_to_sparse_set(s, p, m, n, transpose);
     907                }
     908       
     909                // now we have folded this into the following case
     910                n_sweep = SparseJacobianRev(x, s, jac, work);
     911        }
     912        return n_sweep;
     913}
     914// ===========================================================================
     915/*!
     916Private helper function for vector of bool sparsity pattern cases
     917
     918\param set_type
     919is a \c bool value. This argument is used to dispatch to the proper source
    173920code depending on the value of \c VectorSet::value_type.
    174921
     
    200947        size_t n = Domain();
    201948
    202         // some values
     949        // the value zero
    203950        const Base zero(0);
    204         const Base one(1);
    205951
    206952        // check VectorSet is Simple Vector class with bool elements
     
    210956        CheckSimpleVector<Base, VectorBase>();
    211957
    212         CPPAD_ASSERT_KNOWN(
    213                 x.size() == n,
    214                 "SparseJacobian: size of x not equal domain dimension for f"
    215         );
    216958        CPPAD_ASSERT_KNOWN(
    217959                p.size() == m * n,
     
    219961                " not equal range dimension times domain dimension for f"
    220962        );
     963        CPPAD_ASSERT_UNKNOWN( x.size() == n );
    221964        CPPAD_ASSERT_UNKNOWN(jac.size() == m * n);
    222965
    223         // point at which we are evaluating the Jacobian
    224         Forward(0, x);
    225 
     966        // row index, column index, value representation
     967        // (used to fold problem into user interface case).
     968        size_t K = 0;
     969        for(j = 0; j < n; j++)
     970        {       for(i = 0; i < m; i++)
     971                        if( p[ i * n + j ] )
     972                                K++;
     973        }
     974        VectorBase J(K);
     975        sparse_jacobian_work work;
     976        CppAD::vector<size_t>& row(work.r_sort);
     977        CppAD::vector<size_t>& col(work.c_sort);
     978        CppAD::vector<size_t>& user_k(work.k_sort);
     979
     980        if( n <= m )
     981        {       // forward mode, columns are sorted
     982                row.resize(K); col.resize(K+1); user_k.resize(K);
     983                k = 0;
     984                for(j = 0; j < n; j++)
     985                {       for(i = 0; i < m; i++)
     986                        {       if( p[ i * n + j ] )
     987                                {       row[k]      = i;
     988                                        col[k]      = j;
     989                                        user_k[k] = k;
     990                                        k++;
     991                                }
     992                        }
     993                }
     994                col[K] = n;
     995
     996                // convert transposed sparsity pattern to a sparse_pack object
     997                // so can fold vector of bools and vector of sets into same function
     998                sparse_pack s_transpose;
     999                bool transpose = true;
     1000                vec_bool_to_sparse_pack(s_transpose, p, m, n, transpose);
     1001       
     1002                // now we have folded this into the following case
     1003                SparseJacobianFor(x, s_transpose, J, work);
     1004        }
     1005        else
     1006        {       // reverse mode, rows are sorted
     1007                row.resize(K+1); col.resize(K); user_k.resize(K);
     1008                k = 0;
     1009                for(i = 0; i < m; i++)
     1010                {       for(j = 0; j < n; j++)
     1011                        {       if( p[ i * n + j ] )
     1012                                {       row[k]      = i;
     1013                                        col[k]      = j;
     1014                                        user_k[k] = k;
     1015                                        k++;
     1016                                }
     1017                        }
     1018                }
     1019                row[K] = m;
     1020
     1021                // convert the sparsity pattern to a sparse_pack object
     1022                // so can fold vector of bools and vector of sets into same function
     1023                sparse_pack s;
     1024                bool transpose = false;
     1025                vec_bool_to_sparse_pack(s, p, m, n, transpose);
     1026       
     1027                // now we have folded this into the following case
     1028                SparseJacobianRev(x, s, J, work);
     1029        }
    2261030        // initialize the return value
    2271031        for(i = 0; i < m; i++)
     
    2291033                        jac[i * n + j] = zero;
    2301034
    231         if( n <= m )
    232         {       // use forward mode ----------------------------------------
    233        
    234                 // initial coloring
    235                 SizeVector color(n);
    236                 for(j = 0; j < n; j++)
    237                         color[j] = j;
    238 
    239                 // See GreedyPartialD2Coloring Algorithm Section 3.6.2 of
    240                 // Graph Coloring in Optimization Revisited by
    241                 // Assefaw Gebremedhin, Fredrik Maane, Alex Pothen
    242                 VectorBool forbidden(n);
    243                 for(j = 0; j < n; j++)
    244                 {       // initial all colors as ok for this column
    245                         for(k = 0; k < n; k++)
    246                                 forbidden[k] = false;
    247                         // for each row that is connected to column j
    248                         for(i = 0; i < m; i++) if( p[i * n + j] )
    249                         {       // for each column that is connected to row i
    250                                 for(k = 0; k < n; k++)
    251                                         if( p[i * n + k] & (j != k) )   
    252                                                 forbidden[ color[k] ] = true;
    253                         }
    254                         k = 0;
    255                         while( forbidden[k] && k < n )
    256                         {       k++;
    257                                 CPPAD_ASSERT_UNKNOWN( k < n );
    258                         }
    259                         color[j] = k;
    260                 }
    261                 size_t n_color = 1;
    262                 for(k = 0; k < n; k++)
    263                         n_color = std::max(n_color, color[k] + 1);
    264 
    265                 // direction vector for calls to forward
    266                 VectorBase dx(n);
    267 
    268                 // location for return values from Reverse
    269                 VectorBase dy(m);
    270 
    271                 // loop over colors
    272                 size_t c;
    273                 for(c = 0; c < n_color; c++)
    274                 {       // determine all the colums with this color
    275                         for(j = 0; j < n; j++)
    276                         {       if( color[j] == c )
    277                                         dx[j] = one;
    278                                 else    dx[j] = zero;
    279                         }
    280                         // call forward mode for all these columns at once
    281                         dy = Forward(1, dx);
    282 
    283                         // set the corresponding components of the result
    284                         for(j = 0; j < n; j++) if( color[j] == c )
    285                         {       for(i = 0; i < m; i++)
    286                                         if( p[ i * n + j ] )
    287                                                 jac[i * n + j] = dy[i];
    288                         }
    289                 }
    290         }
    291         else
    292         {       // use reverse mode ----------------------------------------
    293        
    294                 // initial coloring
    295                 SizeVector color(m);
    296                 for(i = 0; i < m; i++)
    297                         color[i] = i;
    298 
    299                 // See GreedyPartialD2Coloring Algorithm Section 3.6.2 of
    300                 // Graph Coloring in Optimization Revisited by
    301                 // Assefaw Gebremedhin, Fredrik Maane, Alex Pothen
    302                 VectorBool forbidden(m);
    303                 for(i = 0; i < m; i++)
    304                 {       // initial all colors as ok for this row
    305                         for(k = 0; k < m; k++)
    306                                 forbidden[k] = false;
    307                         // for each column that is connected to row i
    308                         for(j = 0; j < n; j++) if( p[i * n + j] )
    309                         {       // for each row that is connected to column j
    310                                 for(k = 0; k < m; k++)
    311                                         if( p[k * n + j] & (i != k) )   
    312                                                 forbidden[ color[k] ] = true;
    313                         }
    314                         k = 0;
    315                         while( forbidden[k] && k < m )
    316                         {       k++;
    317                                 CPPAD_ASSERT_UNKNOWN( k < n );
    318                         }
    319                         color[i] = k;
    320                 }
    321                 size_t n_color = 1;
    322                 for(k = 0; k < m; k++)
    323                         n_color = std::max(n_color, color[k] + 1);
    324 
    325                 // weight vector for calls to reverse
    326                 VectorBase w(m);
    327 
    328                 // location for return values from Reverse
    329                 VectorBase dw(n);
    330 
    331                 // loop over colors
    332                 size_t c;
    333                 for(c = 0; c < n_color; c++)
    334                 {       // determine all the rows with this color
    335                         for(i = 0; i < m; i++)
    336                         {       if( color[i] == c )
    337                                         w[i] = one;
    338                                 else    w[i] = zero;
    339                         }
    340                         // call reverse mode for all these rows at once
    341                         dw = Reverse(1, w);
    342 
    343                         // set the corresponding components of the result
    344                         for(i = 0; i < m; i++) if( color[i] == c )
    345                         {       for(j = 0; j < n; j++)
    346                                         if( p[ i * n + j ] )
    347                                                 jac[i * n + j] = dw[j];
    348                         }
    349                 }
    350         }
     1035        // now set the non-zero return values
     1036        for(k = 0; k < K; k++)
     1037                jac[row[k] * n + col[k]] = J[k];
    3511038}
    3521039
    3531040/*!
    354 Private helper function for SparseJacobian(x, p).
    355 
    356 All of the description in the public member function SparseJacobian(x, p)
    357 applies.
     1041Private helper function for vector of std::set<size_t> sparsity pattern cases.
    3581042
    3591043\param set_type
    3601044is a \c std::set<size_t> value.
    361 This argument is used to dispatch to the proper souce
     1045This argument is used to dispatch to the proper source
    3621046code depending on the vlaue of \c VectorSet::value_type.
    3631047
     
    3891073        size_t n = Domain();
    3901074
    391         // some values
     1075        // the value zero
    3921076        const Base zero(0);
    393         const Base one(1);
    3941077
    3951078        // check VectorSet is Simple Vector class with sets for elements
     
    4021085
    4031086        CPPAD_ASSERT_KNOWN(
    404                 x.size() == n,
    405                 "SparseJacobian: size of x not equal domain dimension for f"
    406         );
    407         CPPAD_ASSERT_KNOWN(
    4081087                p.size() == m,
    409                 "SparseJacobian: using sets and size of p "
     1088                "SparseJacobian: using size_t sets for sparsity pattern and p.size() "
    4101089                "not equal range dimension for f"
    4111090        );
     1091        CPPAD_ASSERT_UNKNOWN( x.size() == n );
    4121092        CPPAD_ASSERT_UNKNOWN(jac.size() == m * n);
    4131093
    414         // point at which we are evaluating the Jacobian
    415         Forward(0, x);
    416 
     1094        // row index, column index, value representation
     1095        // (used to fold problem into user interface case).
     1096        size_t K = 0;
     1097        std::set<size_t>::const_iterator itr;
     1098        for(i = 0; i < m; i++)
     1099        {       itr = p[i].begin();
     1100                while( itr != p[i].end() )
     1101                {       itr++;
     1102                        K++;
     1103                }
     1104        }       
     1105        VectorBase J(K);
     1106        sparse_jacobian_work work;
     1107        CppAD::vector<size_t>& row(work.r_sort);
     1108        CppAD::vector<size_t>& col(work.c_sort);
     1109        CppAD::vector<size_t>& user_k(work.k_sort);
     1110
     1111        if( n <= m )
     1112        {       // forward mode, columns are sorted
     1113                row.resize(K); col.resize(K+1); user_k.resize(K);
     1114
     1115                // need a pack_set copies
     1116                sparse_set s_transpose;
     1117                bool transpose = true;
     1118                vec_set_to_sparse_set(s_transpose, p, m, n, transpose);
     1119
     1120                k = 0;
     1121                for(j = 0; j < n; j++)
     1122                {       s_transpose.begin(j);
     1123                        i = s_transpose.next_element();
     1124                        while( i != s_transpose.end() )
     1125                        {       row[k]      = i;
     1126                                col[k]      = j;
     1127                                user_k[k] = k;
     1128                                k++;
     1129                                i    = s_transpose.next_element();
     1130                        }
     1131                }
     1132                col[K] = n;
     1133       
     1134                // now we have folded this into the following case
     1135                SparseJacobianFor(x, s_transpose, J, work);
     1136        }
     1137        else
     1138        {       // reverse mode, rows are sorted
     1139                row.resize(K+1); col.resize(K); user_k.resize(K);
     1140                k = 0;
     1141                for(i = 0; i < m; i++)
     1142                {       itr = p[i].begin();
     1143                        while( itr != p[i].end() )
     1144                        {       j = *itr++;
     1145                                row[k]      = i;
     1146                                col[k]      = j;
     1147                                user_k[k] = k;
     1148                                k++;
     1149                        }
     1150                }
     1151                row[K] = m;
     1152
     1153                // need a pack_set copy of sparsity pattern
     1154                sparse_set s;
     1155                bool transpose = false;
     1156                vec_set_to_sparse_set(s, p, m, n, transpose);
     1157
     1158                // now we have folded this into the following case
     1159                SparseJacobianRev(x, s, J, work);
     1160        }
    4171161        // initialize the return value
    4181162        for(i = 0; i < m; i++)
     
    4201164                        jac[i * n + j] = zero;
    4211165
    422         // create a copy of the transpose sparsity pattern
    423         VectorSet q(n);
    424         std::set<size_t>::const_iterator itr_i, itr_j;
    425         for(i = 0; i < m; i++)
    426         {       itr_j = p[i].begin();
    427                 while( itr_j != p[i].end() )
    428                 {       j = *itr_j++;
    429                         q[j].insert(i);
    430                 }
    431         }       
    432 
    433         if( n <= m )
    434         {       // use forward mode ----------------------------------------
    435        
    436                 // initial coloring
    437                 SizeVector color(n);
    438                 for(j = 0; j < n; j++)
    439                         color[j] = j;
    440 
    441                 // See GreedyPartialD2Coloring Algorithm Section 3.6.2 of
    442                 // Graph Coloring in Optimization Revisited by
    443                 // Assefaw Gebremedhin, Fredrik Maane, Alex Pothen
    444                 VectorBool forbidden(n);
    445                 for(j = 0; j < n; j++)
    446                 {       // initial all colors as ok for this column
    447                         for(k = 0; k < n; k++)
    448                                 forbidden[k] = false;
    449 
    450                         // for each row connected to column j
    451                         itr_i = q[j].begin();
    452                         while( itr_i != q[j].end() )
    453                         {       i = *itr_i++;
    454                                 // for each column connected to row i
    455                                 itr_j = p[i].begin();
    456                                 while( itr_j != p[i].end() )
    457                                 {       // if this is not j, forbid it
    458                                         k = *itr_j++;
    459                                         forbidden[ color[k] ] = (k != j);
    460                                 }
    461                         }
    462                         k = 0;
    463                         while( forbidden[k] && k < n )
    464                         {       k++;
    465                                 CPPAD_ASSERT_UNKNOWN( k < n );
    466                         }
    467                         color[j] = k;
    468                 }
    469                 size_t n_color = 1;
    470                 for(k = 0; k < n; k++)
    471                         n_color = std::max(n_color, color[k] + 1);
    472 
    473                 // direction vector for calls to forward
    474                 VectorBase dx(n);
    475 
    476                 // location for return values from Reverse
    477                 VectorBase dy(m);
    478 
    479                 // loop over colors
    480                 size_t c;
    481                 for(c = 0; c < n_color; c++)
    482                 {       // determine all the colums with this color
    483                         for(j = 0; j < n; j++)
    484                         {       if( color[j] == c )
    485                                         dx[j] = one;
    486                                 else    dx[j] = zero;
    487                         }
    488                         // call forward mode for all these columns at once
    489                         dy = Forward(1, dx);
    490 
    491                         // set the corresponding components of the result
    492                         for(j = 0; j < n; j++) if( color[j] == c )
    493                         {       itr_i = q[j].begin();
    494                                 while( itr_i != q[j].end() )
    495                                 {       i = *itr_i++;
    496                                         jac[i * n + j] = dy[i];
    497                                 }
    498                         }
    499                 }
    500         }
    501         else
    502         {       // use reverse mode ----------------------------------------
    503        
    504                 // initial coloring
    505                 SizeVector color(m);
    506                 for(i = 0; i < m; i++)
    507                         color[i] = i;
    508 
    509                 // See GreedyPartialD2Coloring Algorithm Section 3.6.2 of
    510                 // Graph Coloring in Optimization Revisited by
    511                 // Assefaw Gebremedhin, Fredrik Maane, Alex Pothen
    512                 VectorBool forbidden(m);
    513                 for(i = 0; i < m; i++)
    514                 {       // initial all colors as ok for this row
    515                         for(k = 0; k < m; k++)
    516                                 forbidden[k] = false;
    517 
    518                         // for each column connected to row i
    519                         itr_j = p[i].begin();
    520                         while( itr_j != p[i].end() )
    521                         {       j = *itr_j++;   
    522                                 // for each row connected to column j
    523                                 itr_i = q[j].begin();
    524                                 while( itr_i != q[j].end() )
    525                                 {       // if this is not i, forbid it
    526                                         k = *itr_i++;
    527                                         forbidden[ color[k] ] = (k != i);
    528                                 }
    529                         }
    530                         k = 0;
    531                         while( forbidden[k] && k < m )
    532                         {       k++;
    533                                 CPPAD_ASSERT_UNKNOWN( k < n );
    534                         }
    535                         color[i] = k;
    536                 }
    537                 size_t n_color = 1;
    538                 for(k = 0; k < m; k++)
    539                         n_color = std::max(n_color, color[k] + 1);
    540 
    541                 // weight vector for calls to reverse
    542                 VectorBase w(m);
    543 
    544                 // location for return values from Reverse
    545                 VectorBase dw(n);
    546 
    547                 // loop over colors
    548                 size_t c;
    549                 for(c = 0; c < n_color; c++)
    550                 {       // determine all the rows with this color
    551                         for(i = 0; i < m; i++)
    552                         {       if( color[i] == c )
    553                                         w[i] = one;
    554                                 else    w[i] = zero;
    555                         }
    556                         // call reverse mode for all these rows at once
    557                         dw = Reverse(1, w);
    558 
    559                         // set the corresponding components of the result
    560                         for(i = 0; i < m; i++) if( color[i] == c )
    561                         {       itr_j = p[i].begin();
    562                                 while( itr_j != p[i].end() )
    563                                 {       j = *itr_j++;
    564                                         jac[i * n + j] = dw[j];
    565                                 }
    566                         }
    567                 }
    568         }
     1166        // now set the non-zero return values
     1167        for(k = 0; k < K; k++)
     1168                jac[row[k] * n + col[k]] = J[k];
    5691169}
    570 
     1170// ==========================================================================
     1171// Public Member functions
     1172// ==========================================================================
     1173/*!
     1174Compute user specified subset of a sparse Jacobian using forward mode.
     1175
     1176The C++ source code corresponding to this operation is
     1177\verbatim
     1178        SparceJacobianForward(x, p, row, col, jac, work)
     1179\endverbatim
     1180
     1181\tparam Base
     1182is the base type for the recording that is stored in this
     1183ADFun<Base object.
     1184
     1185\tparam VectorBase
     1186is a simple vector class with elements of type \a Base.
     1187
     1188\tparam VectorSet
     1189is a simple vector class with elements of type
     1190\c bool or \c std::set<size_t>.
     1191
     1192\tparam VectorSize
     1193is a simple vector class with elements of type \c size_t.
     1194
     1195\param x
     1196is a vector specifing the point at which to compute the Jacobian.
     1197
     1198\param p
     1199is the sparsity pattern for the Jacobian that we are calculating.
     1200
     1201\param row
     1202is the vector of row indices for the returned Jacobian values.
     1203
     1204\param col
     1205is the vector of columns indices for the returned Jacobian values.
     1206It must have the same size are r.
     1207
     1208\param jac
     1209is the vector of Jacobian values.
     1210It must have the same size are r.
     1211The return value <code>jac[k]</code> is the partial of the
     1212<code>row[k]</code> component of the function with respect
     1213the the <code>col[k]</code> of its argument.
     1214
     1215\param work
     1216contains information that depends on the function object, sparsity pattern,
     1217\c row, and \c col vector.
     1218If these values are the same, \c work does not need to be recomputed.
     1219To be more specific,
     1220\c r_sort is sorted copy of \c row ,
     1221\c c_sort is sorted copy of \c col ,
     1222<code>k_sort[k]</code> is the original index corresponding to the
     1223values <code>r_sort[k]</code> and <code>c_sort[k]</code>.
     1224The order for the sort is either by rows, for reverse mode calculations,
     1225and by columns, for forward mode calculations.
     1226Let \c m be the range dimension, \c n the domain dimension,
     1227and \c K the size of \c row , \c col , and \c jac.
     1228If sorted by rows, there is one extra entry
     1229in the sorted row array and it has value <code>r_sort[K]=m</code>.
     1230If sorted by columns, there is one extra entry
     1231in the sorted column array and it has value <code>c_sort[K]=n</code>.
     1232The \c color vector is set and used by
     1233\c SparseJacobianFor, for the forward mode case,
     1234and by SparseJacobianRev, for the reverse mode case.
     1235
     1236\return
     1237Is the number of first order forward sweeps used to compute the
     1238requested Jacobian values. The total work, not counting the zero order
     1239forward sweep, or the time to combine computations, is proportional to this
     1240return value.
     1241*/
     1242template<class Base>
     1243template <class VectorBase, class VectorSet, class VectorSize>
     1244size_t ADFun<Base>::SparseJacobianForward(
     1245        const VectorBase&     x    ,
     1246        const VectorSet&      p    ,
     1247        const VectorSize&     row  ,
     1248        const VectorSize&     col  ,
     1249        VectorBase&           jac  ,
     1250        sparse_jacobian_work& work )
     1251{
     1252        size_t n = Domain();
     1253        size_t k, K = jac.size();
     1254        if( work.r_sort.size() == 0 )
     1255        {       // create version of (row, col,k) sorted by columns
     1256                size_t min_bytes = 3 * K * sizeof(size_t);
     1257                size_t cap_bytes;
     1258                void*  v_ptr  = thread_alloc::get_memory(min_bytes, cap_bytes);
     1259                size_t* crk  = reinterpret_cast<size_t*>(v_ptr);
     1260                for(k = 0; k < K; k++)
     1261                {       // must put column first so it is used for sorting
     1262                        crk[3 * k + 0] = col[k];
     1263                        crk[3 * k + 1] = row[k];
     1264                        crk[3 * k + 2] = k;
     1265                }
     1266                std::qsort(crk, K, 3 * sizeof(size_t), std_qsort_compare<size_t>);
     1267                work.c_sort.resize(K+1);
     1268                work.r_sort.resize(K);
     1269                work.k_sort.resize(K);
     1270                for(k = 0; k < K; k++)
     1271                {       work.c_sort[k] = crk[3 * k + 0]; 
     1272                        work.r_sort[k] = crk[3 * k + 1]; 
     1273                        work.k_sort[k] = crk[3 * k + 2]; 
     1274                }
     1275                work.c_sort[K] = n;
     1276                thread_alloc::return_memory(v_ptr);
     1277        }
     1278# ifndef NDEBUG
     1279        size_t m = Range();
     1280        CPPAD_ASSERT_KNOWN(
     1281                x.size() == n ,
     1282                "SparseJacobianForward: size of x not equal domain dimension for f."
     1283        );
     1284        CPPAD_ASSERT_KNOWN(
     1285                row.size() == K && col.size() == K ,
     1286                "SparseJacobianForward: either r or c does not have "
     1287                "the same size as jac."
     1288        );
     1289        CPPAD_ASSERT_KNOWN(
     1290                work.r_sort.size() == K   &&
     1291                work.c_sort.size() == K+1 &&
     1292                work.k_sort.size() == K   ,
     1293                "SparseJacobianForward: invalid value in work."
     1294        );
     1295        CPPAD_ASSERT_KNOWN(
     1296                work.color.size() == 0 || work.color.size() == n,
     1297                "SparseJacobianForward: invalid value in work."
     1298        );
     1299        for(k = 0; k < K; k++)
     1300        {       CPPAD_ASSERT_KNOWN(
     1301                        row[k] < m,
     1302                        "SparseJacobianForward: invalid value in r."
     1303                );
     1304                CPPAD_ASSERT_KNOWN(
     1305                        col[k] < n,
     1306                        "SparseJacobianForward: invalid value in c."
     1307                );
     1308                CPPAD_ASSERT_KNOWN(
     1309                        work.k_sort[k] < K,
     1310                        "SparseJacobianForward: invalid value in work."
     1311                );
     1312                CPPAD_ASSERT_KNOWN(
     1313                        work.r_sort[k] == row[ work.k_sort[k] ],
     1314                        "SparseJacobianForward: invalid value in work."
     1315                );
     1316                CPPAD_ASSERT_KNOWN(
     1317                        work.c_sort[k] == col[ work.k_sort[k] ],
     1318                        "SparseJacobianForward: invalid value in work."
     1319                );
     1320        }
     1321        if( work.color.size() != 0 )
     1322                for(size_t j = 0; j < n; j++) CPPAD_ASSERT_KNOWN(
     1323                        work.color[j] < n,
     1324                        "SparseJacobianForward: invalid value in work."
     1325        );
     1326# endif
     1327 
     1328        typedef typename VectorSet::value_type Set_type;
     1329        size_t n_sweep = SparseJacobianCase(Set_type(), x, p, jac, work);
     1330        return n_sweep;
     1331}
     1332/*!
     1333Compute user specified subset of a sparse Jacobian using reverse mode.
     1334
     1335The C++ source code corresponding to this operation is
     1336\verbatim
     1337        SparceJacobianReverse(x, p, row, col, jac, work)
     1338\endverbatim
     1339
     1340\tparam Base
     1341See \c SparseJacobianForward(x, p, row, col, jac, work).
     1342
     1343\tparam VectorBase
     1344See \c SparseJacobianForward(x, p, row, col, jac, work).
     1345
     1346\tparam VectorSet
     1347is either \c sparse_pack or \c sparse_set.
     1348
     1349\param x
     1350See \c SparseJacobianForward(x, p, row, col, jac, work).
     1351
     1352\param p
     1353Sparsity pattern for the Jacobian of this ADFun<Base> object.
     1354
     1355\param row
     1356See \c SparseJacobianForward(x, p, row, col, jac, work).
     1357
     1358\param col
     1359See \c SparseJacobianForward(x, p, row, col, jac, work).
     1360
     1361\param jac
     1362See \c SparseJacobianForward(x, p, row, col, jac, work).
     1363
     1364\param work
     1365See \c SparseJacobianForward(x, p, row, col, jac, work).
     1366
     1367\return
     1368See \c SparseJacobianForward(x, p, row, col, jac, work).
     1369*/
     1370template<class Base>
     1371template <class VectorBase, class VectorSet, class VectorSize>
     1372size_t ADFun<Base>::SparseJacobianReverse(
     1373        const VectorBase&     x    ,
     1374        const VectorSet&      p    ,
     1375        const VectorSize&     row  ,
     1376        const VectorSize&     col  ,
     1377        VectorBase&           jac  ,
     1378        sparse_jacobian_work& work )
     1379{
     1380        size_t m = Range();
     1381        size_t k, K = jac.size();
     1382        if( work.r_sort.size() == 0 )
     1383        {       // create version of (row, col,k) sorted by rows
     1384                size_t min_bytes = 3 * K * sizeof(size_t);
     1385                size_t cap_bytes;
     1386                void*  v_ptr  = thread_alloc::get_memory(min_bytes, cap_bytes);
     1387                size_t* rck  = reinterpret_cast<size_t*>(v_ptr);
     1388                for(k = 0; k < K; k++)
     1389                {       // must put column first so it is used for sorting
     1390                        rck[3 * k + 0] = row[k];
     1391                        rck[3 * k + 1] = col[k];
     1392                        rck[3 * k + 2] = k;
     1393                }
     1394                std::qsort(rck, K, 3 * sizeof(size_t), std_qsort_compare<size_t>);
     1395                work.c_sort.resize(K);
     1396                work.r_sort.resize(K+1);
     1397                work.k_sort.resize(K);
     1398                for(k = 0; k < K; k++)
     1399                {       work.r_sort[k] = rck[3 * k + 0]; 
     1400                        work.c_sort[k] = rck[3 * k + 1]; 
     1401                        work.k_sort[k] = rck[3 * k + 2]; 
     1402                }
     1403                work.r_sort[K] = m;
     1404                thread_alloc::return_memory(v_ptr);
     1405        }
     1406# ifndef NDEBUG
     1407        size_t n = Domain();
     1408        CPPAD_ASSERT_KNOWN(
     1409                x.size() == n ,
     1410                "SparseJacobianReverse: size of x not equal domain dimension for f."
     1411        );
     1412        CPPAD_ASSERT_KNOWN(
     1413                row.size() == K && col.size() == K ,
     1414                "SparseJacobianReverse: either r or c does not have "
     1415                "the same size as jac."
     1416        );
     1417        CPPAD_ASSERT_KNOWN(
     1418                work.r_sort.size() == K+1 &&
     1419                work.c_sort.size() == K   &&
     1420                work.k_sort.size() == K   ,
     1421                "SparseJacobianReverse: invalid value in work."
     1422        );
     1423        CPPAD_ASSERT_KNOWN(
     1424                work.color.size() == 0 || work.color.size() == m,
     1425                "SparseJacobianReverse: invalid value in work."
     1426        );
     1427        for(k = 0; k < K; k++)
     1428        {       CPPAD_ASSERT_KNOWN(
     1429                        row[k] < m,
     1430                        "SparseJacobianReverse: invalid value in r."
     1431                );
     1432                CPPAD_ASSERT_KNOWN(
     1433                        col[k] < n,
     1434                        "SparseJacobianReverse: invalid value in c."
     1435                );
     1436                CPPAD_ASSERT_KNOWN(
     1437                        work.k_sort[k] < K,
     1438                        "SparseJacobianReverse: invalid value in work."
     1439                );
     1440                CPPAD_ASSERT_KNOWN(
     1441                        work.r_sort[k] == row[ work.k_sort[k] ],
     1442                        "SparseJacobianReverse: invalid value in work."
     1443                );
     1444                CPPAD_ASSERT_KNOWN(
     1445                        work.c_sort[k] == col[ work.k_sort[k] ],
     1446                        "SparseJacobianReverse: invalid value in work."
     1447                );
     1448        }
     1449        if( work.color.size() != 0 )
     1450                for(size_t i = 0; i < m; i++) CPPAD_ASSERT_KNOWN(
     1451                        work.color[i] < m,
     1452                        "SparseJacobianReverse: invalid value in work."
     1453        );
     1454# endif
     1455 
     1456        typedef typename VectorSet::value_type Set_type;
     1457        size_t n_sweep = SparseJacobianCase(Set_type(), x, p, jac, work);
     1458        return n_sweep;
     1459}
    5711460/*!
    5721461Compute a sparse Jacobian.
     
    6071496        VectorBase jac(m * n);
    6081497
     1498        CPPAD_ASSERT_KNOWN(
     1499                x.size() == n,
     1500                "SparseJacobian: size of x not equal domain size for f."
     1501        );
     1502
    6091503        typedef typename VectorSet::value_type Set_type;
    6101504        SparseJacobianCase( Set_type(), x, p, jac);
     
    6461540        VectorBool p(m * n);
    6471541
    648         // return result
    649         VectorBase jac(m * n);
    650 
    6511542        if( n <= m )
    6521543        {       size_t j, k;
     
    6731564                p = RevSparseJac(m, s);
    6741565        }
    675         bool set_type = true; // only used to dispatch compiler to proper case
    676         SparseJacobianCase(set_type, x, p, jac);
    677         return jac;
     1566        return SparseJacobian(x, p);
    6781567}
    6791568
  • trunk/cppad/local/sparse_pack.hpp

    r2178 r2401  
    44
    55/* --------------------------------------------------------------------------
    6 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-11 Bradley M. Bell
     6CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
    77
    88CppAD is distributed under multiple licenses. This distribution is under
     
    2525Vector of sets of postivie integers, each set stored as a packed boolean array.
    2626*/
    27 
    2827
    2928class sparse_pack {
     
    131130                Pack mask = one << k;
    132131                data_[ index * n_pack_ + j] |= mask;
     132        }
     133        // -----------------------------------------------------------------
     134        /*! Is an element of a set.
     135
     136        \param index
     137        is the index for this set in the vector of sets.
     138
     139        \param element
     140        is the element we are checking to see if it is in the set.
     141
     142        \par Checked Assertions
     143        \li index    < n_set_
     144        \li element  < end_
     145        */
     146        bool is_element(size_t index, size_t element)
     147        {       static Pack one(1);
     148                static Pack zero(0);
     149                CPPAD_ASSERT_UNKNOWN( index   < n_set_ );
     150                CPPAD_ASSERT_UNKNOWN( element < end_ );
     151                size_t j  = element / n_bit_;
     152                size_t k  = element - j * n_bit_;
     153                Pack mask = one << k;
     154                return (data_[ index * n_pack_ + j] & mask) != zero;
    133155        }
    134156        // -----------------------------------------------------------------
     
    320342};
    321343
     344/*!
     345Copy a sparsity pattern from a vector of bools into a sparse_pack object.
     346
     347\tparam VectorBool
     348is a simple vector with elements of type bool.
     349
     350\param sparsity
     351The input value of sparisty does not matter.
     352Upon return it contains the same sparsity pattern as \c vec_bool
     353(or the transposed sparsity pattern).
     354
     355\param vec_bool
     356sparsity pattern that we are placing \c sparsity.
     357
     358\param n_row
     359number of rows in the sparsity pattern in \c vec_bool.
     360
     361\param n_col
     362number of columns in the sparsity pattern in \c vec_bool.
     363
     364\param transpose
     365if true, the sparsity pattern in \c sparsity is the transpose
     366of the one in \c vec_bool.
     367Otherwise it is the same sparsity pattern.
     368*/
     369template<class VectorBool>
     370void vec_bool_to_sparse_pack(
     371        sparse_pack&       sparsity  ,
     372        const VectorBool&  vec_bool  ,
     373        size_t             n_row     ,
     374        size_t             n_col     ,
     375        bool               transpose )
     376{       CPPAD_ASSERT_UNKNOWN( n_row * n_col == vec_bool.size() );
     377        size_t i, j;
     378
     379        // transposed pattern case
     380        if( transpose )
     381        {       sparsity.resize(n_col, n_row);
     382                for(j = 0; j < n_col; j++)
     383                {       for(i = 0; i < n_row; i++)
     384                        {       if( vec_bool[ i * n_col + j ] )
     385                                        sparsity.add_element(j, i);
     386                        }
     387                }
     388                return;
     389        }
     390
     391        // same pattern case
     392        sparsity.resize(n_row, n_col);
     393        for(i = 0; i < n_row; i++)
     394        {       for(j = 0; j < n_col; j++)
     395                {       if( vec_bool[ i * n_col + j ] )
     396                                sparsity.add_element(i, j);
     397                }
     398        }
     399        return;
     400}
     401
    322402CPPAD_END_NAMESPACE
    323403# endif
  • trunk/cppad/local/sparse_set.hpp

    r2178 r2401  
    44
    55/* --------------------------------------------------------------------------
    6 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-11 Bradley M. Bell
     6CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
    77
    88CppAD is distributed under multiple licenses. This distribution is under
     
    119119        }
    120120        // -----------------------------------------------------------------
     121        /*! Is an element of a set.
     122
     123        \param index
     124        is the index for this set in the vector of sets.
     125
     126        \param element
     127        is the element we are adding to the set.
     128
     129        \par Checked Assertions
     130        \li index    < n_set_
     131        \li element  < end_
     132        */
     133        bool is_element(size_t index, size_t element)
     134        {       // This routine should use the std::insert operator
     135                // that cashes the iterator of previous insertion for when
     136                // insertions occur in order. We should speed test that this
     137                // actually makes things faster.
     138                CPPAD_ASSERT_UNKNOWN( index   < n_set_ );
     139                CPPAD_ASSERT_UNKNOWN( element < end_ );
     140                std::set<size_t>::iterator itr = data_[ index ].find( element );
     141                return itr != data_[index].end();
     142        }
     143        // -----------------------------------------------------------------
    121144        /*! Begin retrieving elements from one of the sets.
    122145       
     
    272295};
    273296
     297/*!
     298Copy a sparsity pattern from a vector of sets into a sparse_set object.
     299
     300\tparam VectorSet
     301is a simple vector with elements of type \c std::set<size_t>.
     302
     303\param sparsity
     304The input value of sparisty does not matter.
     305Upon return it contains the same sparsity pattern as \c vec_set
     306(or the transposed sparsity pattern).
     307
     308\param vec_set
     309sparsity pattern that we are placing \c sparsity.
     310
     311\param n_row
     312number of rows in the sparsity pattern in \c vec_set.
     313
     314\param n_col
     315number of columns in the sparsity pattern in \c vec_set.
     316
     317\param transpose
     318if true, the sparsity pattern in \c sparsity is the transpose
     319of the one in \c vec_bool.
     320Otherwise it is the same sparsity pattern.
     321*/
     322template<class VectorSet>
     323void vec_set_to_sparse_set(
     324        sparse_set&        sparsity  ,
     325        const VectorSet&   vec_set   ,
     326        size_t             n_row     ,
     327        size_t             n_col     ,
     328        bool               transpose )
     329{       CPPAD_ASSERT_UNKNOWN( n_row == vec_set.size() );
     330        size_t i, j;
     331        std::set<size_t>::const_iterator itr;
     332
     333        // transposed pattern case
     334        if( transpose )
     335        {       sparsity.resize(n_col, n_row);
     336                for(i = 0; i < n_row; i++)
     337                {       itr = vec_set[i].begin();
     338                        while(itr != vec_set[i].end())
     339                        {       j = *itr++;
     340                                CPPAD_ASSERT_UNKNOWN( j < n_col );
     341                                sparsity.add_element(j, i);
     342                        }
     343                }
     344                return;
     345        }
     346
     347        // same pattern case
     348        sparsity.resize(n_row, n_col);
     349        for(i = 0; i < n_row; i++)
     350        {       itr = vec_set[i].begin();
     351                while(itr != vec_set[i].end())
     352                {       j = *itr++;
     353                        CPPAD_ASSERT_UNKNOWN( j < n_col );
     354                        sparsity.add_element(i, j);
     355                }
     356        }
     357        return;
     358}
     359
    274360CPPAD_END_NAMESPACE
    275361# endif
  • trunk/cppad/speed/sparse_evaluate.hpp

    r2269 r2401  
    5454        \DD{f}{x[j[k]]}{x[j[k]]}
    5555\] $$
    56 for some \( k \) between zero and \( \ell-1 \).
     56for some $latex k $$ between zero and $latex \ell-1 $$.
    5757All the other terms of the Hessian are zero.
    5858
  • trunk/example/sparse_hessian.cpp

    r1558 r2401  
    11/* $Id$ */
    22/* --------------------------------------------------------------------------
    3 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-09 Bradley M. Bell
     3CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
    44
    55CppAD is distributed under multiple licenses. This distribution is under
     
    3737        using CppAD::AD;
    3838        using CppAD::NearEqual;
    39         size_t i, j, k;
     39        size_t i, j, k, ell;
     40        typedef CPPAD_TEST_VECTOR< AD<double> >               a_vector;
     41        typedef CPPAD_TEST_VECTOR<double>                     d_vector;
     42        typedef CPPAD_TEST_VECTOR<size_t>                     i_vector;
     43        typedef CPPAD_TEST_VECTOR<bool>                       b_vector;
     44        typedef CPPAD_TEST_VECTOR< std::set<size_t> >         s_vector;
     45        double eps = 10. * CppAD::epsilon<double>();
    4046
    4147        // domain space vector
    42         size_t n = 3;
    43         CPPAD_TEST_VECTOR< AD<double> >  X(n);
    44         for(i = 0; i < n; i++)
    45                 X[i] = AD<double> (0);
     48        size_t n = 12;  // must be greater than or equal 3; see n_sweep below
     49        a_vector a_x(n);
     50        for(j = 0; j < n; j++)
     51                a_x[j] = AD<double> (0);
    4652
    4753        // declare independent variables and starting recording
    48         CppAD::Independent(X);
     54        CppAD::Independent(a_x);
    4955
     56        // range space vector
    5057        size_t m = 1;
    51         CPPAD_TEST_VECTOR< AD<double> >  Y(m);
    52         Y[0] = X[0] * X[0] + X[0] * X[1] + X[1] * X[1] + X[2] * X[2];
     58        a_vector a_y(m);
     59        a_y[0] = a_x[0]*a_x[1];
     60        for(j = 0; j < n; j++)
     61                a_y[0] += a_x[j] * a_x[j] * a_x[j];
    5362
    54         // create f: X -> Y and stop tape recording
    55         CppAD::ADFun<double> f(X, Y);
     63        // create f: x -> y and stop tape recording
     64        // (without executing zero order forward calculation)
     65        CppAD::ADFun<double> f;
     66        f.Dependent(a_x, a_y);
    5667
    57         // new value for the independent variable vector
    58         CPPAD_TEST_VECTOR<double> x(n);
     68        // new value for the independent variable vector, and weighting vector
     69        d_vector w(m), x(n);
     70        for(j = 0; j < n; j++)
     71                x[j] = double(j);
     72        w[0] = 1.0;
     73
     74        // vector used to check the value of the hessian
     75        d_vector check(n * n);
     76        for(ell = 0; ell < n * n; ell++)
     77                check[ell] = 0.0;
     78        ell        = 0 * n + 1;
     79        check[ell] = 1.0;
     80        ell        = 1 * n + 0;
     81        check[ell] = 1.0 ;
     82        for(j = 0; j < n; j++)
     83        {       ell = j * n + j;
     84                check[ell] = 6.0 * x[j];
     85        }
     86
     87        // -------------------------------------------------------------------
     88        // second derivative of y[0] w.r.t x
     89        d_vector hes(n * n);
     90        hes = f.SparseHessian(x, w);
     91        for(ell = 0; ell < n * n; ell++)
     92                ok &=  NearEqual(w[0] * check[ell], hes[ell], eps, eps );
     93
     94        // --------------------------------------------------------------------
     95        // example using vectors of bools to compute sparsity pattern for Hessian
     96        b_vector r_bool(n * n);
    5997        for(i = 0; i < n; i++)
    60                 x[i] = double(i);
     98        {       for(j = 0; j < n; j++)
     99                        r_bool[i * n + j] = false;
     100                r_bool[i * n + i] = true;
     101        }
     102        f.ForSparseJac(n, r_bool);
     103        //
     104        b_vector s_bool(m);
     105        for(i = 0; i < m; i++)
     106                s_bool[i] = w[i] != 0;
     107        b_vector p_bool = f.RevSparseHes(n, s_bool);
    61108
    62         // second derivative of y[1]
    63         CPPAD_TEST_VECTOR<double> w(m);
    64         w[0] = 1.;
    65         CPPAD_TEST_VECTOR<double> h( n * n );
    66         h = f.SparseHessian(x, w);
    67         /*
    68             [ 2 1 0 ]
    69         h = [ 1 2 0 ]
    70             [ 0 0 2 ]
    71         */
    72         CPPAD_TEST_VECTOR<double> check(n * n);
    73         check[0] = 2.; check[1] = 1.; check[2] = 0.;
    74         check[3] = 1.; check[4] = 2.; check[5] = 0.;
    75         check[6] = 0.; check[7] = 0.; check[8] = 2.;
    76         for(k = 0; k < n * n; k++)
    77                 ok &=  NearEqual(check[k], h[k], 1e-10, 1e-10 );
     109        hes = f.SparseHessian(x, w, p_bool);
     110        for(ell = 0; ell < n * n; ell++)
     111                ok &=  NearEqual(w[0] * check[ell], hes[ell], eps, eps );
    78112
    79         // use vectors of bools to compute sparse hessian -------------------
    80         // determine the sparsity pattern p for Hessian of w^T F
    81         CPPAD_TEST_VECTOR<bool> r_bool(n * n);
    82         for(j = 0; j < n; j++)
    83         {       for(k = 0; k < n; k++)
    84                         r_bool[j * n + k] = false;
    85                 r_bool[j * n + j] = true;
    86         }
    87         f.ForSparseJac(n, r_bool);
    88         //
    89         CPPAD_TEST_VECTOR<bool> s_bool(m);
    90         for(i = 0; i < m; i++)
    91                 s_bool[i] = w[i] != 0;
    92         CPPAD_TEST_VECTOR<bool> p_bool = f.RevSparseHes(n, s_bool);
    93 
    94         // test passing sparsity pattern
    95         h = f.SparseHessian(x, w, p_bool);
    96         for(k = 0; k < n * n; k++)
    97                 ok &=  NearEqual(check[k], h[k], 1e-10, 1e-10 );
    98 
    99         // use vectors of sets to compute sparse hessian ------------------
    100         // determine the sparsity pattern p for Hessian of w^T F
    101         CPPAD_TEST_VECTOR< std::set<size_t> > r_set(n);
    102         for(j = 0; j < n; j++)
    103                 r_set[j].insert(j);
    104         f.ForSparseJac(n, r_set);
    105         //
    106         CPPAD_TEST_VECTOR< std::set<size_t> > s_set(1);
    107         for(i = 0; i < m; i++)
     113        // --------------------------------------------------------------------
     114        // example using vectors of sets to compute sparsity pattern for Hessian
     115        s_vector r_set(n);
     116        for(i = 0; i < n; i++)
     117                r_set[i].insert(i);
     118        f.ForSparseJac(n, r_set);
     119        //
     120        s_vector s_set(m);
     121        for(i = 0; i < m; i++)
    108122                if( w[i] != 0. )
    109123                        s_set[0].insert(i);
    110         CPPAD_TEST_VECTOR< std::set<size_t> > p_set = f.RevSparseHes(n, s_set);
     124        s_vector p_set = f.RevSparseHes(n, s_set);
    111125
    112         // test passing sparsity pattern
    113         h = f.SparseHessian(x, w, p_set);
    114         for(k = 0; k < n * n; k++)
    115                 ok &=  NearEqual(check[k], h[k], 1e-10, 1e-10 );
     126        // example passing sparsity pattern to SparseHessian
     127        hes = f.SparseHessian(x, w, p_set);
     128        for(ell = 0; ell < n * n; ell++)
     129                ok &=  NearEqual(w[0] * check[ell], hes[ell], eps, eps );
     130
     131        // --------------------------------------------------------------------
     132        // use row and column indices to specify upper triangle of
     133        // non-zero elements of Hessian
     134        size_t K = n + 1;
     135        i_vector row(K), col(K);
     136        hes.resize(K);
     137        k = 0;
     138        for(j = 0; j < n; j++)
     139        {       // diagonal of Hessian
     140                row[k] = j;
     141                col[k] = j;
     142                k++;
     143        }
     144        // only off diagonal non-zero elemenet in upper triangle
     145        row[k] = 0;
     146        col[k] = 1;
     147        k++;
     148        ok &= k == K;
     149        CppAD::sparse_hessian_work work;
     150
     151        // can use p_set or p_bool.
     152        size_t n_sweep = f.SparseHessian(x, w, p_set, row, col, hes, work);
     153        for(k = 0; k < K; k++)
     154        {       ell = row[k] * n + col[k];
     155                ok &=  NearEqual(w[0] * check[ell], hes[k], eps, eps );
     156        }
     157        ok &= n_sweep == 2;
     158
     159        // now recompute at a different x and w (using work from previous call
     160        w[0]       = 2.0;
     161        x[1]       = 0.5;
     162        ell        = 1 * n + 1;
     163        check[ell] = 6.0 * x[1];
     164        n_sweep    = f.SparseHessian(x, w, p_set, row, col, hes, work);
     165        for(k = 0; k < K; k++)
     166        {       ell = row[k] * n + col[k];
     167                ok &=  NearEqual(w[0] * check[ell], hes[k], eps, eps );
     168        }
     169        ok &= n_sweep == 2;
     170       
     171
    116172
    117173        return ok;
  • trunk/example/sparse_jacobian.cpp

    r1551 r2401  
    11/* $Id$ */
    22/* --------------------------------------------------------------------------
    3 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-09 Bradley M. Bell
     3CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
    44
    55CppAD is distributed under multiple licenses. This distribution is under
     
    3939        using CppAD::AD;
    4040        using CppAD::NearEqual;
    41         size_t i, j, k;
     41        typedef CPPAD_TEST_VECTOR< AD<double> > a_vector;
     42        typedef CPPAD_TEST_VECTOR<double>       d_vector;
     43        typedef CPPAD_TEST_VECTOR<size_t>       i_vector;
     44        size_t i, j, k, ell;
     45        double eps = 10. * CppAD::epsilon<double>();
    4246
    4347        // domain space vector
    4448        size_t n = 4;
    45         CPPAD_TEST_VECTOR< AD<double> >  X(n);
    46         for(j = 0; j < n; j++)
    47                 X[j] = AD<double> (0);
     49        a_vector  a_x(n);
     50        for(j = 0; j < n; j++)
     51                a_x[j] = AD<double> (0);
    4852
    4953        // declare independent variables and starting recording
    50         CppAD::Independent(X);
     54        CppAD::Independent(a_x);
    5155
    5256        size_t m = 3;
    53         CPPAD_TEST_VECTOR< AD<double> >  Y(m);
    54         Y[0] = X[0] + X[1];
    55         Y[1] = X[2] + X[3];
    56         Y[2] = X[0] + X[1] + X[2] + X[3] * X[3] / 2.;
    57 
    58         // create f: X -> Y and stop tape recording
    59         CppAD::ADFun<double> f(X, Y);
     57        a_vector  a_y(m);
     58        a_y[0] = a_x[0] + a_x[1];
     59        a_y[1] = a_x[2] + a_x[3];
     60        a_y[2] = a_x[0] + a_x[1] + a_x[2] + a_x[3] * a_x[3] / 2.;
     61
     62        // create f: x -> y and stop tape recording
     63        CppAD::ADFun<double> f(a_x, a_y);
    6064
    6165        // new value for the independent variable vector
    62         CPPAD_TEST_VECTOR<double> x(n);
     66        d_vector x(n);
    6367        for(j = 0; j < n; j++)
    6468                x[j] = double(j);
    6569
    6670        // Jacobian of y without sparsity pattern
    67         CPPAD_TEST_VECTOR<double> jac(m * n);
     71        d_vector jac(m * n);
    6872        jac = f.SparseJacobian(x);
    6973        /*
     
    7276              [ 1 1 1 x_3]
    7377        */
    74         CPPAD_TEST_VECTOR<double> check(m * n);
     78        d_vector check(m * n);
    7579        check[0] = 1.; check[1] = 1.; check[2]  = 0.; check[3]  = 0.;
    7680        check[4] = 0.; check[5] = 0.; check[6]  = 1.; check[7]  = 1.;
    7781        check[8] = 1.; check[9] = 1.; check[10] = 1.; check[11] = x[3];
    78         for(k = 0; k < 12; k++)
    79                 ok &=  NearEqual(check[k], jac[k], 1e-10, 1e-10 );
     82        for(ell = 0; ell < check.size(); ell++)
     83                ok &=  NearEqual(check[ell], jac[ell], eps, eps );
    8084
    8185        // using packed boolean sparsity patterns
    8286        CppAD::vectorBool s_b(m * m), p_b(m * n);
    8387        for(i = 0; i < m; i++)
    84         {       for(k = 0; k < m; k++)
    85                         s_b[i * m + k] = false;
     88        {       for(ell = 0; ell < m; ell++)
     89                        s_b[i * m + ell] = false;
    8690                s_b[i * m + i] = true;
    8791        }
    8892        p_b   = f.RevSparseJac(m, s_b);
    8993        jac   = f.SparseJacobian(x, p_b);
    90         for(k = 0; k < 12; k++)
    91                 ok &=  NearEqual(check[k], jac[k], 1e-10, 1e-10 );
     94        for(ell = 0; ell < check.size(); ell++)
     95                ok &=  NearEqual(check[ell], jac[ell], eps, eps );
    9296
    9397        // using vector of sets sparsity patterns
     
    97101        p_s   = f.RevSparseJac(m, s_s);
    98102        jac   = f.SparseJacobian(x, p_s);
    99         for(k = 0; k < 12; k++)
    100                 ok &=  NearEqual(check[k], jac[k], 1e-10, 1e-10 );
     103        for(ell = 0; ell < check.size(); ell++)
     104                ok &=  NearEqual(check[ell], jac[ell], eps, eps );
     105
     106        // using row and column indices to compute non-zero in rows 1 and 2
     107        // (skip row 0).
     108        size_t K = 6;
     109        i_vector row(K), col(K);
     110        jac.resize(K);
     111        k = 0;
     112        for(j = 0; j < n; j++)
     113        {       for(i = 1; i < m; i++)
     114                {       ell = i * n + j;
     115                        if( p_b[ell] )
     116                        {       ok &= check[ell] != 0.;
     117                                row[k] = i;
     118                                col[k] = j;
     119                                k++;
     120                        }
     121                }
     122        }
     123        ok &= k == K;
     124
     125        // empty work structure
     126        CppAD::sparse_jacobian_work work;
     127
     128        // could use p_b
     129        size_t n_sweep = f.SparseJacobianReverse(x, p_s, row, col, jac, work);
     130        for(k = 0; k < K; k++)
     131        {    ell = row[k] * n + col[k];
     132                ok &= NearEqual(check[ell], jac[k], eps, eps);
     133        }
     134        ok &= n_sweep == 2;
     135
     136        // now recompute at a different x value (using work from previous call)
     137        check[11] = x[3] = 10.;
     138        n_sweep = f.SparseJacobianReverse(x, p_s, row, col, jac, work);
     139        for(k = 0; k < K; k++)
     140        {    ell = row[k] * n + col[k];
     141                ok &= NearEqual(check[ell], jac[k], eps, eps);
     142        }
     143        ok &= n_sweep == 2;
    101144
    102145        return ok;
     
    107150        using CppAD::AD;
    108151        using CppAD::NearEqual;
    109         size_t j, k;
     152        typedef CPPAD_TEST_VECTOR< AD<double> > a_vector;
     153        typedef CPPAD_TEST_VECTOR<double>       d_vector;
     154        typedef CPPAD_TEST_VECTOR<size_t>       i_vector;
     155        size_t i, j, k, ell;
     156        double eps = 10. * CppAD::epsilon<double>();
    110157
    111158        // domain space vector
    112159        size_t n = 3;
    113         CPPAD_TEST_VECTOR< AD<double> >  X(n);
    114         for(j = 0; j < n; j++)
    115                 X[j] = AD<double> (0);
     160        a_vector  a_x(n);
     161        for(j = 0; j < n; j++)
     162                a_x[j] = AD<double> (0);
    116163
    117164        // declare independent variables and starting recording
    118         CppAD::Independent(X);
     165        CppAD::Independent(a_x);
    119166
    120167        size_t m = 4;
    121         CPPAD_TEST_VECTOR< AD<double> >  Y(m);
    122         Y[0] = X[0] + X[2];
    123         Y[1] = X[0] + X[2];
    124         Y[2] = X[1] + X[2];
    125         Y[3] = X[1] + X[2] * X[2] / 2.;
    126 
    127         // create f: X -> Y and stop tape recording
    128         CppAD::ADFun<double> f(X, Y);
     168        a_vector  a_y(m);
     169        a_y[0] = a_x[0] + a_x[2];
     170        a_y[1] = a_x[0] + a_x[2];
     171        a_y[2] = a_x[1] + a_x[2];
     172        a_y[3] = a_x[1] + a_x[2] * a_x[2] / 2.;
     173
     174        // create f: x -> y and stop tape recording
     175        CppAD::ADFun<double> f(a_x, a_y);
    129176
    130177        // new value for the independent variable vector
    131         CPPAD_TEST_VECTOR<double> x(n);
     178        d_vector x(n);
    132179        for(j = 0; j < n; j++)
    133180                x[j] = double(j);
    134181
    135182        // Jacobian of y without sparsity pattern
    136         CPPAD_TEST_VECTOR<double> jac(m * n);
     183        d_vector jac(m * n);
    137184        jac = f.SparseJacobian(x);
    138185        /*
     
    142189              [ 0 1 x_2 ]
    143190        */
    144         CPPAD_TEST_VECTOR<double> check(m * n);
     191        d_vector check(m * n);
    145192        check[0] = 1.; check[1]  = 0.; check[2]  = 1.;
    146193        check[3] = 1.; check[4]  = 0.; check[5]  = 1.;
    147194        check[6] = 0.; check[7]  = 1.; check[8]  = 1.;
    148195        check[9] = 0.; check[10] = 1.; check[11] = x[2];
    149         for(k = 0; k < 12; k++)
    150                 ok &=  NearEqual(check[k], jac[k], 1e-10, 1e-10 );
     196        for(ell = 0; ell < check.size(); ell++)
     197                ok &=  NearEqual(check[ell], jac[ell], eps, eps );
    151198
    152199        // test using packed boolean vectors for sparsity pattern
    153200        CppAD::vectorBool r_b(n * n), p_b(m * n);
    154201        for(j = 0; j < n; j++)
    155         {       for(k = 0; k < n; k++)
    156                         r_b[j * n + k] = false;
     202        {       for(ell = 0; ell < n; ell++)
     203                        r_b[j * n + ell] = false;
    157204                r_b[j * n + j] = true;
    158205        }
    159206        p_b = f.ForSparseJac(n, r_b);
    160207        jac = f.SparseJacobian(x, p_b);
    161         for(k = 0; k < 12; k++)
    162                 ok &=  NearEqual(check[k], jac[k], 1e-10, 1e-10 );
     208        for(ell = 0; ell < check.size(); ell++)
     209                ok &=  NearEqual(check[ell], jac[ell], eps, eps );
    163210
    164211        // test using vector of sets for sparsity pattern
     
    168215        p_s = f.ForSparseJac(n, r_s);
    169216        jac = f.SparseJacobian(x, p_s);
    170         for(k = 0; k < 12; k++)
    171                 ok &=  NearEqual(check[k], jac[k], 1e-10, 1e-10 );
     217        for(ell = 0; ell < check.size(); ell++)
     218                ok &=  NearEqual(check[ell], jac[ell], eps, eps );
     219
     220        // using row and column indices to compute non-zero elements excluding
     221        // row 0 and column 0.
     222        size_t K = 5;
     223        i_vector row(K), col(K);
     224        jac.resize(K);
     225        k = 0;
     226        for(i = 1; i < m; i++)
     227        {       for(j = 1; j < n; j++)
     228                {       ell = i * n + j;
     229                        if( p_b[ell] )
     230                        {       ok &= check[ell] != 0.;
     231                                row[k] = i;
     232                                col[k] = j;
     233                                k++;
     234                        }
     235                }
     236        }
     237        ok &= k == K;
     238
     239        // empty work structure
     240        CppAD::sparse_jacobian_work work;
     241
     242        // could use p_s
     243        size_t n_sweep = f.SparseJacobianForward(x, p_b, row, col, jac, work);
     244        for(k = 0; k < K; k++)
     245        {    ell = row[k] * n + col[k];
     246                ok &= NearEqual(check[ell], jac[k], eps, eps);
     247        }
     248        ok &= n_sweep == 2;
     249
     250        // now recompute at a different x value (using work from previous call)
     251        check[11] = x[2] = 10.;
     252        n_sweep = f.SparseJacobianForward(x, p_s, row, col, jac, work);
     253        for(k = 0; k < K; k++)
     254        {    ell = row[k] * n + col[k];
     255                ok &= NearEqual(check[ell], jac[k], eps, eps);
     256        }
     257        ok &= n_sweep == 2;
    172258
    173259        return ok;
  • trunk/omh/glossary.omh

    r2333 r2401  
    3131
    3232$head AD Function$$
    33 Given an $xref/ADFun/$$ object $italic f$$
     33Given an $cref ADFun$$ object $icode f$$
    3434there is a corresponding
    35 AD of $italic Base$$ $xref/glossary/Operation/Sequence/operation sequence/1/$$.
     35AD of $icode Base$$ $cref/operation sequence/glossary/Operation/Sequence/$$.
    3636This operation sequence
    3737defines a function
    3838$latex F : B^n \rightarrow B^m $$
    39 where $italic B$$ is the space corresponding to objects of type $italic Base$$.
     39where $icode B$$ is the space corresponding to objects of type $icode Base$$,
     40$icode n$$ is the size of the $cref/domain/seq_property/Domain/$$ space, and
     41$icode m$$ is the size of the $cref/range/seq_property/Range/$$ space.
    4042We refer to $latex F$$ as the AD function corresponding to
    41 the operation sequence and to the object $italic f$$.
    42 (See the $xref/FunCheck/Discussion/FunCheck discussion/$$ for
     43the operation sequence stored in the object $icode f$$.
     44(See the $cref/FunCheck discussion/FunCheck/Discussion/$$ for
    4345possible differences between $latex F(x)$$ and the algorithm that defined
    4446the operation sequence.)
    4547
    4648$head AD of Base$$
    47 An object is called an AD of $italic Base$$ object its type is
    48 either $syntax%AD<%Base%>%$$
     49An object is called an AD of $icode Base$$ object its type is
     50either $codei%AD<%Base%>%$$
    4951(see the default and copy $cref/constructors/ad_ctor/$$
    50 or $syntax%VecAD<%Base%>::reference%$$ (see $xref/VecAD/$$)
    51 for some $italic Base$$ type.
     52or $codei%VecAD<%Base%>::reference%$$ (see $cref VecAD$$)
     53for some $icode Base$$ type.
    5254
    5355$head AD Levels Above Base$$
    5456$index level, AD$$
    5557$index AD, level$$
    56 If $italic Base$$ is a type,
    57 the AD levels above $italic Base$$
     58If $icode Base$$ is a type,
     59the AD levels above $icode Base$$
    5860is the following sequence of types:
    59 $syntax%
     61$codei%
    6062        AD<%Base%> %,% AD< AD<%Base%> > %,% AD< AD< AD<%Base%> > > %,% %...%
    6163%$$
     
    6365$head Base Function$$
    6466A function $latex f : B \rightarrow B$$
    65 is referred to as a $italic Base$$ function,
    66 if $italic Base$$ is a C++ type that represent elements of
    67 the domain and range space of $italic f$$; i.e. elements of $latex B$$.
     67is referred to as a $icode Base$$ function,
     68if $icode Base$$ is a C++ type that represent elements of
     69the domain and range space of $icode f$$; i.e. elements of $latex B$$.
    6870
    6971$head Base Type$$
    70 If $italic x$$ is an $syntax%AD<%Base%>%$$ object,
    71 $italic Base$$ is referred to as the base type for $italic x$$.
     72If $icode x$$ is an $codei%AD<%Base%>%$$ object,
     73$icode Base$$ is referred to as the base type for $icode x$$.
    7274
    7375$head Elementary Vector$$
     
    8486
    8587$subhead Atomic$$
    86 An atomic $italic Type$$ operation is an operation that
    87 has a $italic Type$$ result and is not made up of other
     88An atomic $icode Type$$ operation is an operation that
     89has a $icode Type$$ result and is not made up of other
    8890more basic operations.
    8991
    9092$subhead Sequence$$
    91 A sequence of atomic $italic Type$$ operations
    92 is called a $italic Type$$ operation sequence.
    93 A sequence of atomic $xref/glossary/AD of Base/AD of Base/$$ operations
    94 is referred to as an AD of $italic Base$$ operation sequence.
     93A sequence of atomic $icode Type$$ operations
     94is called a $icode Type$$ operation sequence.
     95A sequence of atomic $cref/AD of Base/glossary/AD of Base/$$ operations
     96is referred to as an AD of $icode Base$$ operation sequence.
    9597The abbreviated notation AD operation sequence is often used
    9698when it is not necessary to specify the base type.
    9799
    98100$subhead Dependent$$
    99 Suppose that $italic x$$ and $italic y$$ are $italic Type$$ objects and
     101Suppose that $icode x$$ and $icode y$$ are $icode Type$$ objects and
    100102the result of
    101 $syntax%
     103$codei%
    102104        %x% < %y%
    103105%$$
    104 has type $code bool$$ (where $italic Type$$ is not the same as $code bool$$).
     106has type $code bool$$ (where $icode Type$$ is not the same as $code bool$$).
    105107If one executes the following code
    106 $syntax%
     108$codei%
    107109        if( %x% < %y% )
    108110                %y% = cos(%x%);
    109111        else    %y% = sin(%x%);
    110112%$$
    111 the choice above depends on the value of $italic x$$ and $italic y$$
    112 and the two choices result in a different $italic Type$$ operation sequence.
    113 In this case, we say that the $italic Type$$ operation sequence depends
    114 on $italic x$$ and $italic y$$.
     113the choice above depends on the value of $icode x$$ and $icode y$$
     114and the two choices result in a different $icode Type$$ operation sequence.
     115In this case, we say that the $icode Type$$ operation sequence depends
     116on $icode x$$ and $icode y$$.
    115117
    116118$subhead Independent$$
    117 Suppose that $italic i$$ and $italic n$$ are $code size_t$$ objects,
    118 and $syntax%%x%[%i%]%$$, $italic y$$ are $italic Type$$ objects,
    119 where $italic Type$$ is different from $code size_t$$.
    120 The $italic Type$$ sequence of operations corresponding to
    121 $syntax%
     119Suppose that $icode i$$ and $icode n$$ are $code size_t$$ objects,
     120and $icode%x%[%i%]%$$, $icode y$$ are $icode Type$$ objects,
     121where $icode Type$$ is different from $code size_t$$.
     122The $icode Type$$ sequence of operations corresponding to
     123$codei%
    122124        %y% = %Type%(0);
    123125        for(%i% = 0; %i% < %n%; %i%++)
    124126                %y% += %x%[%i%];
    125127%$$
    126 does not depend on the value of $italic x$$ or $italic y$$.
    127 In this case, we say that the $italic Type$$ operation sequence
    128 is independent of $italic y$$ and the elements of $italic x$$.
     128does not depend on the value of $icode x$$ or $icode y$$.
     129In this case, we say that the $icode Type$$ operation sequence
     130is independent of $icode y$$ and the elements of $icode x$$.
    129131
    130132$head Parameter$$
    131 All $italic Base$$ objects are parameters.
    132 An $syntax%AD<%Base%>%$$ object $italic u$$ is currently a parameter if
     133All $icode Base$$ objects are parameters.
     134An $codei%AD<%Base%>%$$ object $icode u$$ is currently a parameter if
    133135its value does not depend on the value of
    134 an $xref/Independent/$$ variable vector for an
     136an $cref Independent$$ variable vector for an
    135137$cref/active tape/glossary/Tape/Active/$$.
    136 If $italic u$$ is a parameter, the function
    137 $xref/ParVar//Parameter(u)/$$ returns true
    138 and $xref/ParVar//Variable(u)/$$ returns false.
     138If $icode u$$ is a parameter, the function
     139$cref/Parameter(u)/ParVar/$$ returns true
     140and $cref/Variable(u)/ParVar/$$ returns false.
    139141
    140142$head Sparsity Pattern$$
     
    177179B_{i * n + j} = {\rm true}
    178180\] $$
    179 Given two sparsity patterns $latex B$$ and $italic C$$
    180 for a matrix $italic A$$, we say that $italic B$$ is more efficient than
    181 $italic C$$ if $italic B$$ has fewer true elements than $italic C$$.
     181Given two sparsity patterns $latex B$$ and $icode C$$
     182for a matrix $icode A$$, we say that $icode B$$ is more efficient than
     183$icode C$$ if $icode B$$ has fewer true elements than $icode C$$.
    182184For example, if $latex A$$ is the identity matrix,
    183185$latex \[
     
    196198\; \Rightarrow \; j \in S_i
    197199\] $$
    198 Given two sparsity patterns $latex S$$ and $italic T$$
    199 for a matrix $italic A$$, we say that $italic S$$ is more efficient than
    200 $italic T$$ if $italic S$$ has fewer elements than $italic T$$.
     200Given two sparsity patterns $latex S$$ and $icode T$$
     201for a matrix $icode A$$, we say that $icode S$$ is more efficient than
     202$icode T$$ if $icode S$$ has fewer elements than $icode T$$.
    201203For example, if $latex A$$ is the identity matrix,
    202204$latex \[
     
    209211$subhead Active$$
    210212A new tape is created and becomes active
    211 after each call of the form (see $cref/Independent/$$)
    212 $syntax%
     213after each call of the form (see $cref Independent$$)
     214$codei%
    213215        Independent(%x%)
    214216%$$
    215 All operations that depend on the elements of $italic x$$ are
     217All operations that depend on the elements of $icode x$$ are
    216218recorded on this active tape.
    217219
     
    221223stored in a tape must be transferred to a function object using the syntax
    222224(see $cref/ADFun<Base> f(x, y)/FunConstruct/$$)
    223 $syntax%
     225$codei%
    224226        ADFun<%Base%> %f%( %x%, %y%)
    225227%$$
    226228or using the syntax (see $cref/f.Dependent(x, y)/Dependent/$$)
    227 $syntax%
     229$codei%
    228230        %f%.Dependent( %x%, %y%)
    229231%$$
     
    231233
    232234$subhead Independent Variable$$
    233 While the tape is active, we refer to the elements of $italic x$$
     235While the tape is active, we refer to the elements of $icode x$$
    234236as the independent variables for the tape.
    235237When the tape becomes inactive,
     
    265267
    266268$head Variable$$
    267 An $syntax%AD<%Base%>%$$ object $italic u$$ is a variable if
     269An $codei%AD<%Base%>%$$ object $icode u$$ is a variable if
    268270its value depends on an independent variable vector for
    269271a currently $cref/active tape/glossary/Tape/Active/$$.
    270 If $italic u$$ is a variable,
    271 $xref/ParVar//Variable(u)/$$ returns true and
    272 $xref/ParVar//Parameter(u)/$$ returns false.
     272If $icode u$$ is a variable,
     273$cref/Variable(u)/ParVar/$$ returns true and
     274$cref/Parameter(u)/ParVar/$$ returns false.
    273275For example,
    274276directly after the code sequence
    275 $syntax%
     277$codei%
    276278        Independent(%x%);
    277279        AD<double> %u% = %x%[0];
    278280%$$
    279 the $syntax%AD<double>%$$ object $italic u$$ is currently a variable.
     281the $codei%AD<double>%$$ object $icode u$$ is currently a variable.
    280282Directly after the code sequence
    281 $syntax%
     283$codei%
    282284        Independent(%x%);
    283285        AD<double> %u% = %x%[0];
    284286        %u% = 5;
    285287%$$
    286 $italic u$$  is currently a $xref/glossary/Parameter/parameter/$$
     288$icode u$$  is currently a $cref/parameter/glossary/Parameter/$$
    287289(not a variable).
    288290$pre
     
    290292$$
    291293Note that we often drop the word currently and
    292 just refer to an $syntax%AD<%Base%>%$$ object as a variable
     294just refer to an $codei%AD<%Base%>%$$ object as a variable
    293295or parameter.
    294296
  • trunk/omh/whats_new_12.omh

    r2354 r2401  
    4444The purpose of this section is to
    4545assist you in learning about changes between various versions of CppAD.
     46
     47$head 05-24$$
     48Merged in changes from $code branches/sparse$$:
     49$list number$$
     50A new interface was added to
     51$cref sparse_jacobian$$ and $cref sparse_hessian$$.
     52This interface
     53returns a sparse representation of the corresponding matrices
     54using row and column index vectors.
     55$lnext
     56The examples
     57$cref sparse_jacobian.cpp$$ and
     58$cref sparse_hessian.cpp$$  were improved
     59and extended to include the new interface.
     60$lnext
     61The definition of an
     62$cref/AD function/glossary/AD Function/$$ was improved
     63to include definition of the corresponding $icode n$$ and $icode m$$.
     64$lend
     65
    4666
    4767$head 04-19$$
  • trunk/speed/src/link_sparse_hessian.cpp

    r2268 r2401  
    5151        \DD{f}{x[j[k]]}{x[j[k]]}
    5252\] $$
    53 for some \( k \) between zero and \( \ell-1 \).
     53for some $latex k $$ between zero and $latex \ell-1 $$.
    5454All the other terms of the Hessian are zero.
    5555
  • trunk/test_more/sparse_hessian.cpp

    r1558 r2401  
    11/* $Id$ */
    22/* --------------------------------------------------------------------------
    3 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-09 Bradley M. Bell
     3CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
    44
    55CppAD is distributed under multiple licenses. This distribution is under
     
    1717# include <cppad/cppad.hpp>
    1818namespace { // ---------------------------------------------------------
     19
     20bool rc_tridiagonal(void)
     21{       bool ok = true;
     22        using CppAD::AD;
     23        using CppAD::NearEqual;
     24        size_t i, j, k, ell;
     25        double eps = 10. * CppAD::epsilon<double>();
     26
     27        size_t n = 12; // must be greater than or equal 3; see n_sweep.
     28        size_t m = n - 1;
     29        CPPAD_TEST_VECTOR< AD<double> > a_x(n), a_y(m);
     30        CPPAD_TEST_VECTOR<double> x(n), check(n * n), w(m);
     31        for(j = 0; j < n; j++)
     32                a_x[j] = x[j] = double(j+1);
     33
     34        // declare independent variables and start taping
     35        CppAD::Independent(a_x);
     36
     37        for(ell = 0; ell < n * n; ell++)
     38                check[ell] = 0.0;
     39
     40        for(i = 0; i < m; i++)
     41        {       AD<double> diff = a_x[i+1] - a_x[i];   
     42                a_y[i] = 0.5 * diff * diff / double(i+2);
     43                w[i] = double(i+1);
     44                ell         = i * n + i;
     45                check[ell] += w[i] / double(i+2);
     46                ell         = (i+1) * n + i+1;
     47                check[ell] += w[i] / double(i+2);
     48                ell         = (i+1) * n + i;
     49                check[ell] -= w[i] / double(i+2);
     50                ell         = i * n + i+1;
     51                check[ell] -= w[i] / double(i+2);
     52        }
     53
     54        // create f: x -> y
     55        CppAD::ADFun<double> f(a_x, a_y);
     56                 
     57        // determine the sparsity pattern p for Hessian of w^T f
     58        typedef CppAD::vector< std::set<size_t> > VectorSet;
     59        VectorSet p_r(n);
     60        for(j = 0; j < n; j++)
     61                p_r[j].insert(j);
     62        f.ForSparseJac(n, p_r);
     63        //
     64        VectorSet p_s(1);
     65        for(i = 0; i < m; i++)
     66                if( w[i] != 0 )
     67                        p_s[0].insert(i);
     68        VectorSet p_h = f.RevSparseHes(n, p_s);
     69
     70        // requres the upper triangle of the Hessian
     71        size_t K = 2 * n - 1;
     72        CPPAD_TEST_VECTOR<size_t> r(K), c(K);
     73        CPPAD_TEST_VECTOR<double> hes(K);
     74        k = 0;
     75        for(i = 0; i < n; i++)
     76        {       r[k] = i;
     77                c[k] = i;
     78                k++;
     79                if( i < n-1 )
     80                {       r[k] = i;
     81                        c[k] = i+1;
     82                        k++;
     83                }
     84        }
     85        ok &= k == K;
     86
     87        // test computing sparse Hessian
     88        CppAD::sparse_hessian_work work;
     89        size_t n_sweep = f.SparseHessian(x, w, p_h, r, c, hes, work);
     90        ok &= n_sweep == 3;
     91        for(k = 0; k < K; k++)
     92        {       ell = r[k] * n + c[k];
     93                ok &=  NearEqual(check[ell], hes[k], eps, eps);
     94                ell = c[k] * n + r[k];
     95                ok &=  NearEqual(check[ell], hes[k], eps, eps);
     96        }
     97
     98        return ok;
     99}
     100
     101
    19102template <class VectorBase, class VectorBool>
    20103bool bool_case()
     
    63146
    64147        // determine the sparsity pattern p for Hessian of w^T F
    65         VectorBool r(n * n);
    66         for(j = 0; j < n; j++)
    67         {       for(k = 0; k < n; k++)
    68                         r[j * n + k] = false;
    69                 r[j * n + j] = true;
    70         }
    71         f.ForSparseJac(n, r);
    72         //
    73         VectorBool s(m);
    74         for(i = 0; i < m; i++)
    75                 s[i] = w[i] != 0;
    76         VectorBool p = f.RevSparseHes(n, s);
     148        VectorBool r(n * n);
     149        for(j = 0; j < n; j++)
     150        {       for(k = 0; k < n; k++)
     151                r[j * n + k] = false;
     152                r[j * n + j] = true;
     153        }
     154        f.ForSparseJac(n, r);
     155        //
     156        VectorBool s(m);
     157        for(i = 0; i < m; i++)
     158                s[i] = w[i] != 0;
     159        VectorBool p = f.RevSparseHes(n, s);
    77160
    78161        // test passing sparsity pattern
     
    129212
    130213        // determine the sparsity pattern p for Hessian of w^T F
    131         VectorSet r(n);
    132         for(j = 0; j < n; j++)
     214        VectorSet r(n);
     215        for(j = 0; j < n; j++)
    133216                r[j].insert(j);
    134         f.ForSparseJac(n, r);
    135         //
    136         VectorSet s(1);
    137         for(i = 0; i < m; i++)
    138                 if( w[i] != 0 )
     217        f.ForSparseJac(n, r);
     218        //
     219        VectorSet s(1);
     220        for(i = 0; i < m; i++)
     221                if( w[i] != 0 )
    139222                        s[0].insert(i);
    140         VectorSet p = f.RevSparseHes(n, s);
     223        VectorSet p = f.RevSparseHes(n, s);
    141224
    142225        // test passing sparsity pattern
     
    152235bool sparse_hessian(void)
    153236{       bool ok = true;
     237
     238        ok &= rc_tridiagonal();
     239        // ---------------------------------------------------------------
    154240        // vector of bool cases
    155241        ok &= bool_case< CppAD::vector  <double>, CppAD::vectorBool   >();
  • trunk/test_more/sparse_jacobian.cpp

    r1658 r2401  
    11/* $Id$ */
    22/* --------------------------------------------------------------------------
    3 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-10 Bradley M. Bell
     3CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
    44
    55CppAD is distributed under multiple licenses. This distribution is under
     
    1818namespace { // ---------------------------------------------------------
    1919
     20bool rc_tridiagonal(void)
     21{       bool ok = true;
     22        using CppAD::AD;
     23        using CppAD::NearEqual;
     24        size_t i, j, k, ell;
     25        double eps = 10. * CppAD::epsilon<double>();
     26
     27        // domain space vector
     28        size_t n = 13; // must be greater than or equal 3 (see n_sweep below)
     29        CPPAD_TEST_VECTOR< AD<double> >  X(n);
     30        CPPAD_TEST_VECTOR<double>        x(n);
     31        for(j = 0; j < n; j++)
     32                X[j] = x[j] = double(j+1);
     33
     34        // declare independent variables and starting recording
     35        CppAD::Independent(X);
     36
     37        size_t m = n;
     38        CPPAD_TEST_VECTOR< AD<double> >  Y(m);
     39        CPPAD_TEST_VECTOR<double> check(m * n );
     40        for(ell = 0; ell < m * n; ell++)
     41                check[ell] = 0.0;
     42
     43        size_t K = 0;
     44        for(i = 0; i < n; i++)
     45        {       ell        = i * n + i;
     46                Y[i]       = (ell+1) * 0.5 * X[i] * X[i];
     47                check[ell] = (ell+1) * x[i];
     48                K++;
     49                if( i < n-1 )
     50                {       j          = i + 1;
     51                        ell        = i * n + j;
     52                        Y[i]      += (ell+1) * 0.5 * X[i+1] * X[i+1];
     53                        check[ell] = (ell+1) * x[i+1];
     54                        K++;
     55                }
     56                if(i > 0 )
     57                {       j          = i - 1;
     58                        ell        = i * n + j;
     59                        Y[i]      += (ell+1) * 0.5 * X[i-1] * X[i-1];
     60                        check[ell] = (ell+1) * x[i-1];
     61                }
     62        }
     63
     64        // create f: X -> Y and stop tape recording
     65        CppAD::ADFun<double> f(X, Y);
     66
     67        // sparsity pattern
     68        CppAD::vector< std::set<size_t> > s(m), p(m);
     69        for(i = 0; i < m; i++)
     70                s[i].insert(i);
     71        p   = f.RevSparseJac(m, s);
     72
     73        // Request the upper triangle of the array
     74        CPPAD_TEST_VECTOR<size_t> r(K), c(K);
     75        CPPAD_TEST_VECTOR<double> jac(K);
     76        k = 0;
     77        for(i = 0; i < n; i++)
     78        {       r[k] = i;
     79                c[k] = i;
     80                k++;   
     81                if( i < n-1 )
     82                {       r[k] = i;
     83                        c[k] = i+1;
     84                        k++;
     85                }
     86        }
     87        ok &= K == k;
     88
     89        CppAD::sparse_jacobian_work work;
     90        size_t n_sweep = f.SparseJacobianForward(x, p, r, c, jac, work);
     91        ok &= n_sweep == 3;
     92        for(k = 0; k < K; k++)
     93        {       ell = r[k] * n + c[k];
     94                ok &=  NearEqual(check[ell], jac[k], eps, eps);
     95        }
     96        work.clear();
     97        n_sweep = f.SparseJacobianReverse(x, p, r, c, jac, work);
     98        ok &= n_sweep == 3;
     99        for(k = 0; k < K; k++)
     100        {       ell = r[k] * n + c[k];
     101                ok &=  NearEqual(check[ell], jac[k], eps, eps);
     102        }
     103
     104        return ok;
     105}
     106
     107template <class VectorBase, class VectorSet>
     108bool rc_set(void)
     109{       bool ok = true;
     110        using CppAD::AD;
     111        using CppAD::NearEqual;
     112        size_t i, j, k, ell;
     113        double eps = 10. * CppAD::epsilon<double>();
     114
     115        // domain space vector
     116        size_t n = 4;
     117        CPPAD_TEST_VECTOR< AD<double> >  X(n);
     118        for(j = 0; j < n; j++)
     119                X[j] = AD<double> (0);
     120
     121        // declare independent variables and starting recording
     122        CppAD::Independent(X);
     123
     124        size_t m = 3;
     125        CPPAD_TEST_VECTOR< AD<double> >  Y(m);
     126        Y[0] = 1.0*X[0] + 2.0*X[1];
     127        Y[1] = 3.0*X[2] + 4.0*X[3];
     128        Y[2] = 5.0*X[0] + 6.0*X[1] + 7.0*X[3]*X[3]/2.;
     129
     130        // create f: X -> Y and stop tape recording
     131        CppAD::ADFun<double> f(X, Y);
     132
     133        // new value for the independent variable vector
     134        VectorBase x(n);
     135        for(j = 0; j < n; j++)
     136                x[j] = double(j);
     137
     138        // Jacobian of y
     139        /*
     140              [ 1 2 0 0    ]
     141        jac = [ 0 0 3 4    ]
     142              [ 5 6 0 7*x_3]
     143        */
     144        VectorBase check(m * n);
     145        check[0] = 1.; check[1] = 2.; check[2]  = 0.; check[3]  = 0.;
     146        check[4] = 0.; check[5] = 0.; check[6]  = 3.; check[7]  = 4.;
     147        check[8] = 5.; check[9] = 6.; check[10] = 0.; check[11] = 7.*x[3];
     148
     149        // sparsity pattern
     150        VectorSet s(m), p(m);
     151        for(i = 0; i < m; i++)
     152                s[i].insert(i);
     153        p   = f.RevSparseJac(m, s);
     154
     155        // Use forward mode to compute columns 0 and 2
     156        // (make sure order of rows and columns does not matter)
     157        CPPAD_TEST_VECTOR<size_t> r(3), c(3);
     158        VectorBase jac(3);
     159        r[0] = 2; c[0] = 0;
     160        r[1] = 1; c[1] = 2;
     161        r[2] = 0; c[2] = 0;
     162        CppAD::sparse_jacobian_work work;
     163        size_t n_sweep = f.SparseJacobianForward(x, p, r, c, jac, work);
     164        for(k = 0; k < 3; k++)
     165        {       ell = r[k] * n + c[k];
     166                ok &=  NearEqual(check[ell], jac[k], eps, eps);
     167        }
     168        ok &= (n_sweep == 1);
     169
     170        // Use reverse mode to compute rows 0 and 1
     171        // (make sure order of rows and columns does not matter)
     172        r.resize(4), c.resize(4); jac.resize(4);
     173        r[0] = 0; c[0] = 0;
     174        r[1] = 1; c[1] = 2;
     175        r[2] = 0; c[2] = 1;
     176        r[3] = 1; c[3] = 3;
     177        work.clear();
     178        n_sweep = f.SparseJacobianReverse(x, p, r, c, jac, work);
     179        for(k = 0; k < 4; k++)
     180        {       ell = r[k] * n + c[k];
     181                ok &=  NearEqual(check[ell], jac[k], eps, eps);
     182        }
     183        ok &= (n_sweep == 1);
     184
     185        return ok;
     186}
    20187template <class VectorBase, class VectorBool>
    21 bool reverse_bool()
     188bool rc_bool(void)
     189{       bool ok = true;
     190        using CppAD::AD;
     191        using CppAD::NearEqual;
     192        size_t j, k, ell;
     193        double eps = 10. * CppAD::epsilon<double>();
     194
     195        // domain space vector
     196        size_t n = 4;
     197        CPPAD_TEST_VECTOR< AD<double> >  X(n);
     198        for(j = 0; j < n; j++)
     199                X[j] = AD<double> (0);
     200
     201        // declare independent variables and starting recording
     202        CppAD::Independent(X);
     203
     204        size_t m = 3;
     205        CPPAD_TEST_VECTOR< AD<double> >  Y(m);
     206        Y[0] = 1.0*X[0] + 2.0*X[1];
     207        Y[1] = 3.0*X[2] + 4.0*X[3];
     208        Y[2] = 5.0*X[0] + 6.0*X[1] + 7.0*X[3]*X[3]/2.;
     209
     210        // create f: X -> Y and stop tape recording
     211        CppAD::ADFun<double> f(X, Y);
     212
     213        // new value for the independent variable vector
     214        VectorBase x(n);
     215        for(j = 0; j < n; j++)
     216                x[j] = double(j);
     217
     218        // Jacobian of y
     219        /*
     220              [ 1 2 0 0    ]
     221        jac = [ 0 0 3 4    ]
     222              [ 5 6 0 7*x_3]
     223        */
     224        VectorBase check(m * n);
     225        check[0] = 1.; check[1] = 2.; check[2]  = 0.; check[3]  = 0.;
     226        check[4] = 0.; check[5] = 0.; check[6]  = 3.; check[7]  = 4.;
     227        check[8] = 5.; check[9] = 6.; check[10] = 0.; check[11] = 7.*x[3];
     228        VectorBool s(m * n);
     229        s[0] = true;   s[1] = true;   s[2] = false;   s[3] = false;
     230        s[4] = false;  s[5] = false;  s[6] = true;    s[7] = true;
     231        s[8] = true;   s[9] = true;  s[10] = false;  s[11] = true;
     232
     233        // Use forward mode to compute columns 0 and 2
     234        // (make sure order of rows and columns does not matter)
     235        CPPAD_TEST_VECTOR<size_t> r(3), c(3);
     236        VectorBase jac(3);
     237        r[0] = 2; c[0] = 0;
     238        r[1] = 1; c[1] = 2;
     239        r[2] = 0; c[2] = 0;
     240        CppAD::sparse_jacobian_work work;
     241        size_t n_sweep = f.SparseJacobianForward(x, s, r, c, jac, work);
     242        for(k = 0; k < 3; k++)
     243        {       ell = r[k] * n + c[k];
     244                ok &=  NearEqual(check[ell], jac[k], eps, eps);
     245        }
     246        ok &= (n_sweep == 1);
     247
     248        // Use reverse mode to compute rows 0 and 1
     249        // (make sure order of rows and columns does not matter)
     250        r.resize(4), c.resize(4); jac.resize(4);
     251        r[0] = 0; c[0] = 0;
     252        r[1] = 1; c[1] = 2;
     253        r[2] = 0; c[2] = 1;
     254        r[3] = 1; c[3] = 3;
     255        work.clear();
     256        n_sweep = f.SparseJacobianReverse(x, s, r, c, jac, work);
     257        for(k = 0; k < 4; k++)
     258        {       ell = r[k] * n + c[k];
     259                ok &=  NearEqual(check[ell], jac[k], eps, eps);
     260        }
     261        ok &= (n_sweep == 1);
     262
     263        return ok;
     264}
     265
     266
     267template <class VectorBase, class VectorBool>
     268bool reverse_bool(void)
     269{       bool ok = true;
     270        using CppAD::AD;
     271        using CppAD::NearEqual;
     272        size_t i, j, k;
     273
     274        // domain space vector
     275        size_t n = 4;
     276        CPPAD_TEST_VECTOR< AD<double> >  X(n);
     277        for(j = 0; j < n; j++)
     278                X[j] = AD<double> (0);
     279
     280        // declare independent variables and starting recording
     281        CppAD::Independent(X);
     282
     283        size_t m = 3;
     284        CPPAD_TEST_VECTOR< AD<double> >  Y(m);
     285        Y[0] = 1.0*X[0] + 2.0*X[1];
     286        Y[1] = 3.0*X[2] + 4.0*X[3];
     287        Y[2] = 5.0*X[0] + 6.0*X[1] + 7.0*X[2] + 8.0*X[3]*X[3]/2.;
     288
     289        // create f: X -> Y and stop tape recording
     290        CppAD::ADFun<double> f(X, Y);
     291
     292        // new value for the independent variable vector
     293        VectorBase x(n);
     294        for(j = 0; j < n; j++)
     295                x[j] = double(j);
     296
     297        // Jacobian of y without sparsity pattern
     298        VectorBase jac(m * n);
     299        jac = f.SparseJacobian(x);
     300        /*
     301              [ 1 2 0 0    ]
     302        jac = [ 0 0 3 4    ]
     303              [ 5 6 7 8*x_3]
     304        */
     305        VectorBase check(m * n);
     306        check[0] = 1.; check[1] = 2.; check[2]  = 0.; check[3]  = 0.;
     307        check[4] = 0.; check[5] = 0.; check[6]  = 3.; check[7]  = 4.;
     308        check[8] = 5.; check[9] = 6.; check[10] = 7.; check[11] = 8.*x[3];
     309        for(k = 0; k < 12; k++)
     310                ok &=  NearEqual(check[k], jac[k], 1e-10, 1e-10 );
     311
     312        // test passing sparsity pattern
     313        VectorBool s(m * m);
     314        VectorBool p(m * n);
     315        for(i = 0; i < m; i++)
     316        {       for(k = 0; k < m; k++)
     317                        s[i * m + k] = false;
     318                s[i * m + i] = true;
     319        }
     320        p   = f.RevSparseJac(m, s);
     321        jac = f.SparseJacobian(x);
     322        for(k = 0; k < 12; k++)
     323                ok &=  NearEqual(check[k], jac[k], 1e-10, 1e-10 );
     324
     325        return ok;
     326}
     327
     328template <class VectorBase, class VectorSet>
     329bool reverse_set(void)
    22330{       bool ok = true;
    23331        using CppAD::AD;
     
    64372
    65373        // test passing sparsity pattern
    66         VectorBool s(m * m);
    67         VectorBool p(m * n);
    68         for(i = 0; i < m; i++)
    69         {       for(k = 0; k < m; k++)
    70                         s[i * m + k] = false;
    71                 s[i * m + i] = true;
    72         }
    73         p   = f.RevSparseJac(m, s);
    74         jac = f.SparseJacobian(x);
    75         for(k = 0; k < 12; k++)
    76                 ok &=  NearEqual(check[k], jac[k], 1e-10, 1e-10 );
    77 
    78         return ok;
    79 }
    80 
    81 template <class VectorBase, class VectorSet>
    82 bool reverse_set()
    83 {       bool ok = true;
    84         using CppAD::AD;
    85         using CppAD::NearEqual;
    86         size_t i, j, k;
    87 
    88         // domain space vector
    89         size_t n = 4;
    90         CPPAD_TEST_VECTOR< AD<double> >  X(n);
    91         for(j = 0; j < n; j++)
    92                 X[j] = AD<double> (0);
    93 
    94         // declare independent variables and starting recording
    95         CppAD::Independent(X);
    96 
    97         size_t m = 3;
    98         CPPAD_TEST_VECTOR< AD<double> >  Y(m);
    99         Y[0] = X[0] + X[1];
    100         Y[1] = X[2] + X[3];
    101         Y[2] = X[0] + X[1] + X[2] + X[3] * X[3] / 2.;
    102 
    103         // create f: X -> Y and stop tape recording
    104         CppAD::ADFun<double> f(X, Y);
    105 
    106         // new value for the independent variable vector
    107         VectorBase x(n);
    108         for(j = 0; j < n; j++)
    109                 x[j] = double(j);
    110 
    111         // Jacobian of y without sparsity pattern
    112         VectorBase jac(m * n);
    113         jac = f.SparseJacobian(x);
    114         /*
    115               [ 1 1 0 0  ]
    116         jac = [ 0 0 1 1  ]
    117               [ 1 1 1 x_3]
    118         */
    119         VectorBase check(m * n);
    120         check[0] = 1.; check[1] = 1.; check[2]  = 0.; check[3]  = 0.;
    121         check[4] = 0.; check[5] = 0.; check[6]  = 1.; check[7]  = 1.;
    122         check[8] = 1.; check[9] = 1.; check[10] = 1.; check[11] = x[3];
    123         for(k = 0; k < 12; k++)
    124                 ok &=  NearEqual(check[k], jac[k], 1e-10, 1e-10 );
    125 
    126         // test passing sparsity pattern
    127374        VectorSet s(m), p(m);
    128375        for(i = 0; i < m; i++)
     
    137384
    138385template <class VectorBase, class VectorBool>
    139 bool forward_bool()
     386bool forward_bool(void)
    140387{       bool ok = true;
    141388        using CppAD::AD;
     
    201448
    202449template <class VectorBase, class VectorSet>
    203 bool forward_set()
     450bool forward_set(void)
    204451{       bool ok = true;
    205452        using CppAD::AD;
     
    260507}
    261508
    262 bool multiple_of_n_bit()
     509bool multiple_of_n_bit(void)
    263510{       bool ok = true;
    264511        using CppAD::AD;
     
    308555bool sparse_jacobian(void)
    309556{       bool ok = true;
     557        ok &= rc_tridiagonal();
     558        ok &= multiple_of_n_bit();
    310559        // ---------------------------------------------------------------
    311560        // vector of bool cases
    312         ok &= forward_bool< CppAD::vector<double>, CppAD::vectorBool   >();
    313         ok &= reverse_bool< CppAD::vector<double>, CppAD::vector<bool> >();
     561        ok &=      rc_bool< CppAD::vector<double>, CppAD::vectorBool   >();
     562        ok &= forward_bool< CppAD::vector<double>, CppAD::vector<bool> >();
    314563        //
    315         ok &= forward_bool< std::vector<double>,   std::vector<bool>   >();
    316         ok &= reverse_bool< std::vector<double>,   std::valarray<bool> >();
     564        ok &= reverse_bool< std::vector<double>,   std::vector<bool>   >();
     565        ok &=      rc_bool< std::vector<double>,   std::valarray<bool> >();
    317566        //
    318567        ok &= forward_bool< std::valarray<double>, CppAD::vectorBool   >();
     
    324573        typedef CppAD::vector< std::set<size_t> > cppad_vector_set;
    325574        //
    326         ok &= forward_set< CppAD::vector<double>, std_vector_set   >();
    327         ok &= reverse_set< std::valarray<double>, std_vector_set   >();
     575        ok &=      rc_set< CppAD::vector<double>, std_vector_set   >();
     576        ok &= forward_set< std::valarray<double>, std_vector_set   >();
    328577        //
    329         ok &= forward_set< std::vector<double>,   cppad_vector_set >();
    330         ok &= reverse_set< CppAD::vector<double>, cppad_vector_set >();
     578        ok &= reverse_set< std::vector<double>,   cppad_vector_set >();
     579        ok &=      rc_set< CppAD::vector<double>, cppad_vector_set >();
    331580        //
    332581        // According to section 26.3.2.3 of the 1998 C++ standard
     
    336585        // ---------------------------------------------------------------
    337586        //
    338         ok &= multiple_of_n_bit();
    339         return ok;
    340 }
     587        return ok;
     588}
Note: See TracChangeset for help on using the changeset viewer.