Changeset 1533


Ignore:
Timestamp:
Sep 28, 2009 11:26:52 AM (11 years ago)
Author:
bradbell
Message:

trunk: Convert reverse sparse Hessian to work with vector_set.

rev_sparse_hes.cpp: test vector_set calculation of Hessian sparsity.
whats_new_09.omh: user's view of the changes.
for_sparse_jac.hpp: move member variables to end of calling sequence.
rev_sparse_jac.hpp: move member variables to end of calling sequence.
rev_sparse_hes.hpp: convert to use vector_set when for sparse does.

Location:
trunk
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • trunk/cppad/local/for_sparse_jac.hpp

    r1532 r1533  
    173173template <class Base, class VectorBool, class VectorSet>
    174174void ForSparseJac(
    175         size_t                 total_num_var    ,
    176175        size_t                 q                ,
    177176        const VectorBool&      r                ,
     177        VectorBool&            s                ,
     178        size_t                 total_num_var    ,
    178179        CppAD::vector<size_t>& dep_taddr        ,
    179180        CppAD::vector<size_t>& ind_taddr        ,
    180181        CppAD::player<Base>&   play             ,
    181         VectorBool&            s                ,
    182182        VectorSet&             for_jac_sparsity )
    183183{
     
    256256                // store results in for_jac_sparsity_
    257257                CppAD::ForSparseJac(
    258                         total_num_var_   ,
    259258                        q                ,
    260259                        r                ,
     260                        s                ,
     261                        total_num_var_   ,
    261262                        dep_taddr_       ,
    262263                        ind_taddr_       ,
    263264                        play_            ,
    264                         s                ,
    265265                        for_jac_sparsity_
    266266                );
     
    271271                // store results in for_jac_sparse_set_
    272272                CppAD::ForSparseJac(
    273                         total_num_var_   ,
    274273                        q                ,
    275274                        r                ,
     275                        s                ,
     276                        total_num_var_   ,
    276277                        dep_taddr_       ,
    277278                        ind_taddr_       ,
    278279                        play_            ,
    279                         s                ,
    280280                        for_jac_sparse_set_
    281281                );
  • trunk/cppad/local/rev_sparse_hes.hpp

    r1531 r1533  
    9191It must be the same value as in the previous $xref/ForSparseJac/$$ call
    9292$syntax%
    93         %f%.ForSparseJac(%q%, %r%)
    94 %$$
    95 Note that the memory required for the calculation is proportional
     93        %f%.ForSparseJac(%q%, %r%, %packed%)
     94%$$
     95Note that if $italic packed$$ was true,
     96the memory required for the calculation is proportional
    9697to $latex q$$ times the total number of variables
    9798in the AD operation sequence corresponding to $italic f$$
    9899($xref/SeqProperty/size_var/f.size_var/$$).
     100
     101$subhead packed$$
     102If $italic packed$$ was true in the call to $code ForSparseJac$$,
     103during the sparsity calculation sets of indices are represented
     104as vectors of bits that packed into words and operations are done
     105on multiple bits at a time (the number of bits in a word is unspecified).
     106Otherwise, sets of indices are represents using a sparse structure
     107that only includes the non-zero indices and operations are done
     108one index at a time.
     109$pre
     110
     111$$
     112The default value for $italic packed$$ is true; i.e.,
     113the value used if it is not present.
    99114
    100115$head r$$
     
    194209namespace CppAD {
    195210
    196 template <class Base>
    197 template <class VectorBool>
    198 VectorBool ADFun<Base>::RevSparseHes(size_t q,  const VectorBool &s)
     211template <class Base, class VectorBool, class VectorSet>
     212void RevSparseHes(
     213        size_t                    q                 ,
     214        const VectorBool&         s                 ,
     215        VectorBool&               h                 ,
     216        size_t                    total_num_var     ,
     217        CppAD::vector<size_t>&    dep_taddr         ,
     218        CppAD::vector<size_t>&    ind_taddr         ,
     219        CppAD::player<Base>&      play              ,
     220        VectorSet&                for_jac_sparsity  )
    199221{
    200222        // temporary indices
     
    205227
    206228        // range and domain dimensions for F
    207         size_t m = dep_taddr_.size();
    208         size_t n = ind_taddr_.size();
     229        size_t m = dep_taddr.size();
     230        size_t n = ind_taddr.size();
    209231
    210232        CPPAD_ASSERT_KNOWN(
    211                 q == for_jac_sparsity_.limit(),
     233                q == for_jac_sparsity.limit(),
    212234                "RevSparseHes: q (first argument) is not equal to its value"
    213235                " in the previous call to ForSparseJac with this ADFun object."
     
    219241        );
    220242
    221         // array that will hold packed reverse Jacobian values
     243        // Array that will hold reverse Jacobian dependency flag.
     244        // Initialize as true for the dependent variables.
    222245        bool *RevJac = CPPAD_NULL;
    223         RevJac       = CPPAD_TRACK_NEW_VEC(total_num_var_, RevJac);     
     246        RevJac       = CPPAD_TRACK_NEW_VEC(total_num_var, RevJac);     
     247        for(i = 0; i < total_num_var; i++)
     248                RevJac[i] = false;
     249        for(i = 0; i < m; i++)
     250        {       CPPAD_ASSERT_UNKNOWN( dep_taddr[i] < total_num_var );
     251                RevJac[ dep_taddr[i] ] = s[i];
     252        }
     253
    224254
    225255        // vector of sets that will hold packed reverse Hessain values
    226         vector_pack      rev_hes_sparsity;
    227         rev_hes_sparsity.resize(total_num_var_, q);
    228 
    229         // initialize RevJac matrix to false
    230         for(i = 0; i < total_num_var_; i++)
    231                 RevJac[i] = false;
    232 
    233         for(i = 0; i < m; i++)
    234         {       CPPAD_ASSERT_UNKNOWN( dep_taddr_[i] < total_num_var_ );
    235                 RevJac[ dep_taddr_[i] ] = s[i];
    236         }
     256        VectorSet      rev_hes_sparsity;
     257        rev_hes_sparsity.resize(total_num_var, q);
    237258
    238259        // compute the Hessian sparsity patterns
    239260        RevHesSweep(
    240261                n,
    241                 total_num_var_,
    242                 &play_,
    243                 for_jac_sparsity_,
     262                total_num_var,
     263                &play,
     264                for_jac_sparsity,
    244265                RevJac,
    245266                rev_hes_sparsity
     
    247268
    248269        // return values corresponding to independent variables
    249         VectorBool h(n * q);
    250 
    251         // j is index corresponding to reverse mode martial
     270        CPPAD_ASSERT_UNKNOWN( h.size() == n * q );
     271
     272        // j is index corresponding to reverse mode partial
    252273        for(j = 0; j < n; j++)
    253         {       CPPAD_ASSERT_UNKNOWN( ind_taddr_[j] < total_num_var_ );
    254 
    255                 // ind_taddr_[j] is operator taddr for j-th independent variable
    256                 CPPAD_ASSERT_UNKNOWN( ind_taddr_[j] == j + 1 );
    257                 CPPAD_ASSERT_UNKNOWN( play_.GetOp( ind_taddr_[j] ) == InvOp );
    258 
    259                 // i is index corresponding to forward mode partial
     274        {       CPPAD_ASSERT_UNKNOWN( ind_taddr[j] < total_num_var );
     275
     276                // ind_taddr[j] is operator taddr for j-th independent variable
     277                CPPAD_ASSERT_UNKNOWN( ind_taddr[j] == j + 1 );
     278                CPPAD_ASSERT_UNKNOWN( play.GetOp( ind_taddr[j] ) == InvOp );
     279
     280                // set all bits false
    260281                for(i = 0; i < q; i++)
    261282                        h[ i * n + j ] = false;
    262283
     284                // set bits that are true
    263285                CPPAD_ASSERT_UNKNOWN( rev_hes_sparsity.limit() == q );
    264286                i = rev_hes_sparsity.next_element(j + 1);
     
    272294        CPPAD_TRACK_DEL_VEC(RevJac);
    273295
     296        return;
     297}
     298
     299template <class Base>
     300template <class VectorBool>
     301VectorBool ADFun<Base>::RevSparseHes(size_t q,  const VectorBool &s)
     302{       size_t n = ind_taddr_.size();   
     303        VectorBool h( n * q );
     304
     305        if( for_jac_sparsity_.n_set() > 0 )
     306        {       CPPAD_ASSERT_UNKNOWN( for_jac_sparse_set_.n_set() == 0 );
     307                // use vector_pack for the calculation
     308                CppAD::RevSparseHes(
     309                        q                 ,
     310                        s                 ,
     311                        h                 ,
     312                        total_num_var_    ,
     313                        dep_taddr_        ,
     314                        ind_taddr_        ,
     315                        play_             ,
     316                        for_jac_sparsity_
     317                );
     318        }
     319        else
     320        {       CPPAD_ASSERT_UNKNOWN( for_jac_sparsity_.n_set() == 0 );
     321                // use vector_pack for the calculation
     322                CppAD::RevSparseHes(
     323                        q                    ,
     324                        s                    ,
     325                        h                    ,
     326                        total_num_var_       ,
     327                        dep_taddr_           ,
     328                        ind_taddr_           ,
     329                        play_                ,
     330                        for_jac_sparse_set_
     331                );
     332        }
    274333        return h;
    275334}
  • trunk/cppad/local/rev_sparse_jac.hpp

    r1532 r1533  
    172172template <class Base, class VectorBool, class VectorSet>
    173173void ForSparseJac(
    174         size_t                 total_num_var    ,
    175174        size_t                 p                ,
    176175        const VectorBool&      s                ,
     176        VectorBool&            r                ,
     177        size_t                 total_num_var    ,
    177178        CppAD::vector<size_t>& dep_taddr        ,
    178179        CppAD::vector<size_t>& ind_taddr        ,
    179         CppAD::player<Base>&   play             ,
    180         VectorBool&            r                )
     180        CppAD::player<Base>&   play             )
    181181{
    182182        // temporary indices
     
    253253        if( packed )
    254254        {       CppAD::ForSparseJac<Base, VectorBool, vector_pack>(
    255                         total_num_var_   ,
    256255                        p                ,
    257256                        s                ,
     257                        r                ,
     258                        total_num_var_   ,
    258259                        dep_taddr_       ,
    259260                        ind_taddr_       ,
    260                         play_            ,
    261                         r
     261                        play_
    262262                );
    263263        }
    264264        else
    265265        {       CppAD::ForSparseJac<Base, VectorBool, vector_set>(
    266                         total_num_var_   ,
    267266                        p                ,
    268267                        s                ,
     268                        r                ,
     269                        total_num_var_   ,
    269270                        dep_taddr_       ,
    270271                        ind_taddr_       ,
    271                         play_            ,
    272                         r
     272                        play_
    273273                );
    274274        }
  • trunk/omh/whats_new_09.omh

    r1532 r1533  
    5656trying to read and understand the CppAD source code.)
    5757
     58$head 09-28$$
     59Changed $cref/RevSparseHes/$$ so that it uses a sparse
     60representation when the corresponding call to
     61$cref/ForSparseJac/$$ used a sparse representation.
     62This should have been included with the change on 09-26
     63because Hessian sparsity patters after $code ForSparseJac$$
     64with $icode packed$$ did not work.
     65Thus, this could be considered a bug fix.
     66
    5867$head 09-26$$
    5968Added the $cref/packed/ForSparseJac/packed/$$ parameter to
  • trunk/test_more/rev_sparse_hes.cpp

    r1521 r1533  
    1616namespace { // Begin empty namespace
    1717
    18 bool case_one(void)
     18bool case_one(bool packed)
    1919{       bool ok = true;
    2020        using namespace CppAD;
     
    8080
    8181        // compute sparsity pattern for Jacobian of F(U(x))
    82         F.ForSparseJac(n, Px);
     82        F.ForSparseJac(n, Px, packed);
    8383
    8484        // compute sparsity pattern for Hessian of F_0 ( U(x) )
     
    103103}
    104104
    105 bool case_two(void)
     105bool case_two(bool packed)
    106106{       bool ok = true;
    107107        using namespace CppAD;
     
    152152
    153153        // compute sparsity pattern for Jacobian of F(U(x))
    154         F.ForSparseJac(n, Px);
     154        F.ForSparseJac(n, Px, packed);
    155155
    156156        // compute sparsity pattern for Hessian of F_0 ( U(x) )
     
    167167}
    168168
    169 bool case_three(void)
     169bool case_three(bool packed)
    170170{       bool ok = true;
    171171        using CppAD::AD;
     
    200200
    201201        // compute sparsity pattern for J(x) = F^{(1)} (x)
    202         f.ForSparseJac(n, r);
     202        f.ForSparseJac(n, r, packed);
    203203
    204204        // compute sparsity pattern for H(x) = F_0^{(2)} (x)
     
    216216}
    217217
    218 bool case_four(void)
     218bool case_four(bool packed)
    219219{       bool ok = true;
    220220        using namespace CppAD;
     
    276276                r[ i * n + i ] = true;
    277277        }
    278         F.ForSparseJac(n, r);
     278        F.ForSparseJac(n, r, packed);
    279279
    280280        // compute the reverse Hessian sparsity pattern for F
     
    296296bool rev_sparse_hes(void)
    297297{       bool ok = true;
    298         ok &= case_one();
    299         ok &= case_two();
    300         ok &= case_three();
    301         ok &= case_four();
    302 
    303         return ok;
    304 }
     298
     299        ok &= case_one(true);
     300        ok &= case_two(true);
     301        ok &= case_three(true);
     302        ok &= case_four(true);
     303
     304        ok &= case_one(false);
     305        ok &= case_two(false);
     306        ok &= case_three(false);
     307        ok &= case_four(false);
     308
     309        return ok;
     310}
Note: See TracChangeset for help on using the changeset viewer.