CppAD: A C++ Algorithmic Differentiation Package  20130918
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-14 Bradley M. Bell
00007 
00008 CppAD is distributed under multiple licenses. This distribution is under
00009 the terms of the 
00010                     Eclipse 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 $begin sparse_jacobian$$
00017 $spell
00018      cppad
00019      colpack
00020      cmake
00021      recomputed
00022      valarray
00023      std
00024      CppAD
00025      Bool
00026      jac
00027      Jacobian
00028      const
00029      Taylor
00030 $$
00031 
00032 $section Sparse Jacobian: Easy Driver$$
00033 $index SparseJacobian$$
00034 $index jacobian, sparse$$
00035 
00036 $head Syntax$$
00037 $icode%jac% = %f%.SparseJacobian(%x%)
00038 %jac% = %f%.SparseJacobian(%x%, %p%)
00039 %n_sweep% = %f%.SparseJacobianForward(%x%, %p%, %row%, %col%, %jac%, %work%)
00040 %n_sweep% = %f%.SparseJacobianReverse(%x%, %p%, %row%, %col%, %jac%, %work%)
00041 %$$
00042 
00043 $head Purpose$$
00044 We use $latex n$$ for the $cref/domain/seq_property/Domain/$$ size,
00045 and $latex m$$ for the $cref/range/seq_property/Range/$$ size of $icode f$$.
00046 We use $latex F : \B{R}^n \rightarrow \B{R}^m$$ do denote the
00047 $cref/AD function/glossary/AD Function/$$
00048 corresponding to $icode f$$. 
00049 The syntax above sets $icode jac$$ to the Jacobian 
00050 $latex \[
00051      jac = F^{(1)} (x) 
00052 \] $$
00053 This routine takes advantage of the sparsity of the Jacobian
00054 in order to reduce the amount of computation necessary.
00055 If $icode row$$ and $icode col$$ are present, it also takes
00056 advantage of the reduced set of elements of the Jacobian that
00057 need to be computed.
00058 One can use speed tests (e.g. $cref speed_test$$)
00059 to verify that results are computed faster
00060 than when using the routine $cref Jacobian$$.
00061 
00062 $head f$$
00063 The object $icode f$$ has prototype
00064 $codei%
00065      ADFun<%Base%> %f%
00066 %$$
00067 Note that the $cref ADFun$$ object $icode f$$ is not $code const$$
00068 (see $cref/Uses Forward/sparse_jacobian/Uses Forward/$$ below).
00069 
00070 $head x$$
00071 The argument $icode x$$ has prototype
00072 $codei%
00073      const %VectorBase%& %x%
00074 %$$
00075 (see $cref/VectorBase/sparse_jacobian/VectorBase/$$ below)
00076 and its size 
00077 must be equal to $icode n$$, the dimension of the
00078 $cref/domain/seq_property/Domain/$$ space for $icode f$$.
00079 It specifies
00080 that point at which to evaluate the Jacobian.
00081 
00082 $head p$$
00083 The argument $icode p$$ is optional and has prototype
00084 $codei%
00085      const %VectorSet%& %p%
00086 %$$
00087 (see $cref/VectorSet/sparse_jacobian/VectorSet/$$ below).
00088 If it has elements of type $code bool$$,
00089 its size is $latex m * n$$.
00090 If it has elements of type $code std::set<size_t>$$,
00091 its size is $latex m$$ and all its set elements are between
00092 zero and $latex n - 1$$.
00093 It specifies a 
00094 $cref/sparsity pattern/glossary/Sparsity Pattern/$$ 
00095 for the Jacobian $latex F^{(1)} (x)$$.
00096 $pre
00097 
00098 $$
00099 If this sparsity pattern does not change between calls to 
00100 $codei SparseJacobian$$, it should be faster to calculate $icode p$$ once 
00101 (using $cref ForSparseJac$$ or $cref RevSparseJac$$)
00102 and then pass $icode p$$ to $codei SparseJacobian$$.
00103 In addition,
00104 if you specify $icode p$$, CppAD will use the same
00105 type of sparsity representation 
00106 (vectors of $code bool$$ or vectors of $code std::set<size_t>$$)
00107 for its internal calculations.
00108 Otherwise, the representation
00109 for the internal calculations is unspecified.
00110 
00111 $head row, col$$
00112 The arguments $icode row$$ and $icode col$$ are optional and have prototype
00113 $codei%
00114      const %VectorSize%& %row%
00115      const %VectorSize%& %col%
00116 %$$
00117 (see $cref/VectorSize/sparse_jacobian/VectorSize/$$ below).
00118 They specify which rows and columns of $latex F^{(1)} (x)$$ are
00119 computes and in what order.
00120 Not all the non-zero entries in $latex F^{(1)} (x)$$ need be computed,
00121 but all the entries specified by $icode row$$ and $icode col$$
00122 must be possibly non-zero in the sparsity pattern.
00123 We use $latex K$$ to denote the value $icode%jac%.size()%$$
00124 which must also equal the size of $icode row$$ and $icode col$$.
00125 Furthermore,
00126 for $latex k = 0 , \ldots , K-1$$, it must hold that
00127 $latex row[k] < m$$ and $latex col[k] < n$$.
00128 
00129 $head jac$$
00130 The result $icode jac$$ has prototype
00131 $codei%
00132      %VectorBase%& %jac%
00133 %$$
00134 In the case where the arguments $icode row$$ and $icode col$$ are not present,
00135 the size of $icode jac$$ is $latex m * n$$ and
00136 for $latex i = 0 , \ldots , m-1$$,
00137 $latex j = 0 , \ldots , n-1$$,
00138 $latex \[
00139      jac [ i * n + j ] = \D{ F_i }{ x_j } (x)
00140 \] $$
00141 $pre
00142 
00143 $$
00144 In the case where the arguments $icode row$$ and $icode col$$ are present,
00145 we use $latex K$$ to denote the size of $icode jac$$. 
00146 The input value of its elements does not matter.
00147 Upon return, for $latex k = 0 , \ldots , K - 1$$,
00148 $latex \[
00149      jac [ k ] = \D{ F_i }{ x_j } (x)
00150      \; , \;
00151      \; {\rm where} \;
00152      i = row[k]
00153      \; {\rm and } \;
00154      j = col[k]
00155 \] $$
00156 
00157 $head work$$
00158 If this argument is present, it has prototype
00159 $codei%
00160      sparse_jacobian_work& %work%
00161 %$$
00162 This object can only be used with the routines 
00163 $code SparseJacobianForward$$ and $code SparseJacobianReverse$$.
00164 During its the first use, information is stored in $icode work$$.
00165 This is used to reduce the work done by future calls to the same mode
00166 (forward or reverse), 
00167 the same $icode f$$, $icode p$$, $icode row$$, and $icode col$$.
00168 If a future call is for a different mode,
00169 or any of these values have changed,
00170 you must first call $icode%work%.clear()%$$
00171 to inform CppAD that this information needs to be recomputed.
00172 
00173 $subhead color_method$$
00174 The coloring algorithm determines which columns (forward mode)
00175 or rows (reverse mode) can be computed during the same sweep.
00176 This field has prototype
00177 $codep%
00178      std::string %work%.color_method
00179 %$$
00180 and its default value (after a constructor or $code clear()$$) 
00181 is $code "cppad"$$.
00182 If $cref colpack_prefix$$ is specified on the
00183 $cref/cmake command/cmake/CMake Command/$$ line,
00184 you can set this method to $code "colpack"$$.
00185 This value only matters on the first call to $code sparse_jacobian$$
00186 after the $icode work$$ constructor or a call to $code clear$$.
00187 
00188 $head n_sweep$$
00189 The return value $icode n_sweep$$ has prototype
00190 $codei%
00191      size_t %n_sweep%
00192 %$$
00193 If $code SparseJacobianForward$$ ($code SparseJacobianReverse$$) is used, 
00194 $icode n_sweep$$ is the number of first order forward (reverse) sweeps 
00195 used to compute the requested Jacobian values. 
00196 (This is also the number of colors determined by the coloring method
00197 mentioned above).
00198 This is proportional to the total work that $code SparseJacobian$$ does, 
00199 not counting the zero order forward sweep, 
00200 or the work to combine multiple columns (rows) into a single sweep.
00201 
00202 $head VectorBase$$
00203 The type $icode VectorBase$$ must be a $cref SimpleVector$$ class with
00204 $cref/elements of type/SimpleVector/Elements of Specified Type/$$
00205 $icode Base$$.
00206 The routine $cref CheckSimpleVector$$ will generate an error message
00207 if this is not the case.
00208 
00209 $head VectorSet$$
00210 The type $icode VectorSet$$ must be a $cref SimpleVector$$ class with
00211 $cref/elements of type/SimpleVector/Elements of Specified Type/$$
00212 $code bool$$ or $code std::set<size_t>$$;
00213 see $cref/sparsity pattern/glossary/Sparsity Pattern/$$ for a discussion
00214 of the difference.
00215 The routine $cref CheckSimpleVector$$ will generate an error message
00216 if this is not the case.
00217 
00218 $subhead Restrictions$$
00219 If $icode VectorSet$$ has elements of $code std::set<size_t>$$,
00220 then $icode%p%[%i%]%$$ must return a reference (not a copy) to the 
00221 corresponding set.
00222 According to section 26.3.2.3 of the 1998 C++ standard,
00223 $code std::valarray< std::set<size_t> >$$ does not satisfy
00224 this condition. 
00225 
00226 $head VectorSize$$
00227 The type $icode VectorSize$$ must be a $cref SimpleVector$$ class with
00228 $cref/elements of type/SimpleVector/Elements of Specified Type/$$
00229 $code size_t$$.
00230 The routine $cref CheckSimpleVector$$ will generate an error message
00231 if this is not the case.
00232 
00233 $head Uses Forward$$
00234 After each call to $cref Forward$$,
00235 the object $icode f$$ contains the corresponding 
00236 $cref/Taylor coefficients/glossary/Taylor Coefficient/$$.
00237 After a call to any of the sparse Jacobian routines,
00238 the zero order Taylor coefficients correspond to
00239 $icode%f%.Forward(0, %x%)%$$
00240 and the other coefficients are unspecified.
00241 
00242 After $code SparseJacobian$$,
00243 the previous calls to $cref Forward$$ are undefined.
00244 
00245 $head Example$$
00246 $children%
00247      example/sparse_jacobian.cpp
00248 %$$
00249 The routine
00250 $cref sparse_jacobian.cpp$$
00251 is examples and tests of $code sparse_jacobian$$.
00252 It return $code true$$, if it succeeds and $code false$$ otherwise.
00253 
00254 $end
00255 ==============================================================================
00256 */
00257 # include <cppad/local/std_set.hpp>
00258 # include <cppad/local/color_general.hpp>
00259 
00260 namespace CppAD { // BEGIN_CPPAD_NAMESPACE
00261 /*!
00262 \file sparse_jacobian.hpp
00263 Sparse Jacobian driver routine and helper functions.
00264 */
00265 // ===========================================================================
00266 /*!
00267 class used by SparseJacobian to hold information so it does not need to be 
00268 recomputed.
00269 */
00270 class sparse_jacobian_work {
00271      public:
00272           /// Coloring method: color_general_cppad or color_general_colpack
00273           /// (this field is set by user)
00274           std::string color_method;
00275           /// indices that sort the user row and col arrays by color 
00276           CppAD::vector<size_t> order;
00277           /// results of the coloring algorithm
00278           CppAD::vector<size_t> color;
00279      
00280           /// constructor
00281           sparse_jacobian_work(void) : color_method("cppad")
00282           { }
00283           /// reset coloring method to its default and
00284           /// inform CppAD that color and order need to be recomputed
00285           void clear(void)
00286           {    color_method = "cppad";
00287                order.clear();
00288                color.clear();
00289           }
00290 };
00291 // ===========================================================================
00292 /*!
00293 Private helper function forward mode cases
00294 
00295 \tparam Base
00296 is the base type for the recording that is stored in this
00297 <code>ADFun<Base></code> object.
00298 
00299 \tparam VectorBase
00300 is a simple vector class with elements of type \a Base.
00301 
00302 \tparam VectorSet
00303 is either \c sparse_pack, \c sparse_set or \c sparse_list.
00304 
00305 \tparam VectorSize
00306 is a simple vector class with elements of type \c size_t.
00307 
00308 \param x [in]
00309 is a vector specifing the point at which to compute the Jacobian.
00310 
00311 \param p_transpose [in]
00312 If <code>work.color.size() != 0</code>, 
00313 then \c p_transpose is not used.
00314 Otherwise, it is a
00315 sparsity pattern for the transpose of the Jacobian of this ADFun<Base> object.
00316 Note that we do not change the values in \c p_transpose,
00317 but is not \c const because we use its iterator facility.
00318 
00319 \param row [in]
00320 is the vector of row indices for the returned Jacobian values.
00321 
00322 \param col [in]
00323 is the vector of columns indices for the returned Jacobian values.
00324 It must have the same size as \c row.
00325 
00326 \param jac [out]
00327 is the vector of Jacobian values. We use \c K to denote the size of \c jac.
00328 The return value <code>jac[k]</code> is the partial of the
00329 <code>row[k]</code> range component of the function with respect
00330 the the <code>col[k]</code> domain component of its argument.
00331 
00332 \param work
00333 <code>work.color_method</code> is an input. The rest of
00334 this structure contains information that is computed by \c SparseJacobainFor.
00335 If the sparsity pattern, \c row vector, or \c col vectors
00336 are not the same between calls to \c SparseJacobianFor, 
00337 \c work.clear() must be called to reinitialize \c work.
00338 
00339 \return
00340 Is the number of first order forward sweeps used to compute the
00341 requested Jacobian values. The total work, not counting the zero order
00342 forward sweep, or the time to combine computations, is proportional to this
00343 return value.
00344 */
00345 template<class Base>
00346 template <class VectorBase, class VectorSet, class VectorSize>
00347 size_t ADFun<Base>::SparseJacobianFor(
00348      const VectorBase&            x           ,
00349            VectorSet&             p_transpose ,
00350      const VectorSize&            row         ,
00351      const VectorSize&            col         ,
00352            VectorBase&            jac         ,
00353             sparse_jacobian_work& work        )
00354 {
00355      size_t j, k, ell;
00356 
00357      CppAD::vector<size_t>& order(work.order);
00358      CppAD::vector<size_t>& color(work.color);
00359 
00360      size_t m = Range();
00361      size_t n = Domain();
00362 
00363      // some values
00364      const Base zero(0);
00365      const Base one(1);
00366 
00367      // check VectorBase is Simple Vector class with Base type elements
00368      CheckSimpleVector<Base, VectorBase>();
00369 
00370      CPPAD_ASSERT_UNKNOWN( size_t(x.size()) == n );
00371      CPPAD_ASSERT_UNKNOWN( color.size() == 0 || color.size() == n );
00372 
00373      // number of components of Jacobian that are required
00374      size_t K = size_t(jac.size());
00375      CPPAD_ASSERT_UNKNOWN( row.size() == K );
00376      CPPAD_ASSERT_UNKNOWN( col.size() == K );
00377 
00378      // Point at which we are evaluating the Jacobian
00379      Forward(0, x);
00380 
00381      // check for case where nothing (except Forward above) to do
00382      if( K == 0 )
00383           return 0;
00384 
00385      if( color.size() == 0 )
00386      {
00387           CPPAD_ASSERT_UNKNOWN( p_transpose.n_set() ==  n );
00388           CPPAD_ASSERT_UNKNOWN( p_transpose.end() ==  m );
00389 
00390           // execute coloring algorithm
00391           color.resize(n);
00392           if(  work.color_method == "cppad" )
00393                color_general_cppad(p_transpose, col, row, color);
00394           else if( work.color_method == "colpack" )
00395           {
00396 # if CPPAD_HAS_COLPACK
00397                color_general_colpack(p_transpose, col, row, color);
00398 # else
00399                CPPAD_ASSERT_KNOWN(
00400                     false,
00401                     "SparseJacobianForward: work.color_method = colpack "
00402                     "and colpack_prefix missing from cmake command line."
00403                );
00404 # endif
00405           }
00406           else CPPAD_ASSERT_KNOWN(
00407                false,
00408                "SparseJacobianForward: work.color_method is not valid."
00409           );
00410 
00411           // put sorting indices in color order
00412           VectorSize key(K);
00413           order.resize(K);
00414           for(k = 0; k < K; k++)
00415                key[k] = color[ col[k] ];
00416           index_sort(key, order);
00417      }
00418      size_t n_color = 1;
00419      for(j = 0; j < n; j++) if( color[j] < n )
00420           n_color = std::max(n_color, color[j] + 1);
00421 
00422      // direction vector for calls to forward
00423      VectorBase dx(n);
00424 
00425      // location for return values from forward
00426      VectorBase dy(m);
00427 
00428      // initialize the return value
00429      for(k = 0; k < K; k++)
00430           jac[k] = zero;
00431 
00432      // loop over colors
00433      k = 0;
00434      for(ell = 0; ell < n_color; ell++)
00435      {    CPPAD_ASSERT_UNKNOWN( color[ col[ order[k] ] ] == ell );
00436 
00437           // combine all columns with this color
00438           for(j = 0; j < n; j++)
00439           {    dx[j] = zero;
00440                if( color[j] == ell )
00441                     dx[j] = one;
00442           }
00443           // call forward mode for all these columns at once
00444           dy = Forward(1, dx);
00445 
00446           // set the corresponding components of the result
00447           while( k < K && color[ col[order[k]] ] == ell ) 
00448           {    jac[ order[k] ] = dy[row[order[k]]];
00449                k++;
00450           }
00451      }
00452      return n_color;
00453 }
00454 /*!
00455 Private helper function for reverse mode cases.
00456 
00457 \tparam Base
00458 is the base type for the recording that is stored in this
00459 <code>ADFun<Base></code> object.
00460 
00461 \tparam VectorBase
00462 is a simple vector class with elements of type \a Base.
00463 
00464 \tparam VectorSet
00465 is either \c sparse_pack, \c sparse_set or \c sparse_list.
00466 
00467 \tparam VectorSize
00468 is a simple vector class with elements of type \c size_t.
00469 
00470 \param x [in]
00471 is a vector specifing the point at which to compute the Jacobian.
00472 
00473 \param p [in]
00474 If <code>work.color.size() != 0</code>, then \c p is not used.
00475 Otherwise, it is a
00476 sparsity pattern for the Jacobian of this ADFun<Base> object.
00477 Note that we do not change the values in \c p,
00478 but is not \c const because we use its iterator facility.
00479 
00480 \param row [in]
00481 is the vector of row indices for the returned Jacobian values.
00482 
00483 \param col [in]
00484 is the vector of columns indices for the returned Jacobian values.
00485 It must have the same size as \c row.
00486 
00487 \param jac [out]
00488 is the vector of Jacobian values.
00489 It must have the same size as \c row. 
00490 The return value <code>jac[k]</code> is the partial of the
00491 <code>row[k]</code> range component of the function with respect
00492 the the <code>col[k]</code> domain component of its argument.
00493 
00494 \param work
00495 <code>work.color_method</code> is an input. The rest of
00496 This structure contains information that is computed by \c SparseJacobainRev.
00497 If the sparsity pattern, \c row vector, or \c col vectors
00498 are not the same between calls to \c SparseJacobianRev, 
00499 \c work.clear() must be called to reinitialize \c work.
00500 
00501 \return
00502 Is the number of first order reverse sweeps used to compute the
00503 reverse Jacobian values. The total work, not counting the zero order
00504 forward sweep, or the time to combine computations, is proportional to this
00505 return value.
00506 */
00507 template<class Base>
00508 template <class VectorBase, class VectorSet, class VectorSize>
00509 size_t ADFun<Base>::SparseJacobianRev(
00510      const VectorBase&           x           ,
00511            VectorSet&            p           ,
00512      const VectorSize&           row         ,
00513      const VectorSize&           col         ,
00514            VectorBase&           jac         ,
00515            sparse_jacobian_work& work        )
00516 {
00517      size_t i, k, ell;
00518 
00519      CppAD::vector<size_t>& order(work.order);
00520      CppAD::vector<size_t>& color(work.color);
00521 
00522      size_t m = Range();
00523      size_t n = Domain();
00524 
00525      // some values
00526      const Base zero(0);
00527      const Base one(1);
00528 
00529      // check VectorBase is Simple Vector class with Base type elements
00530      CheckSimpleVector<Base, VectorBase>();
00531 
00532      CPPAD_ASSERT_UNKNOWN( size_t(x.size()) == n );
00533      CPPAD_ASSERT_UNKNOWN (color.size() == m || color.size() == 0 );
00534 
00535      // number of components of Jacobian that are required
00536      size_t K = size_t(jac.size());
00537      CPPAD_ASSERT_UNKNOWN( row.size() == K );
00538      CPPAD_ASSERT_UNKNOWN( col.size() == K );
00539 
00540      // Point at which we are evaluating the Jacobian
00541      Forward(0, x);
00542 
00543      // check for case where nothing (except Forward above) to do
00544      if( K == 0 )
00545           return 0;
00546 
00547      if( color.size() == 0 )
00548      {
00549           CPPAD_ASSERT_UNKNOWN( p.n_set() == m );
00550           CPPAD_ASSERT_UNKNOWN( p.end()   == n );
00551 
00552           // execute the coloring algorithm
00553           color.resize(m);
00554           if(  work.color_method == "cppad" )
00555                color_general_cppad(p, row, col, color);
00556           else if( work.color_method == "colpack" )
00557           {
00558 # if CPPAD_HAS_COLPACK
00559                color_general_colpack(p, row, col, color);
00560 # else
00561                CPPAD_ASSERT_KNOWN(
00562                     false,
00563                     "SparseJacobianReverse: work.color_method = colpack "
00564                     "and colpack_prefix missing from cmake command line."
00565                );
00566 # endif
00567           }
00568           else CPPAD_ASSERT_KNOWN(
00569                false,
00570                "SparseJacobianReverse: work.color_method is not valid."
00571           );
00572 
00573           // put sorting indices in color order
00574           VectorSize key(K);
00575           order.resize(K);
00576           for(k = 0; k < K; k++)
00577                key[k] = color[ row[k] ];
00578           index_sort(key, order);
00579      }
00580      size_t n_color = 1;
00581      for(i = 0; i < m; i++) if( color[i] < m ) 
00582           n_color = std::max(n_color, color[i] + 1);
00583 
00584      // weighting vector for calls to reverse
00585      VectorBase w(m);
00586 
00587      // location for return values from Reverse
00588      VectorBase dw(n);
00589 
00590      // initialize the return value
00591      for(k = 0; k < K; k++)
00592           jac[k] = zero;
00593 
00594      // loop over colors
00595      k = 0;
00596      for(ell = 0; ell < n_color; ell++)
00597      {    CPPAD_ASSERT_UNKNOWN( color[ row[ order[k] ] ] == ell );
00598 
00599           // combine all the rows with this color
00600           for(i = 0; i < m; i++)
00601           {    w[i] = zero;
00602                if( color[i] == ell )
00603                     w[i] = one;
00604           }
00605           // call reverse mode for all these rows at once
00606           dw = Reverse(1, w);
00607 
00608           // set the corresponding components of the result
00609           while( k < K && color[ row[order[k]] ]  == ell ) 
00610           {    jac[ order[k] ] = dw[col[order[k]]];
00611                k++;
00612           }
00613      }
00614      return n_color;
00615 }
00616 // ==========================================================================
00617 // Public Member functions
00618 // ==========================================================================
00619 /*!
00620 Compute user specified subset of a sparse Jacobian using forward mode.
00621 
00622 The C++ source code corresponding to this operation is
00623 \verbatim
00624      SparceJacobianForward(x, p, row, col, jac, work)
00625 \endverbatim
00626 
00627 \tparam Base
00628 is the base type for the recording that is stored in this
00629 <code>ADFun<Base></code> object.
00630 
00631 \tparam VectorBase
00632 is a simple vector class with elements of type \a Base.
00633 
00634 \tparam VectorSet
00635 is a simple vector class with elements of type 
00636 \c bool or \c std::set<size_t>.
00637 
00638 \tparam VectorSize
00639 is a simple vector class with elements of type \c size_t.
00640 
00641 \param x [in]
00642 is a vector specifing the point at which to compute the Jacobian.
00643 
00644 \param p [in]
00645 is the sparsity pattern for the Jacobian that we are calculating.
00646 
00647 \param row [in]
00648 is the vector of row indices for the returned Jacobian values.
00649 
00650 \param col [in]
00651 is the vector of columns indices for the returned Jacobian values.
00652 It must have the same size as \c row.
00653 
00654 \param jac [out]
00655 is the vector of Jacobian values.
00656 It must have the same size as \c row. 
00657 The return value <code>jac[k]</code> is the partial of the
00658 <code>row[k]</code> range component of the function with respect
00659 the the <code>col[k]</code> domain component of its argument.
00660 
00661 \param work [in,out]
00662 this structure contains information that depends on the function object, 
00663 sparsity pattern, \c row vector, and \c col vector.
00664 If they are not the same between calls to \c SparseJacobianForward, 
00665 \c work.clear() must be called to reinitialize them.
00666 
00667 \return
00668 Is the number of first order forward sweeps used to compute the
00669 requested Jacobian values. The total work, not counting the zero order
00670 forward sweep, or the time to combine computations, is proportional to this
00671 return value.
00672 */
00673 template<class Base>
00674 template <class VectorBase, class VectorSet, class VectorSize>
00675 size_t ADFun<Base>::SparseJacobianForward(
00676      const VectorBase&     x    ,
00677      const VectorSet&      p    ,
00678      const VectorSize&     row  ,
00679      const VectorSize&     col  ,
00680      VectorBase&           jac  ,
00681      sparse_jacobian_work& work )
00682 {
00683      size_t n = Domain();
00684      size_t m = Range();
00685 # ifndef NDEBUG
00686      size_t k, K = jac.size();
00687      CPPAD_ASSERT_KNOWN(
00688           size_t(x.size()) == n ,
00689           "SparseJacobianForward: size of x not equal domain dimension for f."
00690      ); 
00691      CPPAD_ASSERT_KNOWN(
00692           size_t(row.size()) == K && size_t(col.size()) == K ,
00693           "SparseJacobianForward: either r or c does not have "
00694           "the same size as jac."
00695      ); 
00696      CPPAD_ASSERT_KNOWN(
00697           work.color.size() == 0 || work.color.size() == n,
00698           "SparseJacobianForward: invalid value in work."
00699      );
00700      for(k = 0; k < K; k++)
00701      {    CPPAD_ASSERT_KNOWN(
00702                row[k] < m,
00703                "SparseJacobianForward: invalid value in r."
00704           );
00705           CPPAD_ASSERT_KNOWN(
00706                col[k] < n,
00707                "SparseJacobianForward: invalid value in c."
00708           );
00709      }
00710      if( work.color.size() != 0 )
00711           for(size_t j = 0; j < n; j++) CPPAD_ASSERT_KNOWN(
00712                work.color[j] <= n,
00713                "SparseJacobianForward: invalid value in work."
00714      );
00715 # endif
00716 
00717      typedef typename VectorSet::value_type Set_type;
00718      typedef typename internal_sparsity<Set_type>::pattern_type Pattern_type;
00719      Pattern_type s_transpose;
00720      if( work.color.size() == 0 )
00721      {    bool transpose = true;
00722           sparsity_user2internal(s_transpose, p, m, n, transpose);
00723      }
00724      size_t n_sweep = SparseJacobianFor(x, s_transpose, row, col, jac, work);
00725      return n_sweep;
00726 }
00727 /*!
00728 Compute user specified subset of a sparse Jacobian using forward mode.
00729 
00730 The C++ source code corresponding to this operation is
00731 \verbatim
00732      SparceJacobianReverse(x, p, row, col, jac, work)
00733 \endverbatim
00734 
00735 \tparam Base
00736 is the base type for the recording that is stored in this
00737 <code>ADFun<Base></code> object.
00738 
00739 \tparam VectorBase
00740 is a simple vector class with elements of type \a Base.
00741 
00742 \tparam VectorSet
00743 is a simple vector class with elements of type 
00744 \c bool or \c std::set<size_t>.
00745 
00746 \tparam VectorSize
00747 is a simple vector class with elements of type \c size_t.
00748 
00749 \param x [in]
00750 is a vector specifing the point at which to compute the Jacobian.
00751 
00752 \param p [in]
00753 is the sparsity pattern for the Jacobian that we are calculating.
00754 
00755 \param row [in]
00756 is the vector of row indices for the returned Jacobian values.
00757 
00758 \param col [in]
00759 is the vector of columns indices for the returned Jacobian values.
00760 It must have the same size as \c row.
00761 
00762 \param jac [out]
00763 is the vector of Jacobian values.
00764 It must have the same size as \c row. 
00765 The return value <code>jac[k]</code> is the partial of the
00766 <code>row[k]</code> range component of the function with respect
00767 the the <code>col[k]</code> domain component of its argument.
00768 
00769 \param work [in,out]
00770 this structure contains information that depends on the function object, 
00771 sparsity pattern, \c row vector, and \c col vector.
00772 If they are not the same between calls to \c SparseJacobianReverse, 
00773 \c work.clear() must be called to reinitialize them.
00774 
00775 \return
00776 Is the number of first order reverse sweeps used to compute the
00777 reverse Jacobian values. The total work, not counting the zero order
00778 forward sweep, or the time to combine computations, is proportional to this
00779 return value.
00780 */
00781 template<class Base>
00782 template <class VectorBase, class VectorSet, class VectorSize>
00783 size_t ADFun<Base>::SparseJacobianReverse(
00784      const VectorBase&     x    ,
00785      const VectorSet&      p    ,
00786      const VectorSize&     row  ,
00787      const VectorSize&     col  ,
00788      VectorBase&           jac  ,
00789      sparse_jacobian_work& work )
00790 {
00791      size_t m = Range();
00792      size_t n = Domain();
00793 # ifndef NDEBUG
00794      size_t k, K = jac.size();
00795      CPPAD_ASSERT_KNOWN(
00796           size_t(x.size()) == n ,
00797           "SparseJacobianReverse: size of x not equal domain dimension for f."
00798      ); 
00799      CPPAD_ASSERT_KNOWN(
00800           size_t(row.size()) == K && size_t(col.size()) == K ,
00801           "SparseJacobianReverse: either r or c does not have "
00802           "the same size as jac."
00803      ); 
00804      CPPAD_ASSERT_KNOWN(
00805           work.color.size() == 0 || work.color.size() == m,
00806           "SparseJacobianReverse: invalid value in work."
00807      );
00808      for(k = 0; k < K; k++)
00809      {    CPPAD_ASSERT_KNOWN(
00810                row[k] < m,
00811                "SparseJacobianReverse: invalid value in r."
00812           );
00813           CPPAD_ASSERT_KNOWN(
00814                col[k] < n,
00815                "SparseJacobianReverse: invalid value in c."
00816           );
00817      }
00818      if( work.color.size() != 0 )
00819           for(size_t i = 0; i < m; i++) CPPAD_ASSERT_KNOWN(
00820                work.color[i] <= m,
00821                "SparseJacobianReverse: invalid value in work."
00822      );
00823 # endif
00824  
00825      typedef typename VectorSet::value_type Set_type;
00826      typedef typename internal_sparsity<Set_type>::pattern_type Pattern_type;
00827      Pattern_type s;
00828      if( work.color.size() == 0 )
00829      {    bool transpose = false;
00830           sparsity_user2internal(s, p, m, n, transpose);
00831      }
00832      size_t n_sweep = SparseJacobianRev(x, s, row, col, jac, work);
00833      return n_sweep;
00834 }
00835 /*!
00836 Compute a sparse Jacobian.
00837 
00838 The C++ source code corresponding to this operation is
00839 \verbatim
00840      jac = SparseJacobian(x, p)
00841 \endverbatim
00842 
00843 \tparam Base
00844 is the base type for the recording that is stored in this
00845 <code>ADFun<Base></code> object.
00846 
00847 \tparam VectorBase
00848 is a simple vector class with elements of type \a Base.
00849 
00850 \tparam VectorSet
00851 is a simple vector class with elements of type 
00852 \c bool or \c std::set<size_t>.
00853 
00854 \param x [in]
00855 is a vector specifing the point at which to compute the Jacobian.
00856 
00857 \param p [in]
00858 is the sparsity pattern for the Jacobian that we are calculating.
00859 
00860 \return
00861 Will be a vector if size \c m * n containing the Jacobian at the
00862 specified point (in row major order).
00863 */
00864 template <class Base>
00865 template <class VectorBase, class VectorSet>
00866 VectorBase ADFun<Base>::SparseJacobian(
00867      const VectorBase& x, const VectorSet& p
00868 )
00869 {    size_t i, j, k;
00870 
00871      size_t m = Range();
00872      size_t n = Domain();
00873      VectorBase jac(m * n);
00874 
00875      CPPAD_ASSERT_KNOWN(
00876           size_t(x.size()) == n,
00877           "SparseJacobian: size of x not equal domain size for f."
00878      );
00879      CheckSimpleVector<Base, VectorBase>();
00880 
00881      typedef typename VectorSet::value_type Set_type;
00882      typedef typename internal_sparsity<Set_type>::pattern_type Pattern_type;
00883 
00884      // initialize the return value as zero
00885      Base zero(0);
00886      for(i = 0; i < m; i++)
00887           for(j = 0; j < n; j++)
00888                jac[i * n + j] = zero;
00889 
00890      sparse_jacobian_work work;
00891      CppAD::vector<size_t> row;
00892      CppAD::vector<size_t> col;
00893      if( n <= m )
00894      {
00895           // need an internal copy of sparsity pattern
00896           Pattern_type s_transpose;
00897           bool transpose = true;
00898           sparsity_user2internal(s_transpose, p, m, n, transpose);
00899 
00900           k = 0;
00901           for(j = 0; j < n; j++)
00902           {    s_transpose.begin(j);
00903                i = s_transpose.next_element();
00904                while( i != s_transpose.end() )
00905                {    row.push_back(i);
00906                     col.push_back(j);
00907                     k++;
00908                     i = s_transpose.next_element();
00909                }
00910           } 
00911           size_t K = k;
00912           VectorBase J(K);
00913      
00914           // now we have folded this into the following case
00915           SparseJacobianFor(x, s_transpose, row, col, J, work);
00916 
00917           // now set the non-zero return values
00918           for(k = 0; k < K; k++)
00919                jac[ row[k] * n + col[k] ] = J[k];
00920      }
00921      else
00922      {
00923           // need an internal copy of sparsity pattern
00924           Pattern_type s;
00925           bool transpose = false;
00926           sparsity_user2internal(s, p, m, n, transpose);
00927 
00928           k = 0;
00929           for(i = 0; i < m; i++)
00930           {    s.begin(i);
00931                j = s.next_element();
00932                while( j != s.end() )
00933                {    row.push_back(i);
00934                     col.push_back(j);
00935                     k++;
00936                     j = s.next_element();
00937                }
00938           } 
00939           size_t K = k;
00940           VectorBase J(K);
00941 
00942           // now we have folded this into the following case
00943           SparseJacobianRev(x, s, row, col, J, work);
00944 
00945           // now set the non-zero return values
00946           for(k = 0; k < K; k++)
00947                jac[ row[k] * n + col[k] ] = J[k];
00948      }
00949 
00950      return jac;
00951 } 
00952 
00953 /*!
00954 Compute a sparse Jacobian.
00955 
00956 The C++ source code corresponding to this operation is
00957 \verbatim
00958      jac = SparseJacobian(x)
00959 \endverbatim
00960 
00961 \tparam Base
00962 is the base type for the recording that is stored in this
00963 <code>ADFun<Base></code> object.
00964 
00965 \tparam VectorBase
00966 is a simple vector class with elements of the \a Base.
00967 
00968 \param x [in]
00969 is a vector specifing the point at which to compute the Jacobian.
00970 
00971 \return
00972 Will be a vector of size \c m * n containing the Jacobian at the
00973 specified point (in row major order).
00974 */
00975 template <class Base>
00976 template <class VectorBase>
00977 VectorBase ADFun<Base>::SparseJacobian( const VectorBase& x )
00978 {    typedef CppAD::vectorBool   VectorBool;
00979 
00980      size_t m = Range();
00981      size_t n = Domain();
00982 
00983      // sparsity pattern for Jacobian
00984      VectorBool p(m * n);
00985 
00986      if( n <= m )
00987      {    size_t j, k;
00988 
00989           // use forward mode 
00990           VectorBool r(n * n);
00991           for(j = 0; j < n; j++)
00992           {    for(k = 0; k < n; k++)
00993                     r[j * n + k] = false;
00994                r[j * n + j] = true;
00995           }
00996           p = ForSparseJac(n, r);
00997      }
00998      else
00999      {    size_t i, k;
01000 
01001           // use reverse mode 
01002           VectorBool s(m * m);
01003           for(i = 0; i < m; i++)
01004           {    for(k = 0; k < m; k++)
01005                     s[i * m + k] = false;
01006                s[i * m + i] = true;
01007           }
01008           p = RevSparseJac(m, s);
01009      }
01010      return SparseJacobian(x, p);
01011 }
01012 
01013 } // END_CPPAD_NAMESPACE
01014 # endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines