Changeset 1532


Ignore:
Timestamp:
Sep 28, 2009 7:55:30 AM (11 years ago)
Author:
bradbell
Message:

trunk: Add option for sparse Jacobians to used vector_set for calculations.

for_sparse_jac.cpp: test forward Jacobian sparsity using vector_set.
rev_sparse_jac.cpp: test reverse Jacobian sparsity using vector_set.
whats_new_09.omh: user's view of the changes.
ad_fun.hpp: add for_jac_sparse_set to private data (move private to front).
for_sparse_jac.hpp: add packed (true or false) option.
rev_sparse_jac.hpp: add packed (true or false) option.

Location:
trunk
Files:
6 edited

Legend:

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

    r1531 r1532  
    6565template <class Base>
    6666class ADFun {
    67 
     67public:
     68// ------------------------------------------------------------
     69private:
     70
     71        // debug checking number of comparision operations that changed
     72        size_t compare_change_;
     73
     74        // number of taylor_ coefficieint per variable (currently stored)
     75        size_t taylor_per_var_;
     76
     77        // number of columns currently allocated for taylor_ array
     78        size_t taylor_col_dim_;
     79
     80        // number of rows (variables) in the recording (play_)
     81        size_t total_num_var_;
     82
     83        // tape address for the independent variables
     84        CppAD::vector<size_t> ind_taddr_;
     85
     86        // tape address and parameter flag for the dependent variables
     87        CppAD::vector<size_t> dep_taddr_;
     88
     89        // which dependent variables are actually parameters
     90        CppAD::vector<bool>   dep_parameter_;
     91
     92        // the operation sequence corresponding to this object
     93        player<Base> play_;
     94
     95        // results of the forward mode calculations
     96        Base *taylor_;
     97
     98        // results of the forward mode Jacobian sparsity calculations
     99        // (n_set() for one of these two will always be zero)
     100        vector_pack      for_jac_sparsity_;
     101        vector_set       for_jac_sparse_set_;
     102
     103        // change the operation sequence corresponding to this object
     104        template <typename ADvector>
     105        void Dependent(ADTape<Base> *tape, const ADvector &y);
     106
     107// ------------------------------------------------------------
    68108public:
    69109        // default constructor
     
    98138
    99139        // forward mode Jacobian sparsity
    100         template <typename VectorBase>
    101         VectorBase ForSparseJac(size_t q, const VectorBase &Px);
    102 
     140        template <typename VectorBool>
     141        VectorBool ForSparseJac(
     142                size_t q, const VectorBool &Px, bool packed = true
     143        );
    103144        // reverse mode Jacobian sparsity
    104145        template <typename VectorBool>
    105         VectorBool RevSparseJac(size_t q, const VectorBool &Py);
    106 
     146        VectorBool RevSparseJac(
     147                size_t q, const VectorBool &Py, bool packed = true
     148        );
    107149        // reverse mode Hessian sparsity
    108150        template <typename VectorBool>
     
    230272        size_t taylor_size(void) const
    231273        {       return taylor_per_var_; }
    232         // ------------------------------------------------------------
    233 private:
    234         // maximum amount of memory required for this function object
    235         // mutable size_t memoryMax;
    236 
    237         // debug checking number of comparision operations that changed
    238         size_t compare_change_;
    239 
    240         // number of taylor_ coefficieint per variable (currently stored)
    241         size_t taylor_per_var_;
    242 
    243         // number of columns currently allocated for taylor_ array
    244         size_t taylor_col_dim_;
    245 
    246         // number of rows (variables) in the recording (play_)
    247         size_t total_num_var_;
    248 
    249         // tape address for the independent variables
    250         CppAD::vector<size_t> ind_taddr_;
    251 
    252         // tape address and parameter flag for the dependent variables
    253         CppAD::vector<size_t> dep_taddr_;
    254         CppAD::vector<bool>   dep_parameter_;
    255 
    256         // the operations corresponding to this function
    257         player<Base> play_;
    258 
    259         // results of the forward mode calculations
    260         Base *taylor_;
    261 
    262         // results of the forward mode Jacobian sparsity calculations
    263         vector_pack      for_jac_sparsity_;
    264 
    265         template <typename ADvector>
    266         void Dependent(ADTape<Base> *tape, const ADvector &y);
    267274};
    268275// ---------------------------------------------------------------------------
  • trunk/cppad/local/for_sparse_jac.hpp

    r1531 r1532  
    3535$head Syntax$$
    3636$syntax%%s% = %f%.ForSparseJac(%q%, %r%)%$$
    37 
     37$pre
     38$$
     39$syntax%%s% = %f%.ForSparseJac(%q%, %r%, %packed%)%$$
    3840
    3941$head Purpose$$
     
    9597        R_{i,j} \neq 0 ; \Rightarrow \; r [ i * q + j ] = {\rm true}
    9698\] $$
     99
     100$head packed$$
     101If $italic packed$$ is true,
     102during the sparsity calculation sets of indices are represented
     103as vectors of bits that packed into words and operations are done
     104on multiple bits at a time (the number of bits in a word is unspecified).
     105Otherwise, sets of indices are represents using a sparse structure
     106that only includes the non-zero indices and operations are done
     107one index at a time.
     108$pre
     109
     110$$
     111The default value for $italic packed$$ is true; i.e.,
     112the value used if it is not present.
    97113
    98114$head s$$
     
    155171namespace CppAD {
    156172
    157 template <class Base>
    158 template <class VectorBool>
    159 VectorBool ADFun<Base>::ForSparseJac(size_t q, const VectorBool &r)
     173template <class Base, class VectorBool, class VectorSet>
     174void ForSparseJac(
     175        size_t                 total_num_var    ,
     176        size_t                 q                ,
     177        const VectorBool&      r                ,
     178        CppAD::vector<size_t>& dep_taddr        ,
     179        CppAD::vector<size_t>& ind_taddr        ,
     180        CppAD::player<Base>&   play             ,
     181        VectorBool&            s                ,
     182        VectorSet&             for_jac_sparsity )
    160183{
    161184        // temporary indices
     
    166189
    167190        // range and domain dimensions for F
    168         size_t m = dep_taddr_.size();
    169         size_t n = ind_taddr_.size();
     191        size_t m = dep_taddr.size();
     192        size_t n = ind_taddr.size();
    170193
    171194        CPPAD_ASSERT_KNOWN(
     
    181204
    182205        // allocate memory for the requested sparsity calculation
    183         for_jac_sparsity_.resize(total_num_var_, q);
     206        for_jac_sparsity.resize(total_num_var, q);
    184207
    185208        // set values corresponding to independent variables
    186209        for(i = 0; i < n; i++)
    187         {       CPPAD_ASSERT_UNKNOWN( ind_taddr_[i] < total_num_var_ );
    188                 // ind_taddr_[i] is operator taddr for i-th independent variable
    189                 CPPAD_ASSERT_UNKNOWN( play_.GetOp( ind_taddr_[i] ) == InvOp );
     210        {       CPPAD_ASSERT_UNKNOWN( ind_taddr[i] < total_num_var );
     211                // ind_taddr[i] is operator taddr for i-th independent variable
     212                CPPAD_ASSERT_UNKNOWN( play.GetOp( ind_taddr[i] ) == InvOp );
    190213
    191214                // set bits that are true
    192215                for(j = 0; j < q; j++) if( r[ i * q + j ] )
    193                         for_jac_sparsity_.add_element( ind_taddr_[i], j);
     216                        for_jac_sparsity.add_element( ind_taddr[i], j);
    194217        }
    195218
     
    197220        ForJacSweep(
    198221                n,
    199                 total_num_var_,
    200                 &play_,
    201                 for_jac_sparsity_
     222                total_num_var,
     223                &play,
     224                for_jac_sparsity
    202225        );
    203226
    204227        // return values corresponding to dependent variables
    205         VectorBool s(m * q);
     228        CPPAD_ASSERT_UNKNOWN( s.size() == m * q );
    206229        for(i = 0; i < m; i++)
    207         {       CPPAD_ASSERT_UNKNOWN( dep_taddr_[i] < total_num_var_ );
     230        {       CPPAD_ASSERT_UNKNOWN( dep_taddr[i] < total_num_var );
    208231
    209232                // set bits
    210233                for(j = 0; j < q; j++)
    211234                        s[ i * q + j ] = false;
    212                 CPPAD_ASSERT_UNKNOWN( for_jac_sparsity_.limit() == q );
    213                 j = for_jac_sparsity_.next_element( dep_taddr_[i] );
     235                CPPAD_ASSERT_UNKNOWN( for_jac_sparsity.limit() == q );
     236                j = for_jac_sparsity.next_element( dep_taddr[i] );
    214237                while( j < q )
    215238                {       s[i * q + j ] = true;
    216                         j = for_jac_sparsity_.next_element( dep_taddr_[i] );
     239                        j = for_jac_sparsity.next_element( dep_taddr[i] );
    217240                }
    218241        }
    219 
     242}
     243
     244template <class Base>
     245template <class VectorBool>
     246VectorBool ADFun<Base>::ForSparseJac(
     247        size_t             q      ,
     248        const VectorBool&  r      ,
     249        bool               packed )
     250{       size_t m = dep_taddr_.size();
     251        VectorBool s( m * q );
     252
     253        if( packed )
     254        {       // free any results stored in for_jac_sparse_set_       
     255                for_jac_sparse_set_.resize(0, 0);
     256                // store results in for_jac_sparsity_
     257                CppAD::ForSparseJac(
     258                        total_num_var_   ,
     259                        q                ,
     260                        r                ,
     261                        dep_taddr_       ,
     262                        ind_taddr_       ,
     263                        play_            ,
     264                        s                ,
     265                        for_jac_sparsity_
     266                );
     267        }
     268        else
     269        {       // free any results stored in for_jac_sparsity_
     270                for_jac_sparsity_.resize(0, 0);
     271                // store results in for_jac_sparse_set_
     272                CppAD::ForSparseJac(
     273                        total_num_var_   ,
     274                        q                ,
     275                        r                ,
     276                        dep_taddr_       ,
     277                        ind_taddr_       ,
     278                        play_            ,
     279                        s                ,
     280                        for_jac_sparse_set_
     281                );
     282        }
    220283        return s;
    221284}
  • trunk/cppad/local/rev_sparse_jac.hpp

    r1531 r1532  
    3636$head Syntax$$
    3737$syntax%%r% = %F%.RevSparseJac(%p%, %s%)%$$
     38$pre
     39$$
     40$syntax%%r% = %F%.RevSparseJac(%p%, %s%, %packed%)%$$
    3841
    3942
     
    9194        S_{i,j} \neq 0 ; \Rightarrow \; s [ i * m + j ] = {\rm true}
    9295\] $$
     96
     97$head packed$$
     98If $italic packed$$ is true,
     99during the sparsity calculation sets of indices are represented
     100as vectors of bits that packed into words and operations are done
     101on multiple bits at a time (the number of bits in a word is unspecified).
     102Otherwise, sets of indices are represents using a sparse structure
     103that only includes the non-zero indices and operations are done
     104one index at a time.
     105$pre
     106
     107$$
     108The default value for $italic packed$$ is true; i.e.,
     109the value used if it is not present.
    93110
    94111$head r$$
     
    153170namespace CppAD {
    154171
    155 template <class Base>
    156 template <class VectorBool>
    157 VectorBool ADFun<Base>::RevSparseJac(size_t p, const VectorBool &s)
     172template <class Base, class VectorBool, class VectorSet>
     173void ForSparseJac(
     174        size_t                 total_num_var    ,
     175        size_t                 p                ,
     176        const VectorBool&      s                ,
     177        CppAD::vector<size_t>& dep_taddr        ,
     178        CppAD::vector<size_t>& ind_taddr        ,
     179        CppAD::player<Base>&   play             ,
     180        VectorBool&            r                )
    158181{
    159182        // temporary indices
     
    164187
    165188        // range and domain dimensions for F
    166         size_t m = dep_taddr_.size();
    167         size_t n = ind_taddr_.size();
     189        size_t m = dep_taddr.size();
     190        size_t n = ind_taddr.size();
    168191
    169192        CPPAD_ASSERT_KNOWN(
     
    179202
    180203        // vector of sets that will hold the results
    181         vector_pack      var_sparsity;
    182         var_sparsity.resize(total_num_var_, p);
     204        VectorSet      var_sparsity;
     205        var_sparsity.resize(total_num_var, p);
    183206
    184207        // The sparsity pattern corresponding to the dependent variables
    185208        for(i = 0; i < m; i++)
    186         {       CPPAD_ASSERT_UNKNOWN( dep_taddr_[i] < total_num_var_ );
     209        {       CPPAD_ASSERT_UNKNOWN( dep_taddr[i] < total_num_var );
    187210
    188211                for(j = 0; j < p; j++) if( s[ i * m + j ] )
    189                         var_sparsity.add_element( dep_taddr_[i], j );
     212                        var_sparsity.add_element( dep_taddr[i], j );
    190213        }
    191214
     
    193216        RevJacSweep(
    194217                n,
    195                 total_num_var_,
    196                 &play_,
     218                total_num_var,
     219                &play,
    197220                var_sparsity
    198221        );
    199222
    200223        // return values corresponding to dependent variables
    201         VectorBool r(p * n);
     224        CPPAD_ASSERT_UNKNOWN( r.size() == p * n );
    202225        for(j = 0; j < n; j++)
    203         {       CPPAD_ASSERT_UNKNOWN( ind_taddr_[j] == (j+1) );
    204 
    205                 // ind_taddr_[j] is operator taddr for j-th independent variable
    206                 CPPAD_ASSERT_UNKNOWN( play_.GetOp( ind_taddr_[j] ) == InvOp );
     226        {       CPPAD_ASSERT_UNKNOWN( ind_taddr[j] == (j+1) );
     227
     228                // ind_taddr[j] is operator taddr for j-th independent variable
     229                CPPAD_ASSERT_UNKNOWN( play.GetOp( ind_taddr[j] ) == InvOp );
    207230
    208231                // set bits
     
    217240                }
    218241        }
    219 
     242}
     243
     244template <class Base>
     245template <class VectorBool>
     246VectorBool ADFun<Base>::RevSparseJac(
     247        size_t              p      ,
     248        const VectorBool&   s      ,
     249        bool                packed )
     250{       size_t n = ind_taddr_.size();
     251        VectorBool r( p * n );
     252
     253        if( packed )
     254        {       CppAD::ForSparseJac<Base, VectorBool, vector_pack>(
     255                        total_num_var_   ,
     256                        p                ,
     257                        s                ,
     258                        dep_taddr_       ,
     259                        ind_taddr_       ,
     260                        play_            ,
     261                        r
     262                );
     263        }
     264        else
     265        {       CppAD::ForSparseJac<Base, VectorBool, vector_set>(
     266                        total_num_var_   ,
     267                        p                ,
     268                        s                ,
     269                        dep_taddr_       ,
     270                        ind_taddr_       ,
     271                        play_            ,
     272                        r
     273                );
     274        }
    220275        return r;
    221276}
  • trunk/omh/whats_new_09.omh

    r1521 r1532  
    4141        cond_exp
    4242        VecAD
     43        Jacobians
     44        Jac
    4345$$
    4446
     
    5355(Comments about developer documentation are only important if you are
    5456trying to read and understand the CppAD source code.)
     57
     58$head 09-26$$
     59Added the $cref/packed/ForSparseJac/packed/$$ parameter to
     60$cref/ForSparseJac/$$ and $cref/RevSparseJac/$$.
     61If $icode packed$$ is false,
     62a sparse instead of packed representation is used
     63during the calculations of sparsity patterns.
     64The sparse representation
     65should be faster, and use less memory, for very large sparse Jacobians.
     66The functions $code ForSparseJac$$ and $code RevSparseJac$$
     67return packed representations.
     68The plan is to eventually provide new member functions that return
     69sparse representations.
    5570
    5671$head 09-20$$
  • trunk/test_more/for_sparse_jac.cpp

    r1475 r1532  
    7171namespace { // Begin empty namespace
    7272
    73 bool case_one(void)
     73bool case_one(bool packed)
    7474{       bool ok = true;
    7575        using namespace CppAD;
     
    160160        // evaluate the dependency matrix for F(X(x))
    161161        CPPAD_TEST_VECTOR< bool > Py(m * n);
    162         Py = F.ForSparseJac(n, Px);
     162        Py = F.ForSparseJac(n, Px, packed);
    163163
    164164        // check values
     
    171171}
    172172
    173 bool case_two(void)
     173bool case_two(bool packed)
    174174{       bool ok = true;
    175175        using namespace CppAD;
     
    246246        // evaluate the dependency matrix for F(X(x))
    247247        CPPAD_TEST_VECTOR< bool > Py(m * n);
    248         Py = F.ForSparseJac(n, Px);
     248        Py = F.ForSparseJac(n, Px, packed);
    249249
    250250        // check values
     
    261261bool for_sparse_jac(void)
    262262{       bool ok = true;
    263         ok &= case_one();
    264         ok &= case_two();
     263
     264        ok &= case_one(true);
     265        ok &= case_two(true);
     266
     267        ok &= case_one(false);
     268        ok &= case_two(false);
     269
    265270        return ok;
    266271}
  • trunk/test_more/rev_sparse_jac.cpp

    r1521 r1532  
    7070namespace { // BEGIN empty namespace
    7171
    72 bool case_one(void)
     72bool case_one(bool packed)
    7373{       bool ok = true;
    7474        using namespace CppAD;
     
    160160        // evaluate the dependency matrix for U(F(x))
    161161        CPPAD_TEST_VECTOR< bool > Px(m * n);
    162         Px = F.RevSparseJac(m, Py);
     162        Px = F.RevSparseJac(m, Py, packed);
    163163
    164164        // check values
     
    171171}
    172172
    173 bool case_two(void)
     173bool case_two(bool packed)
    174174{       bool ok = true;
    175175        using namespace CppAD;
     
    246246        // evaluate the dependency matrix for S [ F(x) ]
    247247        CPPAD_TEST_VECTOR< bool > Px(m * n);
    248         Px = F.RevSparseJac(m, Py);
     248        Px = F.RevSparseJac(m, Py, packed);
    249249
    250250        // check values
     
    262262{       bool ok = true;
    263263
    264         ok &= case_one();
    265         ok &= case_two();
     264        ok &= case_one(true);
     265        ok &= case_two(true);
     266
     267        ok &= case_one(false);
     268        ok &= case_two(false);
    266269
    267270        return ok;
Note: See TracChangeset for help on using the changeset viewer.