CppAD: A C++ Algorithmic Differentiation Package 20110419
sparse_jacobian.hpp
Go to the documentation of this file.
00001 /* $Id$ */
00002 # ifndef CPPAD_SPARSE_JACOBIAN_INCLUDED
00003 # define CPPAD_SPARSE_JACOBIAN_INCLUDED
00004 
00005 /* --------------------------------------------------------------------------
00006 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-10 Bradley M. Bell
00007 
00008 CppAD is distributed under multiple licenses. This distribution is under
00009 the terms of the 
00010                     Common Public License Version 1.0.
00011 
00012 A copy of this license is included in the COPYING file of this distribution.
00013 Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
00014 -------------------------------------------------------------------------- */
00015 
00016 /*
00017 $begin sparse_jacobian$$
00018 $spell
00019         valarray
00020         std
00021         CppAD
00022         Bool
00023         jac
00024         Jacobian
00025         const
00026         Taylor
00027 $$
00028 
00029 $section Sparse Jacobian: Easy Driver$$
00030 $index SparseJacobian$$
00031 $index jacobian, sparse$$
00032 
00033 $head Syntax$$
00034 $codei%%jac% = %f%.SparseJacobian(%x%)
00035 %$$
00036 $codei%%jac% = %f%.SparseJacobian(%x%, %p%)%$$
00037 
00038 $head Purpose$$
00039 We use $latex F : \R^n \rightarrow \R^m$$ do denote the
00040 $cref/AD function/glossary/AD Function/$$
00041 corresponding to $icode f$$. 
00042 The syntax above sets $icode jac$$ to the Jacobian 
00043 $latex \[
00044         jac = F^{(1)} (x) 
00045 \] $$
00046 This routine assumes
00047 that the matrix $latex F^{(1)} (x) \in \R^{m \times n}$$ is sparse
00048 and uses this assumption to reduce the amount of computation necessary.
00049 One should use speed tests (e.g. $cref/speed_test/$$)
00050 to verify that results are computed faster
00051 than when using the routine $cref/Jacobian/$$.
00052 
00053 $head f$$
00054 The object $icode f$$ has prototype
00055 $codei%
00056         ADFun<%Base%> %f%
00057 %$$
00058 Note that the $cref/ADFun/$$ object $icode f$$ is not $code const$$
00059 (see $cref/Uses Forward/sparse_jacobian/Uses Forward/$$ below).
00060 
00061 $head x$$
00062 The argument $icode x$$ has prototype
00063 $codei%
00064         const %VectorBase%& %x%
00065 %$$
00066 (see $cref/VectorBase/sparse_jacobian/VectorBase/$$ below)
00067 and its size 
00068 must be equal to $icode n$$, the dimension of the
00069 $cref/domain/seq_property/Domain/$$ space for $icode f$$.
00070 It specifies
00071 that point at which to evaluate the Jacobian.
00072 
00073 $head p$$
00074 The argument $icode p$$ is optional and has prototype
00075 $syntax%
00076         const %VectorSet%& %p%
00077 %$$
00078 (see $cref/VectorSet/sparse_jacobian/VectorSet/$$ below).
00079 If it has elements of type $code bool$$,
00080 its size is $latex m * n$$.
00081 If it has elements of type $code std::set<size_t>$$,
00082 its size is $latex m$$ and all its set elements are between
00083 zero and $latex n - 1$$.
00084 It specifies a 
00085 $cref/sparsity pattern/glossary/Sparsity Pattern/$$ 
00086 for the Jacobian $latex F^{(1)} (x)$$.
00087 $pre
00088 
00089 $$
00090 If this sparsity pattern does not change between calls to 
00091 $codei SparseJacobian$$, it should be faster to calculate $icode p$$ once 
00092 (using $cref/ForSparseJac/$$ or $cref/RevSparseJac/$$)
00093 and then pass $icode p$$ to $codei SparseJacobian$$.
00094 In addition,
00095 if you specify $icode p$$, CppAD will use the same
00096 type of sparsity representation 
00097 (vectors of $code bool$$ or vectors of $code std::set<size_t>$$)
00098 for its internal calculations.
00099 Otherwise, the representation
00100 for the internal calculations is unspecified.
00101 
00102 $head jac$$
00103 The result $icode jac$$ has prototype
00104 $codei%
00105         %VectorBase% %jac%
00106 %$$
00107 and its size is $latex m * n$$.
00108 For $latex i = 0 , \ldots , m - 1$$,
00109 and $latex j = 0 , \ldots , n - 1 $$ 
00110 $latex \[
00111         jac [ i * n + j ] = \D{ F_i }{ x_j } (x)
00112 \] $$
00113 
00114 $head VectorBase$$
00115 The type $icode VectorBase$$ must be a $cref/SimpleVector/$$ class with
00116 $cref/elements of type/SimpleVector/Elements of Specified Type/$$
00117 $icode Base$$.
00118 The routine $cref/CheckSimpleVector/$$ will generate an error message
00119 if this is not the case.
00120 
00121 $head VectorSet$$
00122 The type $icode VectorSet$$ must be a $xref/SimpleVector/$$ class with
00123 $xref/SimpleVector/Elements of Specified Type/elements of type/$$
00124 $code bool$$ or $code std::set<size_t>$$;
00125 see $cref/sparsity pattern/glossary/Sparsity Pattern/$$ for a discussion
00126 of the difference.
00127 The routine $cref/CheckSimpleVector/$$ will generate an error message
00128 if this is not the case.
00129 
00130 $subhead Restrictions$$
00131 If $icode VectorSet$$ has elements of $code std::set<size_t>$$,
00132 then $icode%p%[%i%]%$$ must return a reference (not a copy) to the 
00133 corresponding set.
00134 According to section 26.3.2.3 of the 1998 C++ standard,
00135 $code std::valarray< std::set<size_t> >$$ does not satisfy
00136 this condition. 
00137 
00138 $head Uses Forward$$
00139 After each call to $cref/Forward/$$,
00140 the object $icode f$$ contains the corresponding 
00141 $cref/Taylor coefficients/glossary/Taylor Coefficient/$$.
00142 After $code SparseJacobian$$,
00143 the previous calls to $xref/Forward/$$ are undefined.
00144 
00145 $head Example$$
00146 $children%
00147         example/sparse_jacobian.cpp
00148 %$$
00149 The routine
00150 $cref/sparse_jacobian.cpp/$$
00151 is examples and tests of $code sparse_jacobian$$.
00152 It return $code true$$, if it succeeds and $code false$$ otherwise.
00153 
00154 $end
00155 -----------------------------------------------------------------------------
00156 */
00157 
00158 CPPAD_BEGIN_NAMESPACE
00159 /*!
00160 \file sparse_jacobian.hpp
00161 Sparse Jacobian driver routine and helper functions.
00162 */
00163 
00164 /*!
00165 Private helper function for SparseJacobian(x, p).
00166 
00167 All of the description in the public member function SparseJacobian(x, p)
00168 applies.
00169 
00170 \param set_type
00171 is a \c bool value. This argument is used to dispatch to the proper souce
00172 code depending on the value of \c VectorSet::value_type.
00173 
00174 \param x
00175 See \c SparseJacobian(x, p).
00176 
00177 \param p
00178 See \c SparseJacobian(x, p).
00179 
00180 \param jac
00181 is the return value for the corresponding call to \c SparseJacobian(x, p).
00182 On input, it must have size equal to the domain times range dimension
00183 for this ADFun<Base> object.
00184 On return, it will contain the Jacobian.
00185 */
00186 template <class Base>
00187 template <class VectorBase, class VectorSet>
00188 void ADFun<Base>::SparseJacobianCase(
00189         bool               set_type        ,
00190         const VectorBase&  x               , 
00191         const VectorSet&   p               ,
00192         VectorBase&        jac             )
00193 {
00194         typedef CppAD::vector<size_t> SizeVector;
00195         typedef CppAD::vectorBool     VectorBool;
00196         size_t i, j, k;
00197 
00198         size_t m = Range();
00199         size_t n = Domain();
00200 
00201         // some values
00202         const Base zero(0);
00203         const Base one(1);
00204 
00205         // check VectorSet is Simple Vector class with bool elements
00206         CheckSimpleVector<bool, VectorSet>();
00207 
00208         // check VectorBase is Simple Vector class with Base type elements
00209         CheckSimpleVector<Base, VectorBase>();
00210 
00211         CPPAD_ASSERT_KNOWN(
00212                 x.size() == n,
00213                 "SparseJacobian: size of x not equal domain dimension for f"
00214         ); 
00215         CPPAD_ASSERT_KNOWN(
00216                 p.size() == m * n,
00217                 "SparseJacobian: using bool values and size of p "
00218                 " not equal range dimension times domain dimension for f"
00219         ); 
00220         CPPAD_ASSERT_UNKNOWN(jac.size() == m * n); 
00221 
00222         // point at which we are evaluating the Jacobian
00223         Forward(0, x);
00224 
00225         // initialize the return value
00226         for(i = 0; i < m; i++)
00227                 for(j = 0; j < n; j++)
00228                         jac[i * n + j] = zero;
00229 
00230         if( n <= m )
00231         {       // use forward mode ----------------------------------------
00232         
00233                 // initial coloring
00234                 SizeVector color(n);
00235                 for(j = 0; j < n; j++)
00236                         color[j] = j;
00237 
00238                 // See GreedyPartialD2Coloring Algorithm Section 3.6.2 of
00239                 // Graph Coloring in Optimization Revisited by
00240                 // Assefaw Gebremedhin, Fredrik Maane, Alex Pothen
00241                 VectorBool forbidden(n);
00242                 for(j = 0; j < n; j++)
00243                 {       // initial all colors as ok for this column
00244                         for(k = 0; k < n; k++)
00245                                 forbidden[k] = false;
00246                         // for each row that is connected to column j
00247                         for(i = 0; i < m; i++) if( p[i * n + j] )
00248                         {       // for each column that is connected to row i
00249                                 for(k = 0; k < n; k++)
00250                                         if( p[i * n + k] & (j != k) )   
00251                                                 forbidden[ color[k] ] = true;
00252                         }
00253                         k = 0;
00254                         while( forbidden[k] && k < n )
00255                         {       k++;
00256                                 CPPAD_ASSERT_UNKNOWN( k < n );
00257                         }
00258                         color[j] = k;
00259                 }
00260                 size_t n_color = 1;
00261                 for(k = 0; k < n; k++) 
00262                         n_color = std::max(n_color, color[k] + 1);
00263 
00264                 // direction vector for calls to forward
00265                 VectorBase dx(n);
00266 
00267                 // location for return values from Reverse
00268                 VectorBase dy(m);
00269 
00270                 // loop over colors
00271                 size_t c;
00272                 for(c = 0; c < n_color; c++)
00273                 {       // determine all the colums with this color
00274                         for(j = 0; j < n; j++)
00275                         {       if( color[j] == c )
00276                                         dx[j] = one;
00277                                 else    dx[j] = zero;
00278                         }
00279                         // call forward mode for all these columns at once
00280                         dy = Forward(1, dx);
00281 
00282                         // set the corresponding components of the result
00283                         for(j = 0; j < n; j++) if( color[j] == c )
00284                         {       for(i = 0; i < m; i++) 
00285                                         if( p[ i * n + j ] )
00286                                                 jac[i * n + j] = dy[i];
00287                         }
00288                 }
00289         }
00290         else
00291         {       // use reverse mode ----------------------------------------
00292         
00293                 // initial coloring
00294                 SizeVector color(m);
00295                 for(i = 0; i < m; i++)
00296                         color[i] = i;
00297 
00298                 // See GreedyPartialD2Coloring Algorithm Section 3.6.2 of
00299                 // Graph Coloring in Optimization Revisited by
00300                 // Assefaw Gebremedhin, Fredrik Maane, Alex Pothen
00301                 VectorBool forbidden(m);
00302                 for(i = 0; i < m; i++)
00303                 {       // initial all colors as ok for this row
00304                         for(k = 0; k < m; k++)
00305                                 forbidden[k] = false;
00306                         // for each column that is connected to row i
00307                         for(j = 0; j < n; j++) if( p[i * n + j] )
00308                         {       // for each row that is connected to column j
00309                                 for(k = 0; k < m; k++)
00310                                         if( p[k * n + j] & (i != k) )   
00311                                                 forbidden[ color[k] ] = true;
00312                         }
00313                         k = 0;
00314                         while( forbidden[k] && k < m )
00315                         {       k++;
00316                                 CPPAD_ASSERT_UNKNOWN( k < n );
00317                         }
00318                         color[i] = k;
00319                 }
00320                 size_t n_color = 1;
00321                 for(k = 0; k < m; k++) 
00322                         n_color = std::max(n_color, color[k] + 1);
00323 
00324                 // weight vector for calls to reverse
00325                 VectorBase w(m);
00326 
00327                 // location for return values from Reverse
00328                 VectorBase dw(n);
00329 
00330                 // loop over colors
00331                 size_t c;
00332                 for(c = 0; c < n_color; c++)
00333                 {       // determine all the rows with this color
00334                         for(i = 0; i < m; i++)
00335                         {       if( color[i] == c )
00336                                         w[i] = one;
00337                                 else    w[i] = zero;
00338                         }
00339                         // call reverse mode for all these rows at once
00340                         dw = Reverse(1, w);
00341 
00342                         // set the corresponding components of the result
00343                         for(i = 0; i < m; i++) if( color[i] == c )
00344                         {       for(j = 0; j < n; j++) 
00345                                         if( p[ i * n + j ] )
00346                                                 jac[i * n + j] = dw[j];
00347                         }
00348                 }
00349         }
00350 }
00351 
00352 /*!
00353 Private helper function for SparseJacobian(x, p).
00354 
00355 All of the description in the public member function SparseJacobian(x, p)
00356 applies.
00357 
00358 \param set_type
00359 is a \c std::set<size_t> value.
00360 This argument is used to dispatch to the proper souce
00361 code depending on the vlaue of \c VectorSet::value_type.
00362 
00363 \param x
00364 See \c SparseJacobian(x, p).
00365 
00366 \param p
00367 See \c SparseJacobian(x, p).
00368 
00369 \param jac
00370 is the return value for the corresponding call to \c SparseJacobian(x, p).
00371 On input it must have size equalt to the domain times range dimension
00372 for this ADFun<Base> object.
00373 On return, it will contain the Jacobian.
00374 */
00375 template <class Base>
00376 template <class VectorBase, class VectorSet>
00377 void ADFun<Base>::SparseJacobianCase(
00378         const std::set<size_t>&  set_type        ,
00379         const VectorBase&        x               , 
00380         const VectorSet&         p               ,
00381         VectorBase&              jac             )
00382 {
00383         typedef CppAD::vector<size_t> SizeVector;
00384         typedef CppAD::vectorBool     VectorBool;
00385         size_t i, j, k;
00386 
00387         size_t m = Range();
00388         size_t n = Domain();
00389 
00390         // some values
00391         const Base zero(0);
00392         const Base one(1);
00393 
00394         // check VectorSet is Simple Vector class with sets for elements
00395         static std::set<size_t> two, three;
00396         if( two.empty() )
00397         {       two.insert(2);
00398                 three.insert(3);
00399         }
00400         CPPAD_ASSERT_UNKNOWN( two.size() == 1 );
00401         CPPAD_ASSERT_UNKNOWN( three.size() == 1 );
00402         CheckSimpleVector<std::set<size_t>, VectorSet>(two, three);
00403 
00404         // check VectorBase is Simple Vector class with Base type elements
00405         CheckSimpleVector<Base, VectorBase>();
00406 
00407         CPPAD_ASSERT_KNOWN(
00408                 x.size() == n,
00409                 "SparseJacobian: size of x not equal domain dimension for f"
00410         ); 
00411         CPPAD_ASSERT_KNOWN(
00412                 p.size() == m,
00413                 "SparseJacobian: using sets and size of p "
00414                 "not equal range dimension for f"
00415         ); 
00416         CPPAD_ASSERT_UNKNOWN(jac.size() == m * n); 
00417 
00418         // point at which we are evaluating the Jacobian
00419         Forward(0, x);
00420 
00421         // initialize the return value
00422         for(i = 0; i < m; i++)
00423                 for(j = 0; j < n; j++)
00424                         jac[i * n + j] = zero;
00425 
00426         // create a copy of the transpose sparsity pattern
00427         VectorSet q(n);
00428         std::set<size_t>::const_iterator itr_i, itr_j;
00429         for(i = 0; i < m; i++)
00430         {       itr_j = p[i].begin();
00431                 while( itr_j != p[i].end() )
00432                 {       j = *itr_j++;
00433                         q[j].insert(i);
00434                 }
00435         }       
00436 
00437         if( n <= m )
00438         {       // use forward mode ----------------------------------------
00439         
00440                 // initial coloring
00441                 SizeVector color(n);
00442                 for(j = 0; j < n; j++)
00443                         color[j] = j;
00444 
00445                 // See GreedyPartialD2Coloring Algorithm Section 3.6.2 of
00446                 // Graph Coloring in Optimization Revisited by
00447                 // Assefaw Gebremedhin, Fredrik Maane, Alex Pothen
00448                 VectorBool forbidden(n);
00449                 for(j = 0; j < n; j++)
00450                 {       // initial all colors as ok for this column
00451                         for(k = 0; k < n; k++)
00452                                 forbidden[k] = false;
00453 
00454                         // for each row connected to column j
00455                         itr_i = q[j].begin();
00456                         while( itr_i != q[j].end() )
00457                         {       i = *itr_i++;
00458                                 // for each column connected to row i
00459                                 itr_j = p[i].begin();
00460                                 while( itr_j != p[i].end() )
00461                                 {       // if this is not j, forbid it
00462                                         k = *itr_j++;
00463                                         forbidden[ color[k] ] = (k != j);
00464                                 }
00465                         }
00466                         k = 0;
00467                         while( forbidden[k] && k < n )
00468                         {       k++;
00469                                 CPPAD_ASSERT_UNKNOWN( k < n );
00470                         }
00471                         color[j] = k;
00472                 }
00473                 size_t n_color = 1;
00474                 for(k = 0; k < n; k++) 
00475                         n_color = std::max(n_color, color[k] + 1);
00476 
00477                 // direction vector for calls to forward
00478                 VectorBase dx(n);
00479 
00480                 // location for return values from Reverse
00481                 VectorBase dy(m);
00482 
00483                 // loop over colors
00484                 size_t c;
00485                 for(c = 0; c < n_color; c++)
00486                 {       // determine all the colums with this color
00487                         for(j = 0; j < n; j++)
00488                         {       if( color[j] == c )
00489                                         dx[j] = one;
00490                                 else    dx[j] = zero;
00491                         }
00492                         // call forward mode for all these columns at once
00493                         dy = Forward(1, dx);
00494 
00495                         // set the corresponding components of the result
00496                         for(j = 0; j < n; j++) if( color[j] == c )
00497                         {       itr_i = q[j].begin();
00498                                 while( itr_i != q[j].end() )
00499                                 {       i = *itr_i++;
00500                                         jac[i * n + j] = dy[i];
00501                                 }
00502                         }
00503                 }
00504         }
00505         else
00506         {       // use reverse mode ----------------------------------------
00507         
00508                 // initial coloring
00509                 SizeVector color(m);
00510                 for(i = 0; i < m; i++)
00511                         color[i] = i;
00512 
00513                 // See GreedyPartialD2Coloring Algorithm Section 3.6.2 of
00514                 // Graph Coloring in Optimization Revisited by
00515                 // Assefaw Gebremedhin, Fredrik Maane, Alex Pothen
00516                 VectorBool forbidden(m);
00517                 for(i = 0; i < m; i++)
00518                 {       // initial all colors as ok for this row
00519                         for(k = 0; k < m; k++)
00520                                 forbidden[k] = false;
00521 
00522                         // for each column connected to row i
00523                         itr_j = p[i].begin();
00524                         while( itr_j != p[i].end() )
00525                         {       j = *itr_j++;   
00526                                 // for each row connected to column j
00527                                 itr_i = q[j].begin();
00528                                 while( itr_i != q[j].end() )
00529                                 {       // if this is not i, forbid it
00530                                         k = *itr_i++;
00531                                         forbidden[ color[k] ] = (k != i);
00532                                 }
00533                         }
00534                         k = 0;
00535                         while( forbidden[k] && k < m )
00536                         {       k++;
00537                                 CPPAD_ASSERT_UNKNOWN( k < n );
00538                         }
00539                         color[i] = k;
00540                 }
00541                 size_t n_color = 1;
00542                 for(k = 0; k < m; k++) 
00543                         n_color = std::max(n_color, color[k] + 1);
00544 
00545                 // weight vector for calls to reverse
00546                 VectorBase w(m);
00547 
00548                 // location for return values from Reverse
00549                 VectorBase dw(n);
00550 
00551                 // loop over colors
00552                 size_t c;
00553                 for(c = 0; c < n_color; c++)
00554                 {       // determine all the rows with this color
00555                         for(i = 0; i < m; i++)
00556                         {       if( color[i] == c )
00557                                         w[i] = one;
00558                                 else    w[i] = zero;
00559                         }
00560                         // call reverse mode for all these rows at once
00561                         dw = Reverse(1, w);
00562 
00563                         // set the corresponding components of the result
00564                         for(i = 0; i < m; i++) if( color[i] == c )
00565                         {       itr_j = p[i].begin();
00566                                 while( itr_j != p[i].end() )
00567                                 {       j = *itr_j++;
00568                                         jac[i * n + j] = dw[j];
00569                                 }
00570                         }
00571                 }
00572         }
00573 }
00574 
00575 /*!
00576 Compute a sparse Jacobian.
00577 
00578 The C++ source code corresponding to this operation is
00579 \verbatim
00580         jac = SparseJacobian(x, p)
00581 \endverbatim
00582 
00583 \tparam Base
00584 is the base type for the recording that is stored in this
00585 ADFun<Base object.
00586 
00587 \tparam VectorBase
00588 is a simple vector class with elements of type \a Base.
00589 
00590 \tparam VectorSet
00591 is a simple vector class with elements of type 
00592 \c bool or \c std::set<size_t>.
00593 
00594 \param x
00595 is a vector specifing the point at which to compute the Jacobian.
00596 
00597 \param p
00598 is the sparsity pattern for the Jacobian that we are calculating.
00599 
00600 \return
00601 Will be a vector if size \c m * n containing the Jacobian at the
00602 specified point (in row major order).
00603 */
00604 template <class Base>
00605 template <class VectorBase, class VectorSet>
00606 VectorBase ADFun<Base>::SparseJacobian(
00607         const VectorBase& x, const VectorSet& p
00608 )
00609 {       size_t m = Range();
00610         size_t n = Domain();
00611         VectorBase jac(m * n);
00612 
00613         typedef typename VectorSet::value_type Set_type;
00614         SparseJacobianCase( Set_type(), x, p, jac);
00615 
00616         return jac;
00617 } 
00618 
00619 /*!
00620 Compute a sparse Jacobian.
00621 
00622 The C++ source code corresponding to this operation is
00623 \verbatim
00624         jac = SparseJacobian(x)
00625 \endverbatim
00626 
00627 \tparam Base
00628 is the base type for the recording that is stored in this
00629 ADFun<Base object.
00630 
00631 \tparam VectorBase
00632 is a simple vector class with elements of the \a Base.
00633 
00634 \param x
00635 is a vector specifing the point at which to compute the Jacobian.
00636 
00637 \return
00638 Will be a vector of size \c m * n containing the Jacobian at the
00639 specified point (in row major order).
00640 */
00641 template <class Base>
00642 template <class VectorBase>
00643 VectorBase ADFun<Base>::SparseJacobian( const VectorBase& x )
00644 {       typedef CppAD::vectorBool   VectorBool;
00645 
00646         size_t m = Range();
00647         size_t n = Domain();
00648 
00649         // sparsity pattern for Jacobian
00650         VectorBool p(m * n);
00651 
00652         // return result
00653         VectorBase jac(m * n);
00654 
00655         if( n <= m )
00656         {       size_t j, k;
00657 
00658                 // use forward mode 
00659                 VectorBool r(n * n);
00660                 for(j = 0; j < n; j++)
00661                 {       for(k = 0; k < n; k++)
00662                                 r[j * n + k] = false;
00663                         r[j * n + j] = true;
00664                 }
00665                 p = ForSparseJac(n, r);
00666         }
00667         else
00668         {       size_t i, k;
00669 
00670                 // use reverse mode 
00671                 VectorBool s(m * m);
00672                 for(i = 0; i < m; i++)
00673                 {       for(k = 0; k < m; k++)
00674                                 s[i * m + k] = false;
00675                         s[i * m + i] = true;
00676                 }
00677                 p = RevSparseJac(m, s);
00678         }
00679         bool set_type = true; // only used to dispatch compiler to proper case
00680         SparseJacobianCase(set_type, x, p, jac);
00681         return jac;
00682 }
00683 
00684 
00685 CPPAD_END_NAMESPACE
00686 # endif