Changeset 1545


Ignore:
Timestamp:
Oct 13, 2009 8:09:15 AM (11 years ago)
Author:
bradbell
Message:

branches/sparsity: Change so user can use packed bool or vectors of sets.

*/for_sparse_jac.cpp: test use of vector of sets for user sparsity pattern.
build.sh: fix problem with test version of doxygen command.
glossary.omh: extend sparsity pattern to include representations.
whats_new_09.omh: user's view of the changes.
ad_fun.hpp: private member helper functions for dispatching calculations.
for_sparse_jac.hpp: include new vector of sets option.

Location:
branches/sparsity
Files:
7 edited

Legend:

Unmodified
Added
Removed
  • branches/sparsity/build.sh

    r1517 r1545  
    448448        if ! rm doxygen.$$
    449449        then
    450                 echo "Error: cannot remove doxygen.err"
     450                echo "Error: cannot remove doxygen.$$"
    451451                exit 1
    452452        fi
    453453        if ! ./check_doxygen.sh
    454454        then
    455                 echo "Warnings of doxygen output."
     455                dir=`pwd`
     456                echo "Error: check_doxygen.sh failed; see $dir/doxygen.log."
    456457                exit 1
    457458        fi
     
    593594        done
    594595        # Test developer documentation ---------------------------------------
    595         echo "doxygen doxyfile >& doxygen.log"
    596         if ! doxygen doxyfile >& doxygen.log
    597         then
    598                 echo "Error: doxygen doxyfile"
     596        # avoid warning in mid sentence of other log output by separating them
     597        echo "doxygen doxyfile 1> doxygen.log 2> doxygen.$$"
     598        if ! doxygen doxyfile 1> doxygen.log 2> doxygen.$$
     599        then
     600                echo "Error: during doxygen"
     601                exit 1
     602        fi
     603        echo "cat doxygen.$$ >> doxygen.log"
     604        if ! cat doxygen.$$ >> doxygen.log
     605        then
     606                echo "Error: cannot add errors and warnings to doxygen.log"
     607                exit 1
     608        fi
     609        echo "rm doxygen.$$"
     610        if ! rm doxygen.$$
     611        then
     612                echo "Error: cannot remove doxygen.$$"
    599613                exit 1
    600614        fi
    601615        if ! ./check_doxygen.sh
    602616        then
    603                 echo "Warnings of doxygen output; see doxygen.log."
     617                dir=`pwd`
     618                echo "Error: check_doxygen.sh failed, see $dir/doxygen.log."
    604619                exit 1
    605620        fi
  • branches/sparsity/cppad/local/ad_fun.hpp

    r1539 r1545  
    7777public:
    7878// ------------------------------------------------------------
     79// Private member variables
    7980private:
    80 
    8181        /// debug checking number of comparision operations that changed
    8282        size_t compare_change_;
     
    116116        sparse_set       for_jac_sparse_set_;
    117117
     118// ------------------------------------------------------------
     119// Private member functions
     120
    118121        /// change the operation sequence corresponding to this object
    119122        template <typename ADvector>
    120123        void Dependent(ADTape<Base> *tape, const ADvector &y);
     124
     125        /// vector of bool version of ForSparseJac
     126        template <class VectorSet>
     127        void ForSparseJacCase(
     128                bool               set_type  ,
     129                size_t             q         ,
     130                const VectorSet&   r         , 
     131                VectorSet&         s
     132        );
     133        /// vector of std::set<size_t> version of ForSparseJac
     134        template <class VectorSet>
     135        void ForSparseJacCase(
     136                const std::set<size_t>&  set_type  ,
     137                size_t                   q         ,
     138                const VectorSet&         r         , 
     139                VectorSet&               s
     140        );
    121141
    122142// ------------------------------------------------------------
     
    155175        // forward mode Jacobian sparsity
    156176        // (see doxygen documentation in for_sparse_jac.hpp)
    157         template <typename VectorBool>
    158         VectorBool ForSparseJac(
    159                 size_t q, const VectorBool &Px, bool packed = true
     177        template <typename VectorSet>
     178        VectorSet ForSparseJac(
     179                size_t q, const VectorSet &r
    160180        );
    161181        // reverse mode Jacobian sparsity
  • branches/sparsity/cppad/local/for_sparse_jac.hpp

    r1539 r1545  
    1717$begin ForSparseJac$$
    1818$spell
     19        std
    1920        var
    2021        Jacobian
     
    3435
    3536$head Syntax$$
    36 $syntax%%s% = %f%.ForSparseJac(%q%, %r%)%$$
    37 $pre
    38 $$
    39 $syntax%%s% = %f%.ForSparseJac(%q%, %r%, %packed%)%$$
     37$icode%s% = %f%.ForSparseJac(%q%, %r%)%$$
    4038
    4139$head Purpose$$
    4240We use $latex F : B^n \rightarrow B^m$$ to denote the
    43 $xref/glossary/AD Function/AD function/$$ corresponding to $italic f$$.
     41$xref/glossary/AD Function/AD function/$$ corresponding to $icode f$$.
    4442For a fixed $latex n \times q$$ matrix $latex R$$,
    4543the Jacobian of $latex F[ x + R * u ]$$
     
    5452
    5553$head f$$
    56 The object $italic f$$ has prototype
    57 $syntax%
     54The object $icode f$$ has prototype
     55$codei%
    5856        ADFun<%Base%> %f%
    5957%$$
    60 Note that the $xref/ADFun/$$ object $italic f$$ is not $code const$$.
     58Note that the $xref/ADFun/$$ object $icode f$$ is not $code const$$.
    6159After this the sparsity pattern
    6260for each of the variables in the operation sequence
    63 is stored in the object $italic f$$.
     61is stored in the object $icode f$$.
    6462
    6563
     
    7068
    7169$head q$$
    72 The argument $italic q$$ has prototype
    73 $syntax%
     70The argument $icode q$$ has prototype
     71$codei%
    7472        size_t %q%
    7573%$$
    76 It specifies the number of columns in the Jacobian $latex J(x)$$.
    77 Note that the memory required for the calculation is proportional
    78 to $latex q$$ times the total number of variables
    79 in the AD operation sequence corresponding to $italic f$$
    80 ($xref/SeqProperty/size_var/f.size_var/$$).
    81 Smaller values for $italic q$$ can be used to
    82 break the sparsity calculation
    83 into groups that do not require to much memory.
     74It specifies the number of columns in
     75$latex R$$ and the Jacobian $latex J(x)$$.
    8476
    8577$head r$$
    86 The argument $italic r$$ has prototype
    87 $syntax%
    88         const %VectorBool% &%r%
     78The argument $icode r$$ has prototype
     79$codei%
     80        const %VectorSet% &%r%
    8981%$$
    90 (see $xref/ForSparseJac/VectorBool/VectorBool/$$ below)
    91 and its size is $latex n * q$$.
     82(see $xref/ForSparseJac/VectorSet/VectorSet/$$ below).
     83If it has elements of type $code bool$$,
     84its size is $latex n * q$$.
     85If it has elements of type $code std::set<size_t>$$,
     86its size is $latex n$$ and all the set elements must be between
     87zero and $icode%q%-1%$$.
    9288It specifies a
    9389$xref/glossary/Sparsity Pattern/sparsity pattern/$$
    94 for the matrix $italic R$$ as follows:
    95 for $latex i = 0 , \ldots , n-1$$ and $latex j = 0 , \ldots , q-1$$.
    96 $latex \[
    97         R_{i,j} \neq 0 ; \Rightarrow \; r [ i * q + j ] = {\rm true}
    98 \] $$
    99 
    100 $head packed$$
    101 If $italic packed$$ is true,
    102 during the sparsity calculation sets of indices are represented
    103 as vectors of bits that packed into words and operations are done
    104 on multiple bits at a time (the number of bits in a word is unspecified).
    105 Otherwise, sets of indices are represents using a sparse structure
    106 that only includes the non-zero indices and operations are done
    107 one index at a time.
    108 $pre
    109 
    110 $$
    111 The default value for $italic packed$$ is true; i.e.,
    112 the value used if it is not present.
     90for the matrix $icode R$$.
    11391
    11492$head s$$
    115 The return value $italic s$$ has prototype
    116 $syntax%
    117         %VectorBool% %s%
     93The return value $icode s$$ has prototype
     94$codei%
     95        %VectorSet% %s%
    11896%$$
    119 (see $xref/ForSparseJac/VectorBool/VectorBool/$$ below)
    120 and its size is $latex m * q$$.
     97(see $xref/ForSparseJac/VectorSet/VectorSet/$$ below).
     98If it has elements of type $code bool$$
     99its size is $latex m * q$$.
     100If it has elements of type $code std::set<size_t>$$,
     101its size if $latex m$$.
    121102It specifies a
    122103$xref/glossary/Sparsity Pattern/sparsity pattern/$$
    123 for the matrix $latex J(x)$$ as follows:
    124 for $latex x \in B^n$$,
    125 for $latex i = 0 , \ldots , m-1$$,
    126 and $latex j = 0 , \ldots , q-1$$
     104for the matrix $latex J(x)$$ where
    127105$latex \[
    128         J(x)_{i,j} \neq 0 ; \Rightarrow \; s [ i * q + j ] = {\rm true}
     106        J(x)_{i,j} = \D{ F_i }{ x_j } (x)
    129107\] $$
    130 
    131 $head VectorBool$$
    132 The type $italic VectorBool$$ must be a $xref/SimpleVector/$$ class with
    133 $xref/SimpleVector/Elements of Specified Type/elements of type bool/$$.
    134 The routine $xref/CheckSimpleVector/$$ will generate an error message
    135 if this is not the case.
    136 In order to save memory,
    137 you may want to use a class that packs multiple elements into one
    138 storage location; for example,
    139 $xref/CppAD_vector/vectorBool/vectorBool/$$.
     108where $latex x \in B^n$$ is any value
     109and $latex F(x)$$ corresponds to the operation sequence
     110currently stored in $icode f$$
     111
     112$head VectorSet$$
     113The type $icode VectorSet$$ must be a $xref/SimpleVector/$$ class with
     114$xref/SimpleVector/Elements of Specified Type/elements of type/$$
     115$code bool$$ or $code std::set<size_t>$$.
    140116
    141117$head Entire Sparsity Pattern$$
    142118Suppose that $latex q = n$$ and
    143 $latex R$$ is the $latex n \times n$$ identity matrix,
    144 If follows that
    145 $latex \[
    146 r [ i * q + j ] = \left\{ \begin{array}{ll}
    147         {\rm true}  & {\rm if} \; i = j \\
    148         {\rm flase} & {\rm otherwise}
    149 \end{array} \right.
    150 \] $$
    151 is an efficient sparsity pattern for $latex R$$;
    152 i.e., the choice for $italic r$$ has as few true values as possible.
     119$latex R$$ is the $latex n \times n$$ identity matrix.
    153120In this case,
    154 the corresponding value for $italic s$$ is a
     121the corresponding value for $icode s$$ is a
    155122sparsity pattern for the Jacobian $latex J(x) = F^{(1)} ( x )$$.
    156123
     
    179146The C++ source code corresponding to this operation is
    180147\verbatim
    181         s = f.ForSparseJac(q, r, packed)
     148        s = f.ForSparseJac(q, r)
    182149\endverbatim
    183150
     
    185152is the base type for this recording.
    186153
    187 \tparam VectorBool
    188 is a simple vector with elements of type bool.
     154\tparam VectorSet
     155is a simple vector class with elements of type \c bool.
    189156
    190157\tparam Sparsity
     
    198165
    199166\param s
    200 the input value of \a s must be a vector with size \c m*q
     167The input value of \a s must be a vector with size \c m*q
    201168where \c m is the number of dependent variables
    202169corresponding to the operation sequence stored in \a play.
    203 The input value of its elements does not matter.
     170The input value of the components of \c s does not matter.
    204171On output, \a s is the sparsity pattern for the matrix
    205172\f[
     
    231198tape (given the sparsity pattern for the independent variables is \f$ R \f$).
    232199*/
    233 template <class Base, class VectorBool, class Sparsity>
    234 void ForSparseJac(
     200template <class Base, class VectorSet, class Sparsity>
     201void ForSparseJacBool(
    235202        size_t                 q                ,
    236         const VectorBool&      r                ,
    237         VectorBool&            s                ,
     203        const VectorSet&       r                ,
     204        VectorSet&             s                ,
    238205        size_t                 total_num_var    ,
    239206        CppAD::vector<size_t>& dep_taddr        ,
     
    245212        size_t i, j;
    246213
    247         // check VectorBool is Simple Vector class with bool elements
    248         CheckSimpleVector<bool, VectorBool>();
     214        // check VectorSet is Simple Vector class with bool elements
     215        CheckSimpleVector<bool, VectorSet>();
    249216
    250217        // range and domain dimensions for F
     
    303270}
    304271
     272/*!
     273Calculate Jacobian sparsity patterns using forward mode.
     274
     275The C++ source code corresponding to this operation is
     276\verbatim
     277        s = f.ForSparseJac(q, r)
     278\endverbatim
     279
     280\tparam Base
     281is the base type for this recording.
     282
     283\tparam VectorSet
     284is a simple vector class with elements of the type \c std::set<size_t>.
     285
     286\tparam Sparsity
     287is either \c sparse_pack or \c sparse_set.
     288
     289\param q
     290is the number of columns in the matrix \f$ R \f$.
     291
     292\param r
     293is a sparsity pattern for the matrix \f$ R \f$.
     294
     295\param s
     296The input value of \a s must be a vector with size \c m
     297where \c m is the number of dependent variables
     298corresponding to the operation sequence stored in \a play.
     299On input, each element of \c s must be an empty set.
     300On output, \a s is the sparsity pattern for the matrix
     301\f[
     302        J(x) = F^{(1)} (x) * R
     303\f]
     304where \f$ F \f$ is the function corresponding to the operation sequence
     305and \a x is any argument value.
     306
     307\param total_num_var
     308is the total number of variable in this recording.
     309
     310\param dep_taddr
     311maps dependendent variable index
     312to the corresponding variable in the tape.
     313
     314\param ind_taddr
     315maps independent variable index
     316to the corresponding variable in the tape.
     317
     318\param play
     319is the recording that defines the function we are computing the sparsity
     320pattern for.
     321
     322\param for_jac_sparsity
     323the input value of \a for_jac_sparsity does not matter.
     324On output, \a for_jac_sparsity.n_set() == \a total_num_var
     325and \a for_jac_sparsity.end() == \a q.
     326It contains the forward sparsity pattern for all of the variables on the
     327tape (given the sparsity pattern for the independent variables is \f$ R \f$).
     328*/
     329
     330template <class Base, class VectorSet, class Sparsity>
     331void ForSparseJacSet(
     332        size_t                  q                ,
     333        const VectorSet&        r                ,
     334        VectorSet&              s                ,
     335        size_t                  total_num_var    ,
     336        CppAD::vector<size_t>&  dep_taddr        ,
     337        CppAD::vector<size_t>&  ind_taddr        ,
     338        CppAD::player<Base>&    play             ,
     339        Sparsity&               for_jac_sparsity )
     340{
     341        // temporary indices
     342        size_t i, j;
     343        std::set<size_t>::iterator itr;
     344
     345        // cannot use CheckSimpleVector when vector elements are sets
     346        // because cannot assign an inteter to a set.
     347        // CheckSimpleVector<std::set<size_t>, VectorSet>();
     348
     349        // range and domain dimensions for F
     350        size_t m = dep_taddr.size();
     351        size_t n = ind_taddr.size();
     352
     353        CPPAD_ASSERT_KNOWN(
     354                q > 0,
     355                "ForSparseJac: q (first arugment) is not greater than zero"
     356        );
     357
     358        CPPAD_ASSERT_KNOWN(
     359                r.size() == n,
     360                "ForSparseJac: r (second argument) length is not equal to\n"
     361                "q (first argument) times domain dimension for ADFun object."
     362        );
     363
     364        // allocate memory for the requested sparsity calculation
     365        for_jac_sparsity.resize(total_num_var, q);
     366
     367        // set values corresponding to independent variables
     368        for(i = 0; i < n; i++)
     369        {       CPPAD_ASSERT_UNKNOWN( ind_taddr[i] < total_num_var );
     370                // ind_taddr[i] is operator taddr for i-th independent variable
     371                CPPAD_ASSERT_UNKNOWN( play.GetOp( ind_taddr[i] ) == InvOp );
     372
     373                // add the elements that are present
     374                for( itr = r[i].begin(); itr != r[i].end(); itr++)
     375                {       j = *itr;
     376                        CPPAD_ASSERT_KNOWN(
     377                                j < q,
     378                                "ForSparseJac: a set element have value "
     379                                " greater than or equal q."
     380                        );
     381                        for_jac_sparsity.add_element( ind_taddr[i], j);
     382                }
     383        }
     384
     385        // evaluate the sparsity patterns
     386        ForJacSweep(
     387                n,
     388                total_num_var,
     389                &play,
     390                for_jac_sparsity
     391        );
     392
     393        // return values corresponding to dependent variables
     394        CPPAD_ASSERT_UNKNOWN( s.size() == m );
     395        for(i = 0; i < m; i++)
     396        {       CPPAD_ASSERT_UNKNOWN( dep_taddr[i] < total_num_var );
     397                CPPAD_ASSERT_UNKNOWN( s[i].empty() );
     398
     399                // add elements bits
     400                CPPAD_ASSERT_UNKNOWN( for_jac_sparsity.end() == q );
     401                for_jac_sparsity.begin( dep_taddr[i] );
     402                j = for_jac_sparsity.next_element();
     403                while( j < q )
     404                {       s[i].insert(j);
     405                        j = for_jac_sparsity.next_element();
     406                }
     407        }
     408}
     409
    305410
    306411/*!
     
    309414The C++ source code corresponding to this operation is
    310415\verbatim
    311         s = f.ForSparseJac(q, r, packed)
     416        s = f.ForSparseJac(q, r)
    312417\endverbatim
    313418
     
    315420is the base type for this recording.
    316421
    317 \tparam VectorBool
    318 is a simple vector with elements of type bool.
     422\tparam VectorSet
     423is a simple vector with elements of type \c bool
     424or \c std::set<size_t>.
    319425
    320426\param q
     
    324430is a sparsity pattern for the matrix \f$ R \f$.
    325431
    326 \param packed
    327 If \a packed is true, the type \c sparse_pack is used for the calculations.
    328 Otherwise the type \c sparse_set is used for the calculations.
    329 
    330432\return
     433If \c VectorSet::value_type is \c bool,
    331434the return value \c s is a vector with size \c m*q
    332435where \c m is the number of dependent variables
    333436corresponding to the operation sequence stored in \c f.
     437If \c VectorSet::value_type is \c std::set<size_t>,
     438the return value \c s is a vector of sets with size \c m
     439and with all its elements between zero and \a q - 1.
    334440The value of \a s is the sparsity pattern for the matrix
    335441\f[
     
    364470*/
    365471template <class Base>
    366 template <class VectorBool>
    367 VectorBool ADFun<Base>::ForSparseJac(
     472template <class VectorSet>
     473VectorSet ADFun<Base>::ForSparseJac(
    368474        size_t             q             ,
    369         const VectorBool&  r             ,
    370         bool               packed        )
    371 {       size_t m = dep_taddr_.size();
    372         VectorBool s( m * q );
    373 
    374         if( packed )
    375         {       // free any results stored in for_jac_sparse_set_       
    376                 for_jac_sparse_set_.resize(0, 0);
    377                 // store results in for_jac_sparse_pack_
    378                 CppAD::ForSparseJac(
    379                         q                ,
    380                         r                ,
    381                         s                ,
    382                         total_num_var_   ,
    383                         dep_taddr_       ,
    384                         ind_taddr_       ,
    385                         play_            ,
    386                         for_jac_sparse_pack_
    387                 );
    388         }
    389         else
    390         {       // free any results stored in for_jac_sparse_pack_
    391                 for_jac_sparse_pack_.resize(0, 0);
    392                 // store results in for_jac_sparse_set_
    393                 CppAD::ForSparseJac(
    394                         q                ,
    395                         r                ,
    396                         s                ,
    397                         total_num_var_   ,
    398                         dep_taddr_       ,
    399                         ind_taddr_       ,
    400                         play_            ,
    401                         for_jac_sparse_set_
    402                 );
    403         }
     475        const VectorSet&   r             )
     476{       VectorSet s;
     477        typedef typename VectorSet::value_type Set_type;
     478
     479        ForSparseJacCase(
     480                Set_type()  ,
     481                q           ,
     482                r           ,
     483                s
     484        );
     485
    404486        return s;
    405487}
    406488
     489/*!
     490Helper function for ForSparseJac(q, r).
     491
     492All of the description in the public member function ForSparseJac(q, r)
     493applies.
     494
     495\param set_type
     496is a bool value. This argument is used to dispatch to the proper source
     497code depending on the value of \c VectroSet::value_type.
     498
     499\param q
     500See ForSparseJac(q, r).
     501
     502\param r
     503See ForSparseJac(q, r).
     504
     505\param s
     506is the return value for the corresponding call to ForSparseJac(q, r).
     507*/
     508
     509template <class Base>
     510template <class VectorSet>
     511void ADFun<Base>::ForSparseJacCase(
     512        bool                set_type      ,
     513        size_t              q             ,
     514        const VectorSet&    r             ,
     515        VectorSet&          s             )
     516{       size_t m = dep_taddr_.size();
     517
     518        // dimension size of result vector
     519        s.resize( m * q );
     520
     521        // free any results stored in for_jac_sparse_set_       
     522        for_jac_sparse_set_.resize(0, 0);
     523
     524        // store results in for_jac_sparse_pack_
     525        CppAD::ForSparseJacBool(
     526                q                ,
     527                r                ,
     528                s                ,
     529                total_num_var_   ,
     530                dep_taddr_       ,
     531                ind_taddr_       ,
     532                play_            ,
     533                for_jac_sparse_pack_
     534        );
     535}
     536
     537
     538/*!
     539Helper function for ForSparseJac(q, r).
     540
     541All of the description in the public member function ForSparseJac(q, r)
     542applies.
     543
     544\param set_type
     545is a std::set<size_t> object.
     546This argument is used to dispatch to the proper source
     547code depending on the value of \c VectroSet::value_type.
     548
     549\param q
     550See ForSparseJac(q, r).
     551
     552\param r
     553See ForSparseJac(q, r).
     554
     555\param s
     556is the return value for the corresponding call to ForSparseJac(q, r).
     557*/
     558template <class Base>
     559template <class VectorSet>
     560void ADFun<Base>::ForSparseJacCase(
     561        const std::set<size_t>&    set_type      ,
     562        size_t                     q             ,
     563        const VectorSet&           r             ,
     564        VectorSet&                 s             )
     565{       size_t m = dep_taddr_.size();
     566
     567        // dimension size of result vector
     568        s.resize( m );
     569
     570        // free any results stored in for_jac_sparse_set_       
     571        for_jac_sparse_pack_.resize(0, 0);
     572
     573        // store results in for_jac_sparse_pack_
     574        CppAD::ForSparseJacSet(
     575                q                ,
     576                r                ,
     577                s                ,
     578                total_num_var_   ,
     579                dep_taddr_       ,
     580                ind_taddr_       ,
     581                play_            ,
     582                for_jac_sparse_set_
     583        );
     584}
     585
     586
    407587CPPAD_END_NAMESPACE
    408588# endif
  • branches/sparsity/example/for_sparse_jac.cpp

    r1370 r1545  
    11/* $Id$ */
    22/* --------------------------------------------------------------------------
    3 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-07 Bradley M. Bell
     3CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-09 Bradley M. Bell
    44
    55CppAD is distributed under multiple licenses. This distribution is under
     
    3535// BEGIN PROGRAM
    3636
     37# include <set>
    3738# include <cppad/cppad.hpp>
     39
    3840namespace { // -------------------------------------------------------------
    39 // define the template function ForSparseJacCases<Vector> in empty namespace
     41// define the template function BoolCases<Vector>
    4042template <typename Vector>
    41 bool ForSparseJacCases(void)
     43bool BoolCases(void)
    4244{       bool ok = true;
    4345        using CppAD::AD;
     
    8587        return ok;
    8688}
     89// define the template function SetCases<Vector>
     90template <typename Vector>
     91bool SetCases(void)
     92{       bool ok = true;
     93        using CppAD::AD;
     94
     95        // domain space vector
     96        size_t n = 2;
     97        CPPAD_TEST_VECTOR< AD<double> > X(n);
     98        X[0] = 0.;
     99        X[1] = 1.;
     100
     101        // declare independent variables and start recording
     102        CppAD::Independent(X);
     103
     104        // range space vector
     105        size_t m = 3;
     106        CPPAD_TEST_VECTOR< AD<double> > Y(m);
     107        Y[0] = X[0];
     108        Y[1] = X[0] * X[1];
     109        Y[2] = X[1];
     110
     111        // create f: X -> Y and stop tape recording
     112        CppAD::ADFun<double> f(X, Y);
     113
     114        // sparsity pattern for the identity matrix
     115        Vector r(n);
     116        size_t i;
     117        for(i = 0; i < n; i++)
     118        {       assert( r[i].empty() );
     119                r[i].insert(i);
     120        }
     121
     122        // sparsity pattern for F'(x)
     123        Vector s(m);
     124        s = f.ForSparseJac(n, r);
     125
     126        // an interator to a standard set
     127        std::set<size_t>::iterator itr;
     128        bool found;
     129
     130        // Y[0] does     depend on X[0]
     131        found = s[0].find(0) != s[0].end();  ok &= ( found == true ); 
     132
     133        // Y[0] does not depend on X[1]
     134        found = s[0].find(1) != s[0].end();  ok &= ( found == false );
     135
     136        // Y[1] does     depend on X[0]
     137        found = s[1].find(0) != s[1].end();  ok &= ( found == true ); 
     138
     139        // Y[1] does     depend on X[1]
     140        found = s[1].find(1) != s[1].end();  ok &= ( found == true ); 
     141
     142        // Y[2] does not depend on X[0]
     143        found = s[2].find(0) != s[2].end();  ok &= ( found == false ); 
     144
     145        // Y[2] does     depend on X[1]
     146        found = s[2].find(1) != s[2].end();  ok &= ( found == true ); 
     147
     148        return ok;
     149}
    87150} // End empty namespace
    88151# include <vector>
     
    92155        // Run with Vector equal to four different cases
    93156        // all of which are Simple Vectors with elements of type bool.
    94         ok &= ForSparseJacCases< CppAD::vectorBool     >();
    95         ok &= ForSparseJacCases< CppAD::vector  <bool> >();
    96         ok &= ForSparseJacCases< std::vector    <bool> >();
    97         ok &= ForSparseJacCases< std::valarray  <bool> >();
     157        ok &= BoolCases< CppAD::vectorBool     >();
     158        ok &= BoolCases< CppAD::vector  <bool> >();
     159        ok &= BoolCases< std::vector    <bool> >();
     160        ok &= BoolCases< std::valarray  <bool> >();
     161
     162        // Run with Vector equal to three different cases all of which are
     163        // Simple Vectors with elements of type std::set<size_t>
     164        typedef std::set<size_t> set;
     165        ok &= SetCases< CppAD::vector  <set> >();
     166        ok &= SetCases< std::vector    <set> >();
     167        ok &= SetCases< std::valarray  <set> >();
    98168
    99169        return ok;
  • branches/sparsity/omh/glossary.omh

    r1370 r1545  
    11$Id$
    22/* --------------------------------------------------------------------------
    3 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-08 Bradley M. Bell
     3CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-09 Bradley M. Bell
    44
    55CppAD is distributed under multiple licenses. This distribution is under
     
    1414$aindex head subhead$$
    1515$spell
     16        std
    1617        Vec
    1718        cos
     
    142143$index pattern, sparsity$$
    143144$index efficient, sparsity$$
     145CppAD describes a sparse matrix as a vector of sets with each
     146vector component corresponding to a row
     147and the elements of the set corresponding to the possibly non-zero columns.
     148A vector of $code bool$$
     149can represent a vector of sets using one bit per element.
     150(Some vectors of $code bool$$ use one byte per element but
     151$cref/vectorBool/CppAD_vector/vectorBool/$$ is an example class that
     152uses one bit per element.)
     153The problem is that this representation uses one bit for both the elements
     154that are there and the ones that are not.
     155$pre
     156
     157$$
     158A vector of $code std::set<size_t>$$ does not
     159represent the elements that are not present,
     160but it uses about three $code size_t$$ values
     161for each element that is present.
     162For example, if $code size_t$$ uses 32 bits,
     163a vector of $code std::set<size_t>$$ uses
     164about 100 bits for each element that is present in the
     165vector of sets.
     166Thus, a vector of $code std::set<size_t>$$ should be more efficient for
     167very sparse matrix representations.
     168
     169$subhead Vector of Boolean$$
    144170Given a matrix
    145171$latex A \in B^{m \times n}$$,
    146 a boolean valued $latex m \times n$$ matrix $latex P$$ is a
    147 sparsity pattern for $latex A$$ if
     172a vector of $code bool$$ $latex B$$ of length $latex m \times n$$
     173is a sparsity pattern for $latex A$$ if
    148174for $latex i = 0, \ldots , m-1$$ and $latex j = 0 , \ldots n-1$$,
    149175$latex \[
    150176A_{i,j} \neq 0 
    151177\; \Rightarrow \;
    152 P_{i,j} = {\rm true}
    153 \] $$
    154 Given two sparsity patterns $latex P$$ and $italic Q$$
    155 for a matrix $italic A$$, we say that $italic P$$ is more efficient than
    156 $italic Q$$ if $italic P$$ has fewer true elements than $italic Q$$.
     178B_{i * n + j} = {\rm true}
     179\] $$
     180Given two sparsity patterns $latex B$$ and $italic C$$
     181for a matrix $italic A$$, we say that $italic B$$ is more efficient than
     182$italic C$$ if $italic B$$ has fewer true elements than $italic C$$.
     183For example, if $latex A$$ is the identity matrix,
     184$latex \[
     185        B_{i * n + j} = (i = j)
     186\] $$
     187defines an efficient sparsity pattern for $latex A$$.
     188
     189$subhead Vector of Sets$$
     190Given a matrix
     191$latex A \in B^{m \times n}$$,
     192a vector of sets $latex S$$ of length $latex m$$ is a
     193sparsity pattern for $latex A$$ if
     194for $latex i = 0, \ldots , m-1$$
     195$latex \[
     196A_{i,j} \neq 0 
     197\; \Rightarrow \; j \in S_i
     198\] $$
     199Given two sparsity patterns $latex S$$ and $italic T$$
     200for a matrix $italic A$$, we say that $italic S$$ is more efficient than
     201$italic T$$ if $italic S$$ has fewer elements than $italic T$$.
     202For example, if $latex A$$ is the identity matrix,
     203$latex \[
     204        S_i =  \{ j \}
     205\] $$
     206defines an efficient sparsity pattern for $latex A$$.
    157207
    158208$head Tape$$
  • branches/sparsity/omh/whats_new_09.omh

    r1538 r1545  
    1313$begin whats_new_09$$
    1414$spell
     15        bool
    1516        Microsoft
    1617        retape
     
    5758trying to read and understand the CppAD source code.)
    5859
     60$head 10-08$$
     61The $icode packed$$ parameter for the sparsity routine
     62$cref/ForSparseJac/$$
     63(introduced on $cref/09-26/whats_new_09/09-26/$$) has been removed.
     64It has been replaced by changing the argument and return values
     65to be more versatile.
     66To be specific, they can now represent sparsity
     67using vectors of sets instead of just vectors of $code bool$$
     68(see $cref/sparsity patterns/glossary/Sparsity Pattern/$$).
     69
    5970$head 10-03$$
    6071The Microsoft Visual Studio project files for
     
    8697
    8798$head 09-26$$
    88 Added the $cref/packed/ForSparseJac/packed/$$ parameter to
     99Added the $code packed$$ parameter to
    89100$cref/ForSparseJac/$$ and $cref/RevSparseJac/$$.
    90101If $icode packed$$ is false,
  • branches/sparsity/test_more/for_sparse_jac.cpp

    r1532 r1545  
    7171namespace { // Begin empty namespace
    7272
    73 bool case_one(bool packed)
     73bool case_one()
    7474{       bool ok = true;
    7575        using namespace CppAD;
     
    149149        ADFun<double> F(X, Y);
    150150
     151        // ---------------------------------------------------------
    151152        // dependency matrix for the identity function W(x) = x
    152153        CPPAD_TEST_VECTOR< bool > Px(n * n);
     
    160161        // evaluate the dependency matrix for F(X(x))
    161162        CPPAD_TEST_VECTOR< bool > Py(m * n);
    162         Py = F.ForSparseJac(n, Px, packed);
     163        Py = F.ForSparseJac(n, Px);
    163164
    164165        // check values
     
    168169        }       
    169170
     171        // ---------------------------------------------------------
     172        // dependency matrix for the identity function W(x) = x
     173        CPPAD_TEST_VECTOR< std::set<size_t> > Sx(n);
     174        for(i = 0; i < n; i++)
     175        {       assert( Sx[i].empty() );
     176                Sx[i].insert(i);
     177        }
     178
     179        // evaluate the dependency matrix for F(X(x))
     180        CPPAD_TEST_VECTOR< std::set<size_t> > Sy(m);
     181        Sy = F.ForSparseJac(n, Sx);
     182
     183        // check values
     184        std::set<size_t>::iterator itr;
     185        bool found;
     186        for(i = 0; i < m; i++)
     187        {       for(j = 0; j < n; j++)
     188                {       found = Sy[i].find(j) != Sy[i].end();
     189                        ok &= (found == Check[i * n + j]);
     190                }
     191        }       
     192
    170193        return ok;
    171194}
    172195
    173 bool case_two(bool packed)
     196bool case_two()
    174197{       bool ok = true;
    175198        using namespace CppAD;
     
    235258        ADFun<double> F(X, Y);
    236259
     260        // -----------------------------------------------------------------
    237261        // dependency matrix for the identity function W(x) = x
    238262        CPPAD_TEST_VECTOR< bool > Px(n * n);
     
    246270        // evaluate the dependency matrix for F(X(x))
    247271        CPPAD_TEST_VECTOR< bool > Py(m * n);
    248         Py = F.ForSparseJac(n, Px, packed);
     272        Py = F.ForSparseJac(n, Px);
    249273
    250274        // check values
     
    254278        }       
    255279
     280        // ---------------------------------------------------------
     281        // dependency matrix for the identity function W(x) = x
     282        CPPAD_TEST_VECTOR< std::set<size_t> > Sx(n);
     283        for(i = 0; i < n; i++)
     284        {       assert( Sx[i].empty() );
     285                Sx[i].insert(i);
     286        }
     287
     288        // evaluate the dependency matrix for F(X(x))
     289        CPPAD_TEST_VECTOR< std::set<size_t> > Sy(m);
     290        Sy = F.ForSparseJac(n, Sx);
     291
     292        // check values
     293        std::set<size_t>::iterator itr;
     294        bool found;
     295        for(i = 0; i < m; i++)
     296        {       for(j = 0; j < n; j++)
     297                {       found = Sy[i].find(j) != Sy[i].end();
     298                        ok &= (found == Check[i * n + j]);
     299                }
     300        }       
     301
    256302        return ok;
    257303}
     
    262308{       bool ok = true;
    263309
    264         ok &= case_one(true);
    265         ok &= case_two(true);
    266 
    267         ok &= case_one(false);
    268         ok &= case_two(false);
     310        ok &= case_one();
     311        ok &= case_two();
    269312
    270313        return ok;
Note: See TracChangeset for help on using the changeset viewer.