Changeset 1551


Ignore:
Timestamp:
Oct 16, 2009 11:23:06 AM (11 years ago)
Author:
bradbell
Message:

trunk: Extend SparseJacobian? to use vector of sets for sparsity.

check_doxygen.sh: add SparseJacobian? to check for warnings.
*/sparse_jacobian.cpp: examples and tests of sparse jacobian calculations.
makefile.am: add makefile entry for old example, now in test_more.
makefile.in: change generated automatically by change in makefile.am.
jacobian.cpp: more testing of sparse jacobians.
test_more.cpp: add new entry for sparse_jacobian.cpp.
whats_new_09.omh: user's view of the changes.
ad_fun.hpp: new private helper functions to split cases.
for_sparse_jac.hpp: minor clean up.
rev_sparse_jac.hpp: minor clean up.
rev_sparse_hes.hpp: minor clean up.
sparse_jacobian.hpp: extend to handle and use vector of sets sparsity patterns.

Location:
trunk
Files:
1 added
11 edited

Legend:

Unmodified
Added
Removed
  • trunk/check_doxygen.sh

    r1539 r1551  
    1818        sparse_pack
    1919        sparse_set
     20"
     21#
     22# member functions names that begin with the following have been documented
     23# (but the corresponding class has not yet been completely documented).
     24member_list="
     25        SparseJacobian
    2026"
    2127# files that have been completely documented
     
    6066        store_op.hpp
    6167        sub_op.hpp
     68        sparse_jacobian.hpp
    6269        sparse_pack.hpp
    6370        sparse_set.hpp
     
    6774do
    6875        if grep -i "warning:.*$class::" doxygen.log
     76        then
     77                echo "Unexpected doxygen error or warning for $file."
     78                exit 1
     79        fi
     80done
     81# --------------------------------------------------------------------------
     82for member in $member_list
     83do
     84        if grep -i "warning:.*member.*$member" doxygen.log
    6985        then
    7086                echo "Unexpected doxygen error or warning for $file."
  • trunk/cppad/local/ad_fun.hpp

    r1550 r1551  
    180180                VectorSet&               h
    181181        );
     182        // ------------------------------------------------------------
     183        // vector of bool version of SparseJacobian
     184        // (see doxygen in sparse_jacobian.hpp)
     185        template <class VectorBase, class VectorSet>
     186        void SparseJacobianCase(
     187                bool                     set_type    ,
     188                const VectorBase&        x           ,
     189                const VectorSet&         p           ,
     190                VectorBase&              jac
     191        );
     192        // vector of std::set<size_t> version of SparseJacobian
     193        // (see doxygen in sparse_jacobian.hpp)
     194        template <class VectorBase, class VectorSet>
     195        void SparseJacobianCase(
     196                const std::set<size_t>&  set_type    ,
     197                const VectorBase&        x           ,
     198                const VectorSet&         p           ,
     199                VectorBase&              jac
     200        );
    182201// ------------------------------------------------------------
    183202public:
     
    318337
    319338        /// calculate sparse Jacobians
    320         template <typename VectorBase, typename VectorBool>
    321         VectorBase SparseJacobian(const VectorBase &x, const VectorBool &p);
     339        template <typename VectorBase, typename VectorSet>
     340        VectorBase SparseJacobian(const VectorBase &x, const VectorSet &p);
    322341
    323342        /// calculate sparse Hessians
  • trunk/cppad/local/for_sparse_jac.hpp

    r1550 r1551  
    3838
    3939$head Purpose$$
    40 We use $latex F : \R^n \rightarrow B^m$$ to denote the
     40We use $latex F : \R^n \rightarrow R^m$$ to denote the
    4141$xref/glossary/AD Function/AD function/$$ corresponding to $icode f$$.
    4242For a fixed $latex n \times q$$ matrix $latex R$$,
     
    109109The type $icode VectorSet$$ must be a $xref/SimpleVector/$$ class with
    110110$xref/SimpleVector/Elements of Specified Type/elements of type/$$
    111 $code bool$$ or $code std::set<size_t>$$.
     111$code bool$$ or $code std::set<size_t>$$;
     112see $cref/sparsity pattern/glossary/Sparsity Pattern/$$ for a discussion
     113of the difference.
    112114
    113115$head Entire Sparsity Pattern$$
     
    492494\param set_type
    493495is a \c bool value. This argument is used to dispatch to the proper source
    494 code depending on the value of \c VectroSet::value_type.
     496code depending on the value of \c VectorSet::value_type.
    495497
    496498\param q
     
    511513        const VectorSet&    r             ,
    512514        VectorSet&          s             )
    513 {       size_t m = dep_taddr_.size();
     515{       size_t m = Range();
    514516
    515517        // dimension size of result vector
     
    560562        const VectorSet&           r             ,
    561563        VectorSet&                 s             )
    562 {       size_t m = dep_taddr_.size();
     564{       size_t m = Range();
    563565
    564566        // dimension size of result vector
  • trunk/cppad/local/rev_sparse_hes.hpp

    r1550 r1551  
    124124The type $icode VectorSet$$ must be a $cref/SimpleVector/$$ class with
    125125$xref/SimpleVector/Elements of Specified Type/elements of type/$$
    126 $code bool$$ or $code std::set<size_t>$$.
     126$code bool$$ or $code std::set<size_t>$$;
     127see $cref/sparsity pattern/glossary/Sparsity Pattern/$$ for a discussion
     128of the difference.
    127129The type of the elements of
    128130$cref/VectorSet/RevSparseHes/VectorSet/$$ must be the
     
    555557        const VectorSet&  s                ,
    556558        VectorSet&        h                )
    557 {       size_t n = ind_taddr_.size();   
     559{       size_t n = Domain();   
    558560        h.resize(q * n );
    559561
  • trunk/cppad/local/rev_sparse_jac.hpp

    r1550 r1551  
    4343
    4444$head Purpose$$
    45 We use $latex F : \R^n \rightarrow B^m$$ to denote the
     45We use $latex F : \R^n \rightarrow R^m$$ to denote the
    4646$xref/glossary/AD Function/AD function/$$ corresponding to $icode f$$.
    4747For a fixed $latex p \times m$$ matrix $latex S$$,
     
    108108The type $icode VectorSet$$ must be a $xref/SimpleVector/$$ class with
    109109$xref/SimpleVector/Elements of Specified Type/elements of type/$$
    110 $code bool$$ or $code std::set<size_t>$$.
     110$code bool$$ or $code std::set<size_t>$$;
     111see $cref/sparsity pattern/glossary/Sparsity Pattern/$$ for a discussion
     112of the difference.
    111113
    112114$head Entire Sparsity Pattern$$
     
    468470        const VectorSet&    s                 ,
    469471        VectorSet&          r                 )
    470 {       size_t n = ind_taddr_.size();
     472{       size_t n = Domain();
    471473
    472474        // dimension of the result vector
  • trunk/cppad/local/sparse_jacobian.hpp

    r1531 r1551  
    1717$begin sparse_jacobian$$
    1818$spell
     19        std
     20        CppAD
    1921        Bool
    2022        jac
     
    3436
    3537$head Purpose$$
    36 We use $latex F : B^n \rightarrow B^m$$ do denote the
     38We use $latex F : \R^n \rightarrow \R^m$$ do denote the
    3739$cref/AD function/glossary/AD Function/$$
    3840corresponding to $icode f$$.
     
    4143        jac = F^{(1)} (x)
    4244\] $$
    43 This is a preliminary implementation of a method for using the fact
    44 that the matrix is sparse to reduce the amount of computation necessary.
    45 One should use speed tests to verify that results are computed faster
     45This routine assumes
     46that the matrix $latex F^{(1)} (x) \in \R^{m \times n}$$
     47is uses this assumption to reduce the amount of computation necessary.
     48One should use speed tests (e.g. $cref/speed_test/$$)
     49to verify that results are computed faster
    4650than when using the routine $cref/Jacobian/$$.
    4751
     
    5761The argument $icode x$$ has prototype
    5862$codei%
    59         const %VectorBase% &%x%
     63        const %VectorBase%& %x%
    6064%$$
    6165(see $cref/VectorBase/sparse_jacobian/VectorBase/$$ below)
     
    6973The argument $icode p$$ is optional and has prototype
    7074$syntax%
    71         const %VectorBool% &%p%
     75        const %VectorSet%& %p%
    7276%$$
    73 (see $cref/VectorBool/sparse_jacobian/VectorBool/$$ below)
    74 and its size is $latex m * n$$.
     77(see $cref/VectorSet/sparse_jacobian/VectorSet/$$ below).
     78If it has elements of type $code bool$$,
     79its size is $latex m * n$$.
     80If it has elements of type $code std::set<size_t>$$,
     81its size is $latex m$$ and all its set elements are between
     82zero and $latex n - 1$$.
    7583It specifies a
    7684$cref/sparsity pattern/glossary/Sparsity Pattern/$$
    77 for the Jacobian; i.e.,
    78 for $latex i = 0 , \ldots , m-1$$ and $latex j = 0 , \ldots , n-1$$.
    79 $latex \[
    80         \D{ F_i }{ x_j } \neq 0 ; \Rightarrow \; p [ i * n + j ] = {\rm true}
    81 \] $$
     85for the Jacobian $latex F^{(1)} (x)$$.
    8286$pre
    8387
    8488$$
    8589If this sparsity pattern does not change between calls to
    86 $codei SparseJacobian$$, it should be faster to calculate $icode p$$ once and
    87 pass this argument to $codei SparseJacobian$$.
     90$codei SparseJacobian$$, it should be faster to calculate $icode p$$ once
     91(using $cref/ForSparseJac/$$ or $cref/RevSparseJac/$$)
     92and then pass $icode p$$ to $codei SparseJacobian$$.
     93In addition,
     94if you specify $icode p$$, CppAD will use the same
     95type of sparsity representation
     96(vectors of $code bool$$ or vectors of $code std::set<size_t>$$)
     97for its internal calculations.
     98Otherwise the representation
     99for the internal calculations is unspecified.
    88100
    89101$head jac$$
     
    96108and $latex j = 0 , \ldots , n - 1 $$
    97109$latex \[
    98         jac [ i * n + j ] = \D{ F_i }{ x_j }
     110        jac [ i * n + j ] = \D{ F_i }{ x_j } (x)
    99111\] $$
    100112
     
    106118if this is not the case.
    107119
    108 $head VectorBool$$
    109 The type $icode VectorBool$$ must be a $xref/SimpleVector/$$ class with
    110 $xref/SimpleVector/Elements of Specified Type/elements of type bool/$$.
    111 The routine $xref/CheckSimpleVector/$$ will generate an error message
    112 if this is not the case.
    113 In order to save memory,
    114 you may want to use a class that packs multiple elements into one
    115 storage location; for example,
    116 $cref/vectorBool/CppAD_vector/vectorBool/$$.
     120$head VectorSet$$
     121The type $icode VectorSet$$ must be a $xref/SimpleVector/$$ class with
     122$xref/SimpleVector/Elements of Specified Type/elements of type/$$
     123$code bool$$ or $code std::set<size_t>$$;
     124see $cref/sparsity pattern/glossary/Sparsity Pattern/$$ for a discussion
     125of the difference.
    117126
    118127$head Uses Forward$$
     
    135144-----------------------------------------------------------------------------
    136145*/
    137 
    138 
    139 
    140 
    141 //  BEGIN CppAD namespace
    142 namespace CppAD {
    143 
     146CPPAD_BEGIN_NAMESPACE
     147
     148/*!
     149\file sparse_jacobian.hpp
     150Computing sparse Jacobians.
     151*/
     152
     153/*!
     154Private helper function for SparseJacobian(x, p).
     155
     156All of the description in the public member function SparseJacobian(x, p)
     157applies.
     158
     159\param set_type
     160is a \c bool value. This argument is used to dispatch to the proper souce
     161code depending on the vlaue of \c VectorSet::value_type.
     162
     163\param x
     164See \c SparseJacobian(x, p).
     165
     166\param p
     167See \c SparseJacobian(x, p).
     168
     169\param jac
     170is the return value for the corresponding call to \c SparseJacobian(x, p).
     171On input it must have size equalt to the domain times range dimension
     172for this ADFun<Base> object.
     173On return, it will contain the Jacobian.
     174*/
    144175template <class Base>
    145 template <class VectorBase>
    146 VectorBase ADFun<Base>::SparseJacobian(const VectorBase &x)
    147 {       typedef CppAD::vector<bool>   VectorBool;
    148 
    149         size_t m = Range();
    150         size_t n = Domain();
    151 
    152         // sparsity pattern for Jacobian
    153         VectorBool p(n * m);
    154 
    155         if( n <= m )
    156         {       size_t j, k;
    157 
    158                 // use forward mode
    159                 VectorBool r(n * n);
    160                 for(j = 0; j < n; j++)
    161                 {       for(k = 0; k < n; k++)
    162                                 r[j * n + k] = false;
    163                         r[j * n + j] = true;
    164                 }
    165                 p = ForSparseJac(n, r);
    166         }
    167         else
    168         {       size_t i, k;
    169 
    170                 // use reverse mode
    171                 VectorBool s(m * m);
    172                 for(i = 0; i < m; i++)
    173                 {       for(k = 0; k < m; k++)
    174                                 s[i * m + k] = false;
    175                         s[i * m + i] = true;
    176                 }
    177                 p = RevSparseJac(m, s);
    178         }
    179         return SparseJacobian(x, p);
    180 }
    181 
    182 template <class Base>
    183 template <class VectorBase, class VectorBool>
    184 VectorBase ADFun<Base>::SparseJacobian(const VectorBase &x, const VectorBool &p)
     176template <class VectorBase, class VectorSet>
     177void ADFun<Base>::SparseJacobianCase(
     178        bool               set_type        ,
     179        const VectorBase&  x               ,
     180        const VectorSet&   p               ,
     181        VectorBase&        jac             )
    185182{
    186183        typedef CppAD::vector<size_t> SizeVector;
     184        typedef CppAD::vectorBool     VectorBool;
    187185        size_t i, j, k;
    188186
     
    194192        const Base one(1);
    195193
    196         // check VectorBool is Simple Vector class with bool elements
    197         CheckSimpleVector<bool, VectorBool>();
     194        // check VectorSet is Simple Vector class with bool elements
     195        CheckSimpleVector<bool, VectorSet>();
    198196
    199197        // check VectorBase is Simple Vector class with Base type elements
     
    204202                "SparseJacobian: size of x not equal domain dimension for f"
    205203        );
    206 
    207204        CPPAD_ASSERT_KNOWN(
    208205                p.size() == m * n,
    209                 "SparseJacobian: size of p not equal domain dimension times "
    210                 " range dimension for f"
     206                "SparseJacobian: size of p not equal range dimension times "
     207                " domain dimension for f"
    211208        );
     209        CPPAD_ASSERT_UNKNOWN(jac.size() == m * n);
    212210
    213211        // point at which we are evaluating the Jacobian
     
    215213
    216214        // initialize the return value
    217         VectorBase jac(m * n);
    218215        for(i = 0; i < m; i++)
    219216                for(j = 0; j < n; j++)
     
    231228                // Graph Coloring in Optimization Revisited by
    232229                // Assefaw Gebremedhin, Fredrik Maane, Alex Pothen
    233                 VectorBool    forbidden(n);
     230                VectorBool forbidden(n);
    234231                for(j = 0; j < n; j++)
    235232                {       // initial all colors as ok for this column
     
    291288                // Graph Coloring in Optimization Revisited by
    292289                // Assefaw Gebremedhin, Fredrik Maane, Alex Pothen
    293                 VectorBool    forbidden(m);
     290                VectorBool forbidden(m);
    294291                for(i = 0; i < m; i++)
    295292                {       // initial all colors as ok for this row
     
    340337                }
    341338        }
     339}
     340
     341/*!
     342Private helper function for SparseJacobian(x, p).
     343
     344All of the description in the public member function SparseJacobian(x, p)
     345applies.
     346
     347\param set_type
     348is a \c std::set<size_t> value.
     349This argument is used to dispatch to the proper souce
     350code depending on the vlaue of \c VectorSet::value_type.
     351
     352\param x
     353See \c SparseJacobian(x, p).
     354
     355\param p
     356See \c SparseJacobian(x, p).
     357
     358\param jac
     359is the return value for the corresponding call to \c SparseJacobian(x, p).
     360On input it must have size equalt to the domain times range dimension
     361for this ADFun<Base> object.
     362On return, it will contain the Jacobian.
     363*/
     364template <class Base>
     365template <class VectorBase, class VectorSet>
     366void ADFun<Base>::SparseJacobianCase(
     367        const std::set<size_t>&  set_type        ,
     368        const VectorBase&        x               ,
     369        const VectorSet&         p               ,
     370        VectorBase&              jac             )
     371{
     372        typedef CppAD::vector<size_t> SizeVector;
     373        typedef CppAD::vectorBool     VectorBool;
     374        size_t i, j, k;
     375
     376        size_t m = Range();
     377        size_t n = Domain();
     378
     379        // some values
     380        const Base zero(0);
     381        const Base one(1);
     382
     383        // cannot use CheckSimpleVector when vector elements are sets
     384        // because cannot assign an integer to a set.
     385        // CheckSimpleVector<std::set<size_t>, VectorSet>();
     386
     387        // check VectorBase is Simple Vector class with Base type elements
     388        CheckSimpleVector<Base, VectorBase>();
     389
     390        CPPAD_ASSERT_KNOWN(
     391                x.size() == n,
     392                "SparseJacobian: size of x not equal domain dimension for f"
     393        );
     394        CPPAD_ASSERT_KNOWN(
     395                p.size() == m,
     396                "SparseJacobian: size of p not equal range dimension for f"
     397        );
     398        CPPAD_ASSERT_UNKNOWN(jac.size() == m * n);
     399
     400        // point at which we are evaluating the Jacobian
     401        Forward(0, x);
     402
     403        // initialize the return value
     404        for(i = 0; i < m; i++)
     405                for(j = 0; j < n; j++)
     406                        jac[i * n + j] = zero;
     407
     408        // create a copy of the transpose sparsity pattern
     409        VectorSet q(n);
     410        std::set<size_t>::iterator itr_i, itr_j;
     411        for(i = 0; i < m; i++)
     412        {       for(itr_j = p[i].begin(); itr_j != p[i].end(); itr_j++)
     413                {       j = *itr_j;
     414                        q[j].insert(i);
     415                }
     416        }       
     417
     418        if( n <= m )
     419        {       // use forward mode ----------------------------------------
     420       
     421                // initial coloring
     422                SizeVector color(n);
     423                for(j = 0; j < n; j++)
     424                        color[j] = j;
     425
     426                // See GreedyPartialD2Coloring Algorithm Section 3.6.2 of
     427                // Graph Coloring in Optimization Revisited by
     428                // Assefaw Gebremedhin, Fredrik Maane, Alex Pothen
     429                VectorBool forbidden(n);
     430                for(j = 0; j < n; j++)
     431                {       // initial all colors as ok for this column
     432                        for(k = 0; k < n; k++)
     433                                forbidden[k] = false;
     434
     435                        // for each row connected to column j
     436                        for(itr_i = q[j].begin(); itr_i != q[j].end(); itr_i++)
     437                        {       i = *itr_i;
     438                                // for each column connected to row i
     439                                for(    itr_j = p[i].begin();
     440                                        itr_j != p[i].end();
     441                                        itr_j++
     442                                )
     443                                {       // if this is not j, forbid it
     444                                        k = *itr_j;
     445                                        forbidden[ color[k] ] = (k != j);
     446                                }
     447                        }
     448                        k = 0;
     449                        while( forbidden[k] && k < n )
     450                        {       k++;
     451                                CPPAD_ASSERT_UNKNOWN( k < n );
     452                        }
     453                        color[j] = k;
     454                }
     455                size_t n_color = 1;
     456                for(k = 0; k < n; k++)
     457                        n_color = std::max(n_color, color[k] + 1);
     458
     459                // direction vector for calls to forward
     460                VectorBase dx(n);
     461
     462                // location for return values from Reverse
     463                VectorBase dy(m);
     464
     465                // loop over colors
     466                size_t c;
     467                for(c = 0; c < n_color; c++)
     468                {       // determine all the colums with this color
     469                        for(j = 0; j < n; j++)
     470                        {       if( color[j] == c )
     471                                        dx[j] = one;
     472                                else    dx[j] = zero;
     473                        }
     474                        // call forward mode for all these columns at once
     475                        dy = Forward(1, dx);
     476
     477                        // set the corresponding components of the result
     478                        for(j = 0; j < n; j++) if( color[j] == c )
     479                        {       for(
     480                                        itr_i = q[j].begin();
     481                                        itr_i != q[j].end();
     482                                        itr_i++
     483                                )
     484                                {       i = *itr_i;
     485                                        jac[i * n + j] = dy[i];
     486                                }
     487                        }
     488                }
     489        }
     490        else
     491        {       // use reverse mode ----------------------------------------
     492       
     493                // initial coloring
     494                SizeVector color(m);
     495                for(i = 0; i < m; i++)
     496                        color[i] = i;
     497
     498                // See GreedyPartialD2Coloring Algorithm Section 3.6.2 of
     499                // Graph Coloring in Optimization Revisited by
     500                // Assefaw Gebremedhin, Fredrik Maane, Alex Pothen
     501                VectorBool forbidden(m);
     502                for(i = 0; i < m; i++)
     503                {       // initial all colors as ok for this row
     504                        for(k = 0; k < m; k++)
     505                                forbidden[k] = false;
     506
     507                        // for each column connected to row i
     508                        for(itr_j = p[i].begin(); itr_j != p[i].end(); itr_j++)
     509                        {       j = *itr_j;     
     510                                // for each row connected to column j
     511                                for(    itr_i = q[j].begin();
     512                                        itr_i != q[j].end();
     513                                        itr_i++
     514                                )
     515                                {       // if this is not i, forbid it
     516                                        k = *itr_i;
     517                                        forbidden[ color[k] ] = (k != i);
     518                                }
     519                        }
     520                        k = 0;
     521                        while( forbidden[k] && k < m )
     522                        {       k++;
     523                                CPPAD_ASSERT_UNKNOWN( k < n );
     524                        }
     525                        color[i] = k;
     526                }
     527                size_t n_color = 1;
     528                for(k = 0; k < m; k++)
     529                        n_color = std::max(n_color, color[k] + 1);
     530
     531                // weight vector for calls to reverse
     532                VectorBase w(m);
     533
     534                // location for return values from Reverse
     535                VectorBase dw(n);
     536
     537                // loop over colors
     538                size_t c;
     539                for(c = 0; c < n_color; c++)
     540                {       // determine all the rows with this color
     541                        for(i = 0; i < m; i++)
     542                        {       if( color[i] == c )
     543                                        w[i] = one;
     544                                else    w[i] = zero;
     545                        }
     546                        // call reverse mode for all these rows at once
     547                        dw = Reverse(1, w);
     548
     549                        // set the corresponding components of the result
     550                        for(i = 0; i < m; i++) if( color[i] == c )
     551                        {       for(
     552                                        itr_j = p[i].begin();
     553                                        itr_j != p[i].end();
     554                                        itr_j++
     555                                )
     556                                {       j = *itr_j;
     557                                        jac[i * n + j] = dw[j];
     558                                }
     559                        }
     560                }
     561        }
     562}
     563
     564/*!
     565Compute a sparse Jacobian.
     566
     567The C++ source code corresponding to this operation is
     568\verbatim
     569        jac = SparseJacobian(x, p)
     570\endverbatim
     571
     572\tparam Base
     573is the base type for the recording that is stored in this
     574ADFun<Base object.
     575
     576\tparam VectorBase
     577is a simple vector class with elements of the \a Base.
     578
     579\tparam VectorSet
     580is a simple vector class with elements of type
     581\c bool or \c std::set<size_t>.
     582
     583\param x
     584is a vector specifing the point at which to compute the Jacobian.
     585
     586\param p
     587is the sparsity pattern for the Jacobian that we are calculating.
     588
     589\return
     590Will be a vector if size \c m * n containing the Jacobian at the
     591specified point (in row major order).
     592*/
     593template <class Base>
     594template <class VectorBase, class VectorSet>
     595VectorBase ADFun<Base>::SparseJacobian(
     596        const VectorBase& x, const VectorSet& p
     597)
     598{       size_t m = Range();
     599        size_t n = Domain();
     600        VectorBase jac(m * n);
     601
     602        typedef typename VectorSet::value_type Set_type;
     603        SparseJacobianCase( Set_type(), x, p, jac);
     604
     605        return jac;
     606}
     607
     608/*!
     609Compute a sparse Jacobian.
     610
     611The C++ source code corresponding to this operation is
     612\verbatim
     613        jac = SparseJacobian(x)
     614\endverbatim
     615
     616\tparam Base
     617is the base type for the recording that is stored in this
     618ADFun<Base object.
     619
     620\tparam VectorBase
     621is a simple vector class with elements of the \a Base.
     622
     623\param x
     624is a vector specifing the point at which to compute the Jacobian.
     625
     626\return
     627Will be a vector if size \c m * n containing the Jacobian at the
     628specified point (in row major order).
     629*/
     630template <class Base>
     631template <class VectorBase>
     632VectorBase ADFun<Base>::SparseJacobian( const VectorBase& x )
     633{       typedef CppAD::vector<bool>   VectorBool;
     634
     635        size_t m = Range();
     636        size_t n = Domain();
     637
     638        // sparsity pattern for Jacobian
     639        VectorBool p(m * n);
     640
     641        // return result
     642        VectorBase jac(m * n);
     643
     644        if( n <= m )
     645        {       size_t j, k;
     646
     647                // use forward mode
     648                VectorBool r(n * n);
     649                for(j = 0; j < n; j++)
     650                {       for(k = 0; k < n; k++)
     651                                r[j * n + k] = false;
     652                        r[j * n + j] = true;
     653                }
     654                p = ForSparseJac(n, r);
     655        }
     656        else
     657        {       size_t i, k;
     658
     659                // use reverse mode
     660                VectorBool s(m * m);
     661                for(i = 0; i < m; i++)
     662                {       for(k = 0; k < m; k++)
     663                                s[i * m + k] = false;
     664                        s[i * m + i] = true;
     665                }
     666                p = RevSparseJac(m, s);
     667        }
     668        bool set_type = true; // only type is used (to dispatch to proper case)
     669        SparseJacobianCase(set_type, x, p, jac);
    342670        return jac;
    343671}
    344672
    345 } // END CppAD namespace
    346 
     673
     674CPPAD_END_NAMESPACE
    347675# endif
  • trunk/example/sparse_jacobian.cpp

    r1531 r1551  
    3535# include <cppad/cppad.hpp>
    3636namespace { // ---------------------------------------------------------
    37 // define the template function in empty namespace
    38 template <class VectorBase, class VectorBool>
    39 bool reverse_case()
     37bool reverse()
    4038{       bool ok = true;
    4139        using CppAD::AD;
     
    6260
    6361        // new value for the independent variable vector
    64         VectorBase x(n);
     62        CPPAD_TEST_VECTOR<double> x(n);
    6563        for(j = 0; j < n; j++)
    6664                x[j] = double(j);
    6765
    6866        // Jacobian of y without sparsity pattern
    69         VectorBase jac(m * n);
     67        CPPAD_TEST_VECTOR<double> jac(m * n);
    7068        jac = f.SparseJacobian(x);
    7169        /*
     
    7472              [ 1 1 1 x_3]
    7573        */
    76         VectorBase check(m * n);
     74        CPPAD_TEST_VECTOR<double> check(m * n);
    7775        check[0] = 1.; check[1] = 1.; check[2]  = 0.; check[3]  = 0.;
    7876        check[4] = 0.; check[5] = 0.; check[6]  = 1.; check[7]  = 1.;
     
    8179                ok &=  NearEqual(check[k], jac[k], 1e-10, 1e-10 );
    8280
    83         // test passing sparsity pattern
    84         VectorBool s(m * m);
    85         VectorBool p(m * n);
     81        // using packed boolean sparsity patterns
     82        CppAD::vectorBool s_b(m * m), p_b(m * n);
    8683        for(i = 0; i < m; i++)
    8784        {       for(k = 0; k < m; k++)
    88                         s[i * m + k] = false;
    89                 s[i * m + i] = true;
     85                        s_b[i * m + k] = false;
     86                s_b[i * m + i] = true;
    9087        }
    91         p   = f.RevSparseJac(m, s);
    92         jac = f.SparseJacobian(x);
     88        p_b   = f.RevSparseJac(m, s_b);
     89        jac   = f.SparseJacobian(x, p_b);
     90        for(k = 0; k < 12; k++)
     91                ok &=  NearEqual(check[k], jac[k], 1e-10, 1e-10 );
     92
     93        // using vector of sets sparsity patterns
     94        std::vector< std::set<size_t> > s_s(m),  p_s(m);
     95        for(i = 0; i < m; i++)
     96                s_s[i].insert(i);
     97        p_s   = f.RevSparseJac(m, s_s);
     98        jac   = f.SparseJacobian(x, p_s);
    9399        for(k = 0; k < 12; k++)
    94100                ok &=  NearEqual(check[k], jac[k], 1e-10, 1e-10 );
     
    97103}
    98104
    99 template <class VectorBase, class VectorBool>
    100 bool forward_case()
     105bool forward()
    101106{       bool ok = true;
    102107        using CppAD::AD;
     
    124129
    125130        // new value for the independent variable vector
    126         VectorBase x(n);
     131        CPPAD_TEST_VECTOR<double> x(n);
    127132        for(j = 0; j < n; j++)
    128133                x[j] = double(j);
    129134
    130135        // Jacobian of y without sparsity pattern
    131         VectorBase jac(m * n);
     136        CPPAD_TEST_VECTOR<double> jac(m * n);
    132137        jac = f.SparseJacobian(x);
    133138        /*
     
    137142              [ 0 1 x_2 ]
    138143        */
    139         VectorBase check(m * n);
     144        CPPAD_TEST_VECTOR<double> check(m * n);
    140145        check[0] = 1.; check[1]  = 0.; check[2]  = 1.;
    141146        check[3] = 1.; check[4]  = 0.; check[5]  = 1.;
     
    145150                ok &=  NearEqual(check[k], jac[k], 1e-10, 1e-10 );
    146151
    147         // test passing sparsity pattern
    148         VectorBool r(n * n);
    149         VectorBool p(m * n);
     152        // test using packed boolean vectors for sparsity pattern
     153        CppAD::vectorBool r_b(n * n), p_b(m * n);
    150154        for(j = 0; j < n; j++)
    151155        {       for(k = 0; k < n; k++)
    152                         r[j * n + k] = false;
    153                 r[j * n + j] = true;
     156                        r_b[j * n + k] = false;
     157                r_b[j * n + j] = true;
    154158        }
    155         p   = f.ForSparseJac(n, r);
    156         jac = f.SparseJacobian(x);
     159        p_b = f.ForSparseJac(n, r_b);
     160        jac = f.SparseJacobian(x, p_b);
     161        for(k = 0; k < 12; k++)
     162                ok &=  NearEqual(check[k], jac[k], 1e-10, 1e-10 );
     163
     164        // test using vector of sets for sparsity pattern
     165        std::vector< std::set<size_t> > r_s(n), p_s(m);
     166        for(j = 0; j < n; j++)
     167                r_s[j].insert(j);
     168        p_s = f.ForSparseJac(n, r_s);
     169        jac = f.SparseJacobian(x, p_s);
    157170        for(k = 0; k < 12; k++)
    158171                ok &=  NearEqual(check[k], jac[k], 1e-10, 1e-10 );
     
    161174}
    162175} // End empty namespace
    163 # include <vector>
    164 # include <valarray>
     176
    165177bool sparse_jacobian(void)
    166178{       bool ok = true;
    167         // Run with VectorBase equal to three different cases
    168         // all of which are Simple Vectors with elements of type double.
    169         // Also vary the type of vector for VectorBool.
    170         ok &= forward_case< CppAD::vector<double>, CppAD::vectorBool   >();
    171         ok &= reverse_case< CppAD::vector<double>, CppAD::vector<bool> >();
    172         //
    173         ok &= forward_case< std::vector<double>,   std::vector<bool>   >();
    174         ok &= reverse_case< std::vector<double>,   std::valarray<bool> >();
    175         //
    176         ok &= forward_case< std::valarray<double>, CppAD::vectorBool   >();
    177         ok &= reverse_case< std::valarray<double>, CppAD::vector<bool> >();
    178         //
     179        ok &= forward();
     180        ok &= reverse();
     181
    179182        return ok;
    180183}
  • trunk/omh/whats_new_09.omh

    r1550 r1551  
    5959trying to read and understand the CppAD source code.)
    6060
     61$head 10-16$$
     62The $cref/sparse_jacobian/$$ routine was extended so the user can now choose
     63between vectors of sets and boolean vectors for representing
     64$cref/sparsity patterns/glossary/Sparsity Pattern/$$.
     65
    6166$head 10-14$$
    6267The $icode packed$$ parameter for the sparsity routines
  • trunk/test_more/jacobian.cpp

    r1381 r1551  
    5050{       bool ok = true;
    5151        using CppAD::vector;
    52         size_t j, k;
     52        size_t i, j, k;
    5353
    5454        size_t n = 4;
     
    6969        check = eval_jac_g(x);
    7070
     71        // regular jacobian
    7172        vector<double> jac_g = fun_g.Jacobian(x);
    72 
    7373        for(k = 0; k < m *n; k++)
    7474                ok &= CppAD::NearEqual(jac_g[k], check[k], 1e-10, 1e-10);
    7575
     76        // one argument sparse jacobian
    7677        jac_g = fun_g.SparseJacobian(x);
     78        for(k = 0; k < m *n; k++)
     79                ok &= CppAD::NearEqual(jac_g[k], check[k], 1e-10, 1e-10);
    7780
     81        // sparse jacobian using bool vectors
     82        CPPAD_TEST_VECTOR<bool> p_b(m * n) , r_b(n * n);
     83        for(i = 0; i < n; i++)
     84                for(j = 0; j < n; j++)
     85                        r_b[i * n + j] = (i == j);
     86        p_b = fun_g.ForSparseJac(n, r_b);
     87        jac_g = fun_g.SparseJacobian(x, p_b);
     88        for(k = 0; k < m *n; k++)
     89                ok &= CppAD::NearEqual(jac_g[k], check[k], 1e-10, 1e-10);
     90
     91        // sparse jacobian using vectors of sets
     92        std::vector< std::set<size_t> > p_s(m) , r_s(n);
     93        for(i = 0; i < n; i++)
     94                for(j = 0; j < n; j++)
     95                        r_s[i].insert(j);
     96        p_s = fun_g.ForSparseJac(n, r_s);
     97        jac_g = fun_g.SparseJacobian(x, p_s);
    7898        for(k = 0; k < m *n; k++)
    7999                ok &= CppAD::NearEqual(jac_g[k], check[k], 1e-10, 1e-10);
  • trunk/test_more/makefile.am

    r1497 r1551  
    119119        sin_cos.cpp \
    120120        sinh.cpp \
     121        sparse_jacobian.cpp \
    121122        sparse_vec_ad.cpp \
    122123        sqrt.cpp \
  • trunk/test_more/test_more.cpp

    r1497 r1551  
    7171extern bool SinCos(void);
    7272extern bool Sinh(void);
     73extern bool sparse_jacobian(void);
    7374extern bool sparse_vec_ad(void);
    7475extern bool Sqrt(void);
     
    165166        ok &= Run( SinCos,          "SinCos"         );
    166167        ok &= Run( Sinh,            "Sinh"           );
     168        ok &= Run( sparse_jacobian, "sparse_jacobian");
    167169        ok &= Run( sparse_vec_ad,   "sparse_vec_ad"  );
    168170        ok &= Run( Sqrt,            "Sqrt"           );
Note: See TracChangeset for help on using the changeset viewer.