Changeset 1558


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

trunk: Extend sparse Hessian driver to optionally use set representation of sparsity.

check_doxygen.sh: add sparse_hessian.hpp to list of documented files.
*/sparse_hessian.cpp: example and test of using sets for sparse calculations.
makefile.am: add sparse_hessian to list of tests.
makefile.in: changes automatically generated by change in makefile.am.
test_more.cpp: add sparse_hessian to list of tests.
sparse_jacobina.cpp: minor edits.
whats_new_09.omh: user's view of the changes.
ad_fun.hpp: add new sparse Hessian private helper functions to class.
sparse_jacobian.hpp: minor edits.
sparse_hessian.cpp: change to allow for vectors of sets sparsity patterns.

Location:
trunk
Files:
1 added
9 edited

Legend:

Unmodified
Added
Removed
  • trunk/check_doxygen.sh

    r1551 r1558  
    6767        sub_op.hpp
    6868        sparse_jacobian.hpp
     69        sparse_hessian.hpp
    6970        sparse_pack.hpp
    7071        sparse_set.hpp
  • trunk/cppad/local/ad_fun.hpp

    r1551 r1558  
    199199                VectorBase&              jac
    200200        );
     201        // ------------------------------------------------------------
     202        // vector of bool version of SparseHessian
     203        // (see doxygen in sparse_hessian.hpp)
     204        template <class VectorBase, class VectorSet>
     205        void SparseHessianCase(
     206                bool                     set_type    ,
     207                const VectorBase&        x           ,
     208                const VectorBase&        w           ,
     209                const VectorSet&         p           ,
     210                VectorBase&              hes
     211        );
     212        // vector of std::set<size_t> version of SparseHessian
     213        // (see doxygen in sparse_hessian.hpp)
     214        template <class VectorBase, class VectorSet>
     215        void SparseHessianCase(
     216                const std::set<size_t>&  set_type    ,
     217                const VectorBase&        x           ,
     218                const VectorBase&        w           ,
     219                const VectorSet&         p           ,
     220                VectorBase&              hes
     221        );
    201222// ------------------------------------------------------------
    202223public:
  • trunk/cppad/local/sparse_hessian.hpp

    r1531 r1558  
    1717$begin sparse_hessian$$
    1818$spell
     19        CppAD
     20        valarray
     21        std
    1922        Bool
    2023        hes
     
    3336
    3437$head Purpose$$
    35 We use $latex F : B^n \rightarrow B^m$$ do denote the
     38We use $latex F : \R^n \rightarrow \R^m$$ do denote the
    3639$cref/AD function/glossary/AD Function/$$
    3740corresponding to $icode f$$.
     
    4043        hes = \dpow{2}{x} \sum_{i=1}^m w_i F_i (x)
    4144\] $$
    42 This is a preliminary implementation of a method for using the fact
    43 that the matrix is sparse to reduce the amount of computation necessary.
     45This routine assumes  that the matrix $icode hes$$ is sparse
     46and uses this assumption sparse to reduce the amount of computation necessary.
    4447One should use speed tests to verify that results are computed faster
    4548than when using the routine $cref/Hessian/$$.
     
    5659The argument $icode x$$ has prototype
    5760$codei%
    58         const %VectorBase% &%x%
     61        const %VectorBase%& %x%
    5962%$$
    6063(see $cref/VectorBase/sparse_hessian/VectorBase/$$ below)
     
    6871The argument $icode w$$ has prototype
    6972$codei%
    70         const %VectorBase% &%w%
     73        const %VectorBase%& %w%
    7174%$$
    7275and size $latex m$$.
     
    7477for $icode hes$$.
    7578The more components of $latex w$$ that are identically zero,
    76 the more spares the resulting Hessian may be (and hence the more efficient
     79the more sparse the resulting Hessian may be (and hence the more efficient
    7780the calculation of $icode hes$$ may be).
    7881
    7982$head p$$
    8083The argument $icode p$$ is optional and has prototype
    81 $syntax%
    82         const %VectorBool% &%p%
    83 %$$
    84 (see $cref/VectorBool/sparse_hessian/VectorBool/$$ below)
    85 and its size is $latex n * n$$.
    86 It specifies a
     84$codei%
     85        const %VectorSet%& %p%
     86%$$
     87(see $cref/VectorSet/sparse_hessian/VectorSet/$$ below)
     88If it has elements of type $code bool$$,
     89its size is $latex n * n$$.
     90If it has elements of type $code std::set<size_t>$$,
     91its size is $latex n$$ and all its set elements are between
     92zero and $latex n - 1$$.
     93It specifies a
    8794$cref/sparsity pattern/glossary/Sparsity Pattern/$$
    88 for the Hessian; i.e.,
    89 for $latex j = 0 , \ldots , n-1$$ and $latex k = 0 , \ldots , n-1$$.
    90 $latex \[
    91         \sum_i w_i \DD{ F_i }{ x_j }{ x_k } \neq 0
    92                 ; \Rightarrow \; p [ j * n + k ] = {\rm true}
    93 \] $$
     95for the Hessian $icode hes$$.
    9496$pre
    9597
     
    98100$codei SparseHessian$$, it should be faster to calculate $icode p$$ once and
    99101pass this argument to $codei SparseHessian$$.
     102In addition,
     103if you specify $icode p$$, CppAD will use the same
     104type of sparsity representation
     105(vectors of $code bool$$ or vectors of $code std::set<size_t>$$)
     106for its internal calculations.
     107Otherwise, the representation
     108for the internal calculations is unspecified.
    100109
    101110$head hes$$
     
    118127if this is not the case.
    119128
    120 $head VectorBool$$
    121 The type $icode VectorBool$$ must be a $xref/SimpleVector/$$ class with
    122 $xref/SimpleVector/Elements of Specified Type/elements of type bool/$$.
     129$head VectorSet$$
     130The type $icode VectorSet$$ must be a $xref/SimpleVector/$$ class with
     131$xref/SimpleVector/Elements of Specified Type/elements of type/$$
     132$code bool$$ or $code std::set<size_t>$$;
     133see $cref/sparsity pattern/glossary/Sparsity Pattern/$$ for a discussion
     134of the difference.
    123135The routine $xref/CheckSimpleVector/$$ will generate an error message
    124136if this is not the case.
    125 In order to save memory,
    126 you may want to use a class that packs multiple elements into one
    127 storage location; for example,
    128 $cref/vectorBool/CppAD_vector/vectorBool/$$.
     137
     138$subhead Restrictions$$
     139If $icode VectorSet$$ has elements of $code std::set<size_t>$$,
     140then $icode%p%[%i%]%$$ must return a reference (not a copy) to the
     141corresponding set.
     142According to section 26.3.2.3 of the 1998 C++ standard,
     143$code std::valarray< std::set<size_t> >$$ does not satisfy
     144this condition.
    129145
    130146$head Uses Forward$$
     
    148164-----------------------------------------------------------------------------
    149165*/
    150 
    151 
    152 
    153 
    154 //  BEGIN CppAD namespace
    155 namespace CppAD {
    156 
    157 template <typename Base>
    158 template <typename VectorBase>
     166CPPAD_BEGIN_NAMESPACE
     167
     168#
     169/*!
     170\file sparse_hessian.hpp
     171Sparse Hessian driver routine and helper functions.
     172*/
     173
     174/*!
     175Private helper function for SparseHessian(x, w, p).
     176
     177All of the description in the public member function SparseHessian(x, w, p)
     178applies.
     179
     180\param set_type
     181is a \c bool value. This argument is used to dispatch to the proper source
     182code depending on the value of \c VectorSet::value_type.
     183
     184\param x
     185See \c SparseHessian(x, w, p).
     186
     187\param w
     188See \c SparseHessian(x, w, p).
     189
     190\param p
     191See \c SparseHessian(x, w, p).
     192
     193\param hes
     194is the return value for the corresponding call to \c SparseHessian(x, w, p).
     195On input, it must have size equal to the domain times range dimension
     196for this ADFun<Base> object.
     197On return, it will contain the Hessian.
     198*/
     199
     200template <class Base>
     201template <class VectorBase, class VectorSet>
     202void ADFun<Base>::SparseHessianCase(
     203        bool               set_type         ,
     204        const VectorBase&  x                ,
     205        const VectorBase&  w                ,
     206        const VectorSet&   p                ,
     207        VectorBase&        hes              )
     208{
     209        typedef CppAD::vector<size_t> SizeVector;
     210        typedef CppAD::vectorBool     VectorBool;
     211        size_t j, k, l;
     212
     213        size_t n = Domain();
     214
     215        // check Vector is Simple Vector class with bool elements
     216        CheckSimpleVector<bool, VectorSet>();
     217
     218        // check Vector is Simple Vector class with Base type elements
     219        CheckSimpleVector<Base, VectorBase>();
     220
     221        CPPAD_ASSERT_KNOWN(
     222                x.size() == n,
     223                "SparseHessian: size of x not equal domain dimension for f"
     224        );
     225        CPPAD_ASSERT_KNOWN(
     226                w.size() == Range(),
     227                "SparseHessian: size of w not equal range dimension for f"
     228        );
     229        CPPAD_ASSERT_KNOWN(
     230                p.size() == n * n,
     231                "SparseHessian: using bool values and size of p "
     232                "not equal square of domain dimension for f"
     233        );
     234       
     235        // initial coloring
     236        SizeVector color(n);
     237        for(j = 0; j < n; j++)
     238                color[j] = j;
     239
     240        // See GreedyPartialD2Coloring Algorithm Section 3.6.2 of
     241        // Graph Coloring in Optimization Revisited by
     242        // Assefaw Gebremedhin, Fredrik Maane, Alex Pothen
     243        VectorBool forbidden(n);
     244        for(j = 0; j < n; j++)
     245        {       // initial all colors as ok for this column
     246                for(k = 0; k < n; k++)
     247                        forbidden[k] = false;
     248                // for each row that is connected to column j
     249                for(l = 0; l < n; l++) if( p[l * n + j] )
     250                {       // for each column that is connected to row l
     251                        for(k = 0; k < n; k++) if( p[l * n + k] & (j != k) )   
     252                                forbidden[ color[k] ] = true;
     253                }
     254                k = 0;
     255                while( forbidden[k] && k < n )
     256                {       k++;
     257                        CPPAD_ASSERT_UNKNOWN( k < n );
     258                }
     259                color[j] = k;
     260        }
     261        size_t n_color = 1;
     262        for(k = 0; k < n; k++) n_color = std::max(n_color, color[k] + 1);
     263
     264        // some values
     265        const Base zero(0);
     266        const Base one(1);
     267
     268        // point at which we are evaluating the Hessian
     269        Forward(0, x);
     270
     271        // initialize the return value
     272        for(j = 0; j < n; j++)
     273                for(k = 0; k < n; k++)
     274                        hes[j * n + k] = zero;
     275
     276        // direction vector for calls to forward
     277        VectorBase u(n);
     278
     279        // location for return values from Reverse
     280        VectorBase ddw(n * 2);
     281
     282        // loop over colors
     283        size_t c;
     284        for(c = 0; c < n_color; c++)
     285        {       // determine all the colums with this color
     286                for(j = 0; j < n; j++)
     287                {       if( color[j] == c )
     288                                u[j] = one;
     289                        else    u[j] = zero;
     290                }
     291                // call forward mode for all these columns at once
     292                Forward(1, u);
     293
     294                // evaluate derivative of w^T * F'(x) * u
     295                ddw = Reverse(2, w);
     296
     297                // set the corresponding components of the result
     298                for(j = 0; j < n; j++) if( color[j] == c )
     299                {       for(l = 0; l < n; l++)
     300                                if( p[ l * n + j ] )
     301                                        hes[l * n + j] = ddw[l * 2 + 1];
     302                }
     303        }
     304}
     305/*!
     306Private helper function for SparseHessian(x, w, p).
     307
     308All of the description in the public member function SparseHessian(x, w, p)
     309applies.
     310
     311\param set_type
     312is a \c std::set<size_t> value.
     313This argument is used to dispatch to the proper source
     314code depending on the value of \c VectorSet::value_type.
     315
     316\param x
     317See \c SparseHessian(x, w, p).
     318
     319\param w
     320See \c SparseHessian(x, w, p).
     321
     322\param p
     323See \c SparseHessian(x, w, p).
     324
     325\param hes
     326is the return value for the corresponding call to \c SparseHessian(x, w, p).
     327On input, it must have size equal to the domain times range dimension
     328for this ADFun<Base> object.
     329On return, it will contain the Hessian.
     330*/
     331
     332template <class Base>
     333template <class VectorBase, class VectorSet>
     334void ADFun<Base>::SparseHessianCase(
     335        const std::set<size_t>&  set_type         ,
     336        const VectorBase&        x                ,
     337        const VectorBase&        w                ,
     338        const VectorSet&         p                ,
     339        VectorBase&              hes              )
     340{
     341        typedef CppAD::vector<size_t> SizeVector;
     342        typedef CppAD::vectorBool     VectorBool;
     343        size_t j, k, l;
     344
     345        size_t n = Domain();
     346
     347        // check VectorSet is Simple Vector class with sets for elements
     348        static std::set<size_t> two, three;
     349        if( two.empty() )
     350        {       two.insert(2);
     351                three.insert(3);
     352        }
     353        CPPAD_ASSERT_UNKNOWN( two.size() == 1 );
     354        CPPAD_ASSERT_UNKNOWN( three.size() == 1 );
     355        CheckSimpleVector<std::set<size_t>, VectorSet>(two, three);
     356
     357        // check Vector is Simple Vector class with Base type elements
     358        CheckSimpleVector<Base, VectorBase>();
     359
     360        CPPAD_ASSERT_KNOWN(
     361                x.size() == n,
     362                "SparseHessian: size of x not equal domain dimension for f"
     363        );
     364        CPPAD_ASSERT_KNOWN(
     365                w.size() == Range(),
     366                "SparseHessian: size of w not equal range dimension for f"
     367        );
     368        CPPAD_ASSERT_KNOWN(
     369                p.size() == n,
     370                "SparseHessian: using set values and size of p "
     371                "not equal domain dimension for f"
     372        );
     373       
     374        // initial coloring
     375        SizeVector color(n);
     376        for(j = 0; j < n; j++)
     377                color[j] = j;
     378
     379        // See GreedyPartialD2Coloring Algorithm Section 3.6.2 of
     380        // Graph Coloring in Optimization Revisited by
     381        // Assefaw Gebremedhin, Fredrik Maane, Alex Pothen
     382        VectorBool forbidden(n);
     383        std::set<size_t>::const_iterator itr_k, itr_l;
     384        for(j = 0; j < n; j++)
     385        {       // initial all colors as ok for this column
     386                for(k = 0; k < n; k++)
     387                        forbidden[k] = false;
     388
     389                // for each row that is connected to column j
     390                itr_l = p[j].begin();
     391                while( itr_l != p[j].end() )
     392                {       l = *itr_l++;
     393                        // for each column that is connected to row l
     394                        // (same as each row that is connect to column l)
     395                        itr_k = p[l].begin();
     396                        while( itr_k != p[l].end() )
     397                        {       k = *itr_k++;
     398                                forbidden[ color[k] ] = (j != k);
     399                        }
     400                }
     401                k = 0;
     402                while( forbidden[k] && k < n )
     403                {       k++;
     404                        CPPAD_ASSERT_UNKNOWN( k < n );
     405                }
     406                color[j] = k;
     407        }
     408        size_t n_color = 1;
     409        for(k = 0; k < n; k++) n_color = std::max(n_color, color[k] + 1);
     410
     411        // some values
     412        const Base zero(0);
     413        const Base one(1);
     414
     415        // point at which we are evaluating the Hessian
     416        Forward(0, x);
     417
     418        // initialize the return value
     419        for(j = 0; j < n; j++)
     420                for(k = 0; k < n; k++)
     421                        hes[j * n + k] = zero;
     422
     423        // direction vector for calls to forward
     424        VectorBase u(n);
     425
     426        // location for return values from Reverse
     427        VectorBase ddw(n * 2);
     428
     429        // loop over colors
     430        size_t c;
     431        for(c = 0; c < n_color; c++)
     432        {       // determine all the colums with this color
     433                for(j = 0; j < n; j++)
     434                {       if( color[j] == c )
     435                                u[j] = one;
     436                        else    u[j] = zero;
     437                }
     438                // call forward mode for all these columns at once
     439                Forward(1, u);
     440
     441                // evaluate derivative of w^T * F'(x) * u
     442                ddw = Reverse(2, w);
     443
     444                // set the corresponding components of the result
     445                for(j = 0; j < n; j++) if( color[j] == c )
     446                {       itr_l = p[j].begin();
     447                        while( itr_l != p[j].end() )
     448                        {       l = *itr_l++;
     449                                hes[l * n + j] = ddw[l * 2 + 1];
     450                        }
     451                }
     452        }
     453}
     454
     455/*!
     456Compute a sparse Hessian.
     457
     458The C++ source code coresponding to this operation is
     459\verbatim
     460        hes = SparseHessian(x, w)
     461\endverbatim
     462
     463
     464\tparam Base
     465is the base type for the recording that is stored in this
     466ADFun<Base object.
     467
     468\tparam VectorBase
     469is a simple vector class with elements of the \a Base.
     470
     471\tparam VectorSet
     472is a simple vector class with elements of type
     473\c bool or \c std::set<size_t>.
     474
     475\param x
     476is a vector specifing the point at which to compute the Hessian.
     477
     478\param w
     479The Hessian is computed for a weighted sum of the components
     480of the function corresponding to this ADFun<Base> object.
     481The argument \a w specifies the weights for each component.
     482It must have size equal to the range dimension for this ADFun<Base> object.
     483
     484\param p
     485is a sparsity pattern for the Hessian.
     486
     487\return
     488Will be a vector of size \c n * n containing the Hessian of
     489at the point specified by \a x
     490(where \c n is the domain dimension for this ADFun<Base> object).
     491*/
     492template <class Base>
     493template <class VectorBase, class VectorSet>
     494VectorBase ADFun<Base>::SparseHessian(
     495        const VectorBase& x, const VectorBase& w, const VectorSet& p
     496)
     497{       size_t n = Domain();
     498        VectorBase hes(n * n);
     499
     500        typedef typename VectorSet::value_type Set_type;
     501        SparseHessianCase(Set_type(), x, w, p, hes);
     502
     503        return hes;
     504}
     505
     506/*!
     507Compute a sparse Hessian
     508
     509The C++ source code coresponding to this operation is
     510\verbatim
     511        hes = SparseHessian(x, w)
     512\endverbatim
     513
     514
     515\tparam Base
     516is the base type for the recording that is stored in this
     517ADFun<Base object.
     518
     519\tparam VectorBase
     520is a simple vector class with elements of the \a Base.
     521
     522\param x
     523is a vector specifing the point at which to compute the Hessian.
     524
     525\param w
     526The Hessian is computed for a weighted sum of the components
     527of the function corresponding to this ADFun<Base> object.
     528The argument \a w specifies the weights for each component.
     529It must have size equal to the range dimension for this ADFun<Base> object.
     530
     531\return
     532Will be a vector of size \c n * n containing the Hessian of
     533at the point specified by \a x
     534(where \c n is the domain dimension for this ADFun<Base> object).
     535*/
     536template <class Base>
     537template <class VectorBase>
    159538VectorBase ADFun<Base>::SparseHessian(const VectorBase &x, const VectorBase &w)
    160539{       size_t i, j, k;
    161 
    162         typedef CppAD::vector<bool>   VectorBool;
     540        typedef CppAD::vectorBool VectorBool;
    163541
    164542        size_t m = Range();
     
    180558
    181559        // compute sparse Hessian
    182         return SparseHessian(x, w, p);
     560        VectorBase hes(n * n);
     561        bool set_type = true; // only used to dispatch compiler to proper case
     562        SparseHessianCase(set_type, x, w, p, hes);
     563
     564        return hes;
    183565}
    184566
    185 template <typename Base>
    186 template <typename VectorBase, typename VectorBool>
    187 VectorBase ADFun<Base>::SparseHessian(
    188         const VectorBase &x ,
    189         const VectorBase &w ,
    190         const VectorBool &p )
    191 {
    192         typedef CppAD::vector<size_t> SizeVector;
    193         size_t j, k, l;
    194 
    195         size_t n = Domain();
    196 
    197         // check Vector is Simple Vector class with bool elements
    198         CheckSimpleVector<bool, VectorBool>();
    199 
    200         // check Vector is Simple Vector class with Base type elements
    201         CheckSimpleVector<Base, VectorBase>();
    202 
    203         CPPAD_ASSERT_KNOWN(
    204                 x.size() == n,
    205                 "SparseHessian: size of x not equal domain dimension for f"
    206         );
    207         CPPAD_ASSERT_KNOWN(
    208                 w.size() == Range(),
    209                 "SparseHessian: size of w not equal range dimension for f"
    210         );
    211         CPPAD_ASSERT_KNOWN(
    212                 p.size() == n * n,
    213                 "SparseHessian: size of p not equal square of"
    214                 " domain dimension for f"
    215         );
    216        
    217         // initial coloring
    218         SizeVector color(n);
    219         for(j = 0; j < n; j++)
    220                 color[j] = j;
    221 
    222         // See GreedyPartialD2Coloring Algorithm Section 3.6.2 of
    223         // Graph Coloring in Optimization Revisited by
    224         // Assefaw Gebremedhin, Fredrik Maane, Alex Pothen
    225         VectorBool    forbidden(n);
    226         for(j = 0; j < n; j++)
    227         {       // initial all colors as ok for this column
    228                 for(k = 0; k < n; k++)
    229                         forbidden[k] = false;
    230                 // for each row that is connected to column j
    231                 for(l = 0; l < n; l++) if( p[l * n + j] )
    232                 {       // for each column that is connected to row l
    233                         for(k = 0; k < n; k++) if( p[l * n + k] & (j != k) )   
    234                                 forbidden[ color[k] ] = true;
    235                 }
    236                 k = 0;
    237                 while( forbidden[k] && k < n )
    238                 {       k++;
    239                         CPPAD_ASSERT_UNKNOWN( k < n );
    240                 }
    241                 color[j] = k;
    242         }
    243         size_t n_color = 1;
    244         for(k = 0; k < n; k++) n_color = std::max(n_color, color[k] + 1);
    245 
    246         // some values
    247         const Base zero(0);
    248         const Base one(1);
    249 
    250         // point at which we are evaluating the Hessian
    251         Forward(0, x);
    252 
    253         // define the return value
    254         VectorBase h(n * n);
    255         for(j = 0; j < n; j++)
    256                 for(k = 0; k < n; k++)
    257                         h[j * n + k] = zero;
    258 
    259         // direction vector for calls to forward
    260         VectorBase u(n);
    261 
    262         // location for return values from Reverse
    263         VectorBase ddw(n * 2);
    264 
    265         // loop over colors
    266         size_t c;
    267         for(c = 0; c < n_color; c++)
    268         {       // determine all the colums with this color
    269                 for(j = 0; j < n; j++)
    270                 {       if( color[j] == c )
    271                                 u[j] = one;
    272                         else    u[j] = zero;
    273                 }
    274                 // call forward mode for all these columns at once
    275                 Forward(1, u);
    276 
    277                 // evaluate derivative of w^T * F'(x) * u
    278                 ddw = Reverse(2, w);
    279 
    280                 // set the corresponding components of the result
    281                 for(j = 0; j < n; j++) if( color[j] == c )
    282                 {       for(l = 0; l < n; l++)
    283                                 if( p[ l * n + j ] )
    284                                         h[l * n + j] = ddw[l * 2 + 1];
    285                 }
    286         }
    287         return h;
    288 }
    289 
    290 } // END CppAD namespace
    291 
     567CPPAD_END_NAMESPACE
    292568# endif
  • trunk/cppad/local/sparse_jacobian.hpp

    r1556 r1558  
    1717$begin sparse_jacobian$$
    1818$spell
     19        valarray
    1920        std
    2021        CppAD
     
    4445\] $$
    4546This routine assumes
    46 that the matrix $latex F^{(1)} (x) \in \R^{m \times n}$$
    47 is uses this assumption to reduce the amount of computation necessary.
     47that the matrix $latex F^{(1)} (x) \in \R^{m \times n}$$ is sparse
     48and uses this assumption to reduce the amount of computation necessary.
    4849One should use speed tests (e.g. $cref/speed_test/$$)
    4950to verify that results are computed faster
     
    9697(vectors of $code bool$$ or vectors of $code std::set<size_t>$$)
    9798for its internal calculations.
    98 Otherwise the representation
     99Otherwise, the representation
    99100for the internal calculations is unspecified.
    100101
     
    127128if this is not the case.
    128129
     130$subhead Restrictions$$
     131If $icode VectorSet$$ has elements of $code std::set<size_t>$$,
     132then $icode%p%[%i%]%$$ must return a reference (not a copy) to the
     133corresponding set.
     134According to section 26.3.2.3 of the 1998 C++ standard,
     135$code std::valarray< std::set<size_t> >$$ does not satisfy
     136this condition.
     137
    129138$head Uses Forward$$
    130139After each call to $cref/Forward/$$,
     
    150159/*!
    151160\file sparse_jacobian.hpp
    152 Computing sparse Jacobians.
     161Sparse Jacobian driver routine and helper functions.
    153162*/
    154163
     
    161170\param set_type
    162171is a \c bool value. This argument is used to dispatch to the proper souce
    163 code depending on the vlaue of \c VectorSet::value_type.
     172code depending on the value of \c VectorSet::value_type.
    164173
    165174\param x
     
    171180\param jac
    172181is the return value for the corresponding call to \c SparseJacobian(x, p).
    173 On input it must have size equalt to the domain times range dimension
     182On input, it must have size equal to the domain times range dimension
    174183for this ADFun<Base> object.
    175184On return, it will contain the Jacobian.
     
    206215        CPPAD_ASSERT_KNOWN(
    207216                p.size() == m * n,
    208                 "SparseJacobian: size of p not equal range dimension times "
    209                 " domain dimension for f"
     217                "SparseJacobian: using bool values and size of p "
     218                " not equal range dimension times domain dimension for f"
    210219        );
    211220        CPPAD_ASSERT_UNKNOWN(jac.size() == m * n);
     
    402411        CPPAD_ASSERT_KNOWN(
    403412                p.size() == m,
    404                 "SparseJacobian: size of p not equal range dimension for f"
     413                "SparseJacobian: using sets and size of p "
     414                "not equal range dimension for f"
    405415        );
    406416        CPPAD_ASSERT_UNKNOWN(jac.size() == m * n);
     
    576586
    577587\tparam VectorBase
    578 is a simple vector class with elements of the \a Base.
     588is a simple vector class with elements of type \a Base.
    579589
    580590\tparam VectorSet
     
    626636
    627637\return
    628 Will be a vector if size \c m * n containing the Jacobian at the
     638Will be a vector of size \c m * n containing the Jacobian at the
    629639specified point (in row major order).
    630640*/
     
    632642template <class VectorBase>
    633643VectorBase ADFun<Base>::SparseJacobian( const VectorBase& x )
    634 {       typedef CppAD::vector<bool>   VectorBool;
     644{       typedef CppAD::vectorBool   VectorBool;
    635645
    636646        size_t m = Range();
     
    667677                p = RevSparseJac(m, s);
    668678        }
    669         bool set_type = true; // only type is used (to dispatch to proper case)
     679        bool set_type = true; // only used to dispatch compiler to proper case
    670680        SparseJacobianCase(set_type, x, p, jac);
    671681        return jac;
  • trunk/example/sparse_hessian.cpp

    r1531 r1558  
    3232*/
    3333// BEGIN PROGRAM
    34 
    3534# include <cppad/cppad.hpp>
    36 namespace { // ---------------------------------------------------------
    37 // define the template function in empty namespace
    38 template <class VectorBase, class VectorBool>
    39 bool Case()
     35bool sparse_hessian(void)
    4036{       bool ok = true;
    4137        using CppAD::AD;
     
    6056
    6157        // new value for the independent variable vector
    62         VectorBase x(n);
     58        CPPAD_TEST_VECTOR<double> x(n);
    6359        for(i = 0; i < n; i++)
    6460                x[i] = double(i);
    6561
    6662        // second derivative of y[1]
    67         VectorBase w(m);
     63        CPPAD_TEST_VECTOR<double> w(m);
    6864        w[0] = 1.;
    69         VectorBase h( n * n );
     65        CPPAD_TEST_VECTOR<double> h( n * n );
    7066        h = f.SparseHessian(x, w);
    7167        /*
     
    7470            [ 0 0 2 ]
    7571        */
    76         VectorBase check(n * n);
     72        CPPAD_TEST_VECTOR<double> check(n * n);
    7773        check[0] = 2.; check[1] = 1.; check[2] = 0.;
    7874        check[3] = 1.; check[4] = 2.; check[5] = 0.;
     
    8177                ok &=  NearEqual(check[k], h[k], 1e-10, 1e-10 );
    8278
     79        // use vectors of bools to compute sparse hessian -------------------
    8380        // determine the sparsity pattern p for Hessian of w^T F
    84         VectorBool r(n * n);
     81        CPPAD_TEST_VECTOR<bool> r_bool(n * n);
    8582        for(j = 0; j < n; j++)
    8683        {       for(k = 0; k < n; k++)
    87                         r[j * n + k] = false;
    88                 r[j * n + j] = true;
     84                        r_bool[j * n + k] = false;
     85                r_bool[j * n + j] = true;
    8986        }
    90         f.ForSparseJac(n, r);
     87        f.ForSparseJac(n, r_bool);
    9188        //
    92         VectorBool s(m);
     89        CPPAD_TEST_VECTOR<bool> s_bool(m);
    9390        for(i = 0; i < m; i++)
    94                 s[i] = w[i] != 0;
    95         VectorBool p = f.RevSparseHes(n, s);
     91                s_bool[i] = w[i] != 0;
     92        CPPAD_TEST_VECTOR<bool> p_bool = f.RevSparseHes(n, s_bool);
    9693
    9794        // test passing sparsity pattern
    98         h = f.SparseHessian(x, w, p);
     95        h = f.SparseHessian(x, w, p_bool);
     96        for(k = 0; k < n * n; k++)
     97                ok &=  NearEqual(check[k], h[k], 1e-10, 1e-10 );
     98
     99        // use vectors of sets to compute sparse hessian ------------------
     100        // determine the sparsity pattern p for Hessian of w^T F
     101        CPPAD_TEST_VECTOR< std::set<size_t> > r_set(n);
     102        for(j = 0; j < n; j++)
     103                r_set[j].insert(j);
     104        f.ForSparseJac(n, r_set);
     105        //
     106        CPPAD_TEST_VECTOR< std::set<size_t> > s_set(1);
     107        for(i = 0; i < m; i++)
     108                if( w[i] != 0. )
     109                        s_set[0].insert(i);
     110        CPPAD_TEST_VECTOR< std::set<size_t> > p_set = f.RevSparseHes(n, s_set);
     111
     112        // test passing sparsity pattern
     113        h = f.SparseHessian(x, w, p_set);
    99114        for(k = 0; k < n * n; k++)
    100115                ok &=  NearEqual(check[k], h[k], 1e-10, 1e-10 );
     
    102117        return ok;
    103118}
    104 } // End empty namespace
    105 # include <vector>
    106 # include <valarray>
    107 bool sparse_hessian(void)
    108 {       bool ok = true;
    109         // Run with VectorBase equal to three different cases
    110         // all of which are Simple Vectors with elements of type double.
    111         // Also vary the type of vector for VectorBool.
    112         ok &= Case< CppAD::vector  <double>, CppAD::vectorBool   >();
    113         ok &= Case< std::vector    <double>, CppAD::vector<bool> >();
    114         ok &= Case< std::valarray  <double>, std::vector<bool>   >();
    115         return ok;
    116 }
    117119// END PROGRAM
  • trunk/omh/whats_new_09.omh

    r1556 r1558  
    5959trying to read and understand the CppAD source code.)
    6060
     61$head 10-23$$
     62The $cref/sparse_hessian/$$ 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-21$$
    6267The $cref/CheckSimpleVector/$$ function was extended so that
  • trunk/test_more/makefile.am

    r1556 r1558  
    120120        sin_cos.cpp \
    121121        sinh.cpp \
     122        sparse_hessian.cpp \
    122123        sparse_jacobian.cpp \
    123124        sparse_vec_ad.cpp \
  • trunk/test_more/sparse_jacobian.cpp

    r1554 r1558  
    264264bool sparse_jacobian(void)
    265265{       bool ok = true;
    266         typedef std::vector< std::set<size_t> >   std_vector_set;
    267         typedef std::valarray< std::set<size_t> > std_valarray_set;
    268         typedef CppAD::vector< std::set<size_t> > cppad_vector_set;
    269266        // ---------------------------------------------------------------
    270267        // vector of bool cases
     
    279276        // ---------------------------------------------------------------
    280277        // vector of set cases
     278        typedef std::vector< std::set<size_t> >   std_vector_set;
     279        typedef std::valarray< std::set<size_t> > std_valarray_set;
     280        typedef CppAD::vector< std::set<size_t> > cppad_vector_set;
     281        //
    281282        ok &= forward_set< CppAD::vector<double>, std_vector_set   >();
    282283        ok &= reverse_set< std::valarray<double>, std_vector_set   >();
     
    285286        ok &= reverse_set< CppAD::vector<double>, cppad_vector_set >();
    286287        //
    287 # ifndef _MSC_VER
    288         // there appears to be a bug in valarray on Microsoft compilers
    289         ok &= forward_set< std::valarray<double>, std_valarray_set >();
    290         ok &= reverse_set< std::valarray<double>, std_valarray_set >();
    291 # endif
     288        // According to section 26.3.2.3 of the 1998 C++ standard
     289        // a const valarray does not return references to its elements.
     290        // ok &= forward_set< std::valarray<double>, std_valarray_set >();
     291        // ok &= reverse_set< std::valarray<double>, std_valarray_set >();
    292292        // ---------------------------------------------------------------
    293293        return ok;
  • trunk/test_more/test_more.cpp

    r1556 r1558  
    7272extern bool SinCos(void);
    7373extern bool Sinh(void);
     74extern bool sparse_hessian(void);
    7475extern bool sparse_jacobian(void);
    7576extern bool sparse_vec_ad(void);
     
    168169        ok &= Run( SinCos,          "SinCos"         );
    169170        ok &= Run( Sinh,            "Sinh"           );
     171        ok &= Run( sparse_hessian,  "sparse_hessian" );
    170172        ok &= Run( sparse_jacobian, "sparse_jacobian");
    171173        ok &= Run( sparse_vec_ad,   "sparse_vec_ad"  );
Note: See TracChangeset for help on using the changeset viewer.