CppAD: A C++ Algorithmic Differentiation Package  20130102
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-12 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      recomputed
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 $icode%jac% = %f%.SparseJacobian(%x%)
00035 %jac% = %f%.SparseJacobian(%x%, %p%)
00036 %n_sweep% = %f%.SparseJacobianForward(%x%, %p%, %row%, %col%, %jac%, %work%)
00037 %n_sweep% = %f%.SparseJacobianReverse(%x%, %p%, %row%, %col%, %jac%, %work%)
00038 %$$
00039 
00040 $head Purpose$$
00041 We use $latex n$$ for the $cref/domain/seq_property/Domain/$$ size,
00042 and $latex m$$ for the $cref/range/seq_property/Range/$$ size of $icode f$$.
00043 We use $latex F : \B{R}^n \rightarrow \B{R}^m$$ do denote the
00044 $cref/AD function/glossary/AD Function/$$
00045 corresponding to $icode f$$. 
00046 The syntax above sets $icode jac$$ to the Jacobian 
00047 $latex \[
00048      jac = F^{(1)} (x) 
00049 \] $$
00050 This routine takes advantage of the sparsity of the Jacobian
00051 in order to reduce the amount of computation necessary.
00052 If $icode row$$ and $icode col$$ are present, it also takes
00053 advantage of the reduced set of elements of the Jacobian that
00054 need to be computed.
00055 One can use speed tests (e.g. $cref speed_test$$)
00056 to verify that results are computed faster
00057 than when using the routine $cref Jacobian$$.
00058 
00059 $head f$$
00060 The object $icode f$$ has prototype
00061 $codei%
00062      ADFun<%Base%> %f%
00063 %$$
00064 Note that the $cref ADFun$$ object $icode f$$ is not $code const$$
00065 (see $cref/Uses Forward/sparse_jacobian/Uses Forward/$$ below).
00066 
00067 $head x$$
00068 The argument $icode x$$ has prototype
00069 $codei%
00070      const %VectorBase%& %x%
00071 %$$
00072 (see $cref/VectorBase/sparse_jacobian/VectorBase/$$ below)
00073 and its size 
00074 must be equal to $icode n$$, the dimension of the
00075 $cref/domain/seq_property/Domain/$$ space for $icode f$$.
00076 It specifies
00077 that point at which to evaluate the Jacobian.
00078 
00079 $head p$$
00080 The argument $icode p$$ is optional and has prototype
00081 $codei%
00082      const %VectorSet%& %p%
00083 %$$
00084 (see $cref/VectorSet/sparse_jacobian/VectorSet/$$ below).
00085 If it has elements of type $code bool$$,
00086 its size is $latex m * n$$.
00087 If it has elements of type $code std::set<size_t>$$,
00088 its size is $latex m$$ and all its set elements are between
00089 zero and $latex n - 1$$.
00090 It specifies a 
00091 $cref/sparsity pattern/glossary/Sparsity Pattern/$$ 
00092 for the Jacobian $latex F^{(1)} (x)$$.
00093 $pre
00094 
00095 $$
00096 If this sparsity pattern does not change between calls to 
00097 $codei SparseJacobian$$, it should be faster to calculate $icode p$$ once 
00098 (using $cref ForSparseJac$$ or $cref RevSparseJac$$)
00099 and then pass $icode p$$ to $codei SparseJacobian$$.
00100 In addition,
00101 if you specify $icode p$$, CppAD will use the same
00102 type of sparsity representation 
00103 (vectors of $code bool$$ or vectors of $code std::set<size_t>$$)
00104 for its internal calculations.
00105 Otherwise, the representation
00106 for the internal calculations is unspecified.
00107 
00108 $head row, col$$
00109 The arguments $icode row$$ and $icode col$$ are optional and have prototype
00110 $codei%
00111      const %VectorSize%& %row%
00112      const %VectorSize%& %col%
00113 %$$
00114 (see $cref/VectorSize/sparse_jacobian/VectorSize/$$ below).
00115 They specify which rows and columns of $latex F^{(1)} (x)$$ are
00116 returned and in what order.
00117 We use $latex K$$ to denote the value $icode%jac%.size()%$$
00118 which must also equal the size of $icode row$$ and $icode col$$.
00119 Furthermore,
00120 for $latex k = 0 , \ldots , K-1$$, it must hold that
00121 $latex row[k] < m$$ and $latex col[k] < n$$.
00122 In addition, 
00123 all of the $latex (row[k], col[k])$$ pairs must correspond to a true value
00124 in the sparsity pattern $icode p$$.
00125 
00126 $head jac$$
00127 The result $icode jac$$ has prototype
00128 $codei%
00129      %VectorBase%& %jac%
00130 %$$
00131 In the case where the arguments $icode row$$ and $icode col$$ are not present,
00132 the size of $icode jac$$ is $latex m * n$$ and
00133 for $latex i = 0 , \ldots , m-1$$,
00134 $latex j = 0 , \ldots , n-1$$,
00135 $latex \[
00136      jac [ i * n + j ] = \D{ F_i }{ x_j } (x)
00137 \] $$
00138 $pre
00139 
00140 $$
00141 In the case where the arguments $icode row$$ and $icode col$$ are present,
00142 we use $latex K$$ to denote the size of $icode jac$$. 
00143 The input value of its elements does not matter.
00144 Upon return, for $latex k = 0 , \ldots , K - 1$$,
00145 $latex \[
00146      jac [ k ] = \D{ F_i }{ x_j } (x)
00147      \; , \;
00148      \; {\rm where} \;
00149      i = row[k]
00150      \; {\rm and } \;
00151      j = col[k]
00152 \] $$
00153 
00154 $head work$$
00155 If this argument is present, it has prototype
00156 $icode%
00157      sparse_jacobian_work& %work%
00158 %$$
00159 This object can only be used with the routines 
00160 $code SparseJacobianForward$$ and $code SparseJacobianReverse$$.
00161 During its the first use, information is stored in $icode work$$.
00162 This is used to reduce the work done by future calls to the same mode
00163 (forward or reverse), 
00164 the same $icode f$$, $icode p$$, $icode row$$, and $icode col$$.
00165 If a future call is for a different mode,
00166 or any of these values have changed,
00167 you must first call $icode%work%.clear()%$$
00168 to inform CppAD that this information needs to be recomputed.
00169 
00170 $head n_sweep$$
00171 The return value $icode n_sweep$$ has prototype
00172 $codei%
00173      size_t %n_sweep%
00174 %$$
00175 If $code SparseJacobianForward$$ ($code SparseJacobianReverse$$) is used, 
00176 $icode n_sweep$$ is the number of first order forward (reverse) sweeps 
00177 used to compute the requested Jacobian values. 
00178 This is proportional to the total work that $code SparseJacobian$$ does, 
00179 not counting the zero order forward sweep, 
00180 or the work to combine multiple columns (rows) into a single sweep.
00181 
00182 $head VectorBase$$
00183 The type $icode VectorBase$$ must be a $cref SimpleVector$$ class with
00184 $cref/elements of type/SimpleVector/Elements of Specified Type/$$
00185 $icode Base$$.
00186 The routine $cref CheckSimpleVector$$ will generate an error message
00187 if this is not the case.
00188 
00189 $head VectorSet$$
00190 The type $icode VectorSet$$ must be a $cref SimpleVector$$ class with
00191 $cref/elements of type/SimpleVector/Elements of Specified Type/$$
00192 $code bool$$ or $code std::set<size_t>$$;
00193 see $cref/sparsity pattern/glossary/Sparsity Pattern/$$ for a discussion
00194 of the difference.
00195 The routine $cref CheckSimpleVector$$ will generate an error message
00196 if this is not the case.
00197 
00198 $subhead Restrictions$$
00199 If $icode VectorSet$$ has elements of $code std::set<size_t>$$,
00200 then $icode%p%[%i%]%$$ must return a reference (not a copy) to the 
00201 corresponding set.
00202 According to section 26.3.2.3 of the 1998 C++ standard,
00203 $code std::valarray< std::set<size_t> >$$ does not satisfy
00204 this condition. 
00205 
00206 $head VectorSize$$
00207 The type $icode VectorSize$$ must be a $cref SimpleVector$$ class with
00208 $cref/elements of type/SimpleVector/Elements of Specified Type/$$
00209 $code size_t$$.
00210 The routine $cref CheckSimpleVector$$ will generate an error message
00211 if this is not the case.
00212 
00213 $head Uses Forward$$
00214 After each call to $cref Forward$$,
00215 the object $icode f$$ contains the corresponding 
00216 $cref/Taylor coefficients/glossary/Taylor Coefficient/$$.
00217 After a call to any of the sparse Jacobian routines,
00218 the zero order Taylor coefficients correspond to
00219 $icode%f%.Forward(0, %x%)%$$
00220 and the other coefficients are unspecified.
00221 
00222 After $code SparseJacobian$$,
00223 the previous calls to $cref Forward$$ are undefined.
00224 
00225 $head Example$$
00226 $children%
00227      example/sparse_jacobian.cpp
00228 %$$
00229 The routine
00230 $cref sparse_jacobian.cpp$$
00231 is examples and tests of $code sparse_jacobian$$.
00232 It return $code true$$, if it succeeds and $code false$$ otherwise.
00233 
00234 $end
00235 ==============================================================================
00236 */
00237 # include <cppad/local/std_set.hpp>
00238 
00239 CPPAD_BEGIN_NAMESPACE
00240 /*!
00241 \defgroup sparse_jacobian_hpp sparse_jacobian.hpp
00242 \{
00243 \file sparse_jacobian.hpp
00244 Sparse Jacobian driver routine and helper functions.
00245 */
00246 // ===========================================================================
00247 /*!
00248 class used by SparseJacobian to hold information so it does not need to be 
00249 recomputed.
00250 */
00251 class sparse_jacobian_work {
00252      public:
00253           /// version of user row array with the extra value m at end
00254           CppAD::vector<size_t> user_row;
00255           /// version of user col array with the extra value n at end
00256           CppAD::vector<size_t> user_col;
00257           /// indices that sort the user arrays by row 
00258           /// with the extra value K at the end
00259           CppAD::vector<size_t> sort_row;
00260           /// indices that sort the user arrays by column 
00261           /// with the extra value K at the end
00262           CppAD::vector<size_t> sort_col;
00263           /// results of the coloring algorithm
00264           CppAD::vector<size_t> color;
00265           /// inform CppAD that this information needs to be recomputed
00266           void clear(void)
00267           {    user_row.clear();
00268                user_col.clear();
00269                sort_row.clear();
00270                sort_col.clear();
00271                color.clear();
00272           }
00273 };
00274 // ===========================================================================
00275 /*!
00276 Private helper function forward mode cases
00277 
00278 \tparam Base
00279 See \c SparseJacobianForward(x, p, row, col, jac, work).
00280 
00281 \tparam VectorBase
00282 See \c SparseJacobianForward(x, p, row, col, jac, work).
00283 
00284 \tparam VectorSet
00285 is either \c sparse_pack, \c sparse_set or \c sparse_list.
00286 
00287 \param x
00288 See \c SparseJacobianForward(x, p, row, col, jac, work).
00289 
00290 \param p_transpose
00291 If <code>work.color.size() != 0</code>, 
00292 then \c p_transpose is not used.
00293 Otherwise, it is a
00294 sparsity pattern for the transpose of the Jacobian of this ADFun<Base> object.
00295 Note that we do not change the values in \c p_transpose,
00296 but is not \c const because we use its iterator facility.
00297 
00298 
00299 \param jac
00300 See \c SparseJacobianForward(x, p, row, col, jac, work).
00301 
00302 \param work
00303 See \c SparseJacobianForward(x, p, row, col, jac, work).
00304 In addition, \c work.color contains information that is computed by
00305 \c SparseJacobainFor.
00306 If the sparsity pattern, \c row vector, and \c col vectors
00307 are not the same between calls, \c work.clear() must be called
00308 to reinitialize \c work.color.
00309 
00310 \return
00311 See \c SparseJacobianForward(x, p, row, col, jac, work).
00312 */
00313 template<class Base>
00314 template <class VectorBase, class VectorSet>
00315 size_t ADFun<Base>::SparseJacobianFor(
00316      const VectorBase&     x           ,
00317      VectorSet&            p_transpose ,
00318      VectorBase&           jac         ,
00319      sparse_jacobian_work& work        )
00320 {
00321      using   CppAD::vectorBool;
00322      size_t i, j, k, ell;
00323 
00324      CppAD::vector<size_t>& row(work.user_row);
00325      CppAD::vector<size_t>& col(work.user_col);
00326      CppAD::vector<size_t>& sort_col(work.sort_col);
00327      CppAD::vector<size_t>& color(work.color);
00328 
00329      size_t m = Range();
00330      size_t n = Domain();
00331 
00332      // some values
00333      const Base zero(0);
00334      const Base one(1);
00335 
00336      // check VectorBase is Simple Vector class with Base type elements
00337      CheckSimpleVector<Base, VectorBase>();
00338 
00339      CPPAD_ASSERT_UNKNOWN( size_t(x.size()) == n );
00340      CPPAD_ASSERT_UNKNOWN( color.size() == 0 || color.size() == n );
00341 
00342      // number of components of Jacobian that are required
00343      size_t K = size_t(jac.size());
00344      CPPAD_ASSERT_UNKNOWN( row.size() == K+1 );
00345      CPPAD_ASSERT_UNKNOWN( col.size() == K+1 );
00346      CPPAD_ASSERT_UNKNOWN( row[K] == m );
00347      CPPAD_ASSERT_UNKNOWN( col[K] == n );
00348 
00349      // Point at which we are evaluating the Jacobian
00350      Forward(0, x);
00351 
00352      if( color.size() == 0 )
00353      {
00354           CPPAD_ASSERT_UNKNOWN( p_transpose.n_set() ==  n );
00355           CPPAD_ASSERT_UNKNOWN( p_transpose.end() ==  m );
00356 
00357           // rows and columns that are in the returned jacobian
00358           VectorSet r_used, c_used;
00359           r_used.resize(n, m);
00360           c_used.resize(m, n);
00361           k = 0;
00362           while( k < K )
00363           {    CPPAD_ASSERT_UNKNOWN( 
00364                     row[sort_col[k]] < m && col[sort_col[k]] < n 
00365                );
00366                CPPAD_ASSERT_UNKNOWN( 
00367                     k == 0 || col[sort_col[k-1]] <= col[sort_col[k]] 
00368                );
00369                CPPAD_ASSERT_KNOWN(
00370                     p_transpose.is_element(col[sort_col[k]], row[sort_col[k]]) ,
00371                     "SparseJacobianForward: "
00372                     "an (row, col) pair is not in sparsity pattern."
00373                );
00374                r_used.add_element(col[sort_col[k]], row[sort_col[k]]);
00375                c_used.add_element(row[sort_col[k]], col[sort_col[k]]);
00376                k++;
00377           }
00378      
00379           // given a row index, which columns are non-zero and not used
00380           VectorSet not_used;
00381           not_used.resize(m, n);
00382           for(j = 0; j < n; j++)
00383           {    p_transpose.begin(j);
00384                i = p_transpose.next_element();
00385                while( i != p_transpose.end() )
00386                {    if( ! c_used.is_element(i, j) )
00387                          not_used.add_element(i, j);
00388                     i = p_transpose.next_element();
00389                }
00390           }
00391 
00392           // initial coloring
00393           color.resize(n);
00394           for(j = 0; j < n; j++)
00395                color[j] = j;
00396      
00397           // See GreedyPartialD2Coloring Algorithm Section 3.6.2 of
00398           // Graph Coloring in Optimization Revisited by
00399           // Assefaw Gebremedhin, Fredrik Maane, Alex Pothen
00400           vectorBool forbidden(n);
00401           for(j = 1; j < n; j++)
00402           {
00403                // initial all colors as ok for this column
00404                // (value of forbidden for ell > j does not matter)
00405                for(ell = 0; ell <= j; ell++)
00406                     forbidden[ell] = false;
00407      
00408                // for each row that is non-zero for this column
00409                p_transpose.begin(j);
00410                i = p_transpose.next_element();
00411                while( i != p_transpose.end() )
00412                {    // for each column that this row uses
00413                     c_used.begin(i);
00414                     ell = c_used.next_element();
00415                     while( ell != c_used.end() )
00416                     {    // if this is not the same column, forbid its color
00417                          if( ell < j )
00418                               forbidden[ color[ell] ] = true;
00419                          ell = c_used.next_element();
00420                     }
00421                     i = p_transpose.next_element();
00422                }
00423      
00424                // for each row that this column uses
00425                r_used.begin(j);
00426                i = r_used.next_element();
00427                while( i != r_used.end() )
00428                {    // For each column that is non-zero for this row
00429                     // (the used columns have already been checked above).
00430                     not_used.begin(i);
00431                     ell = not_used.next_element();
00432                     while( ell != not_used.end() )
00433                     {    // if this is not the same column, forbid its color
00434                          if( ell < j )
00435                               forbidden[ color[ell] ] = true;
00436                          ell = not_used.next_element();
00437                     }
00438                     i = r_used.next_element();
00439                }
00440 
00441                // pick the color with smallest index
00442                ell = 0;
00443                while( forbidden[ell] )
00444                {    ell++;
00445                     CPPAD_ASSERT_UNKNOWN( ell <= j );
00446                }
00447                color[j] = ell;
00448           }
00449      }
00450      size_t n_color = 1;
00451      for(ell = 0; ell < n; ell++) 
00452           n_color = std::max(n_color, color[ell] + 1);
00453 
00454      // direction vector for calls to forward
00455      VectorBase dx(n);
00456 
00457      // location for return values from forward
00458      VectorBase dy(m);
00459 
00460      // initialize the return value
00461      for(k = 0; k < K; k++)
00462           jac[k] = zero;
00463 
00464      // loop over colors
00465      size_t n_sweep = 0;
00466      for(ell = 0; ell < n_color; ell++)
00467      {    bool any = false;
00468           k = 0;
00469           for(j = 0; j < n; j++) if( color[j] == ell )
00470           {    // find first k such that col[sort_col[k]] has color ell
00471                if( ! any )
00472                {    while( col[sort_col[k]] < j )
00473                          k++;
00474                     any = col[sort_col[k]] == j;
00475                }
00476           }
00477           if( any )
00478           {    n_sweep++;
00479                // combine all columns with this color
00480                for(j = 0; j < n; j++)
00481                {    dx[j] = zero;
00482                     if( color[j] == ell )
00483                          dx[j] = one;
00484                }
00485                // call forward mode for all these columns at once
00486                dy = Forward(1, dx);
00487 
00488                // set the corresponding components of the result
00489                for(j = 0; j < n; j++) if( color[j] == ell )
00490                {    // find first index in c for this column
00491                     while( col[sort_col[k]] < j )
00492                          k++;
00493                     // extract the row results for this column
00494                     while( col[sort_col[k]] == j ) 
00495                     {    jac[ sort_col[k] ] = dy[row[sort_col[k]]];
00496                          k++;
00497                     }
00498                }
00499           }
00500      }
00501      return n_sweep;
00502 }
00503 /*!
00504 Private helper function for reverse mode cases.
00505 
00506 \tparam Base
00507 See \c SparseJacobianForward(x, p, row, col, jac, work).
00508 
00509 \tparam VectorBase
00510 See \c SparseJacobianForward(x, p, row, col, jac, work).
00511 
00512 \tparam VectorSet
00513 is either \c sparse_pack, \c sparse_set or \c sparse_list.
00514 
00515 \param x
00516 See \c SparseJacobianForward(x, p, row, col, jac, work).
00517 
00518 \param p
00519 If <code>work.color.size() != 0</code>, then \c p is not used.
00520 Otherwise, it is a
00521 sparsity pattern for the Jacobian of this ADFun<Base> object.
00522 Note that we do not change the values in \c p_transpose,
00523 but is not \c const because we use its iterator facility.
00524 
00525 \param jac
00526 See \c SparseJacobianForward(x, p, row, col, jac, work).
00527 
00528 \param work
00529 See \c SparseJacobianReverse(x, p, row, col, jac, work).
00530 In addition, \c work.color contains information that is computed by
00531 \c SparseJacobainRev.
00532 If the sparsity pattern, \c row vector, and \c col vectors
00533 are not the same between calls, \c work.clear() must be called
00534 to reinitialize \c work.color.
00535 
00536 \return
00537 See \c SparseJacobianReverse(x, p, row, col, jac, work).
00538 */
00539 template<class Base>
00540 template <class VectorBase, class VectorSet>
00541 size_t ADFun<Base>::SparseJacobianRev(
00542      const VectorBase&     x           ,
00543      VectorSet&            p           ,
00544      VectorBase&           jac         ,
00545      sparse_jacobian_work& work        )
00546 {
00547      using   CppAD::vectorBool;
00548      size_t i, j, k, ell;
00549 
00550      CppAD::vector<size_t>& row(work.user_row);
00551      CppAD::vector<size_t>& col(work.user_col);
00552      CppAD::vector<size_t>& sort_row(work.sort_row);
00553      CppAD::vector<size_t>& color(work.color);
00554 
00555      size_t m = Range();
00556      size_t n = Domain();
00557 
00558      // some values
00559      const Base zero(0);
00560      const Base one(1);
00561 
00562      // check VectorBase is Simple Vector class with Base type elements
00563      CheckSimpleVector<Base, VectorBase>();
00564 
00565      CPPAD_ASSERT_UNKNOWN( size_t(x.size()) == n );
00566      CPPAD_ASSERT_UNKNOWN (color.size() == m || color.size() == 0 );
00567 
00568      // number of components of Jacobian that are required
00569      size_t K = size_t(jac.size());
00570      CPPAD_ASSERT_UNKNOWN( row.size() == K+1 );
00571      CPPAD_ASSERT_UNKNOWN( col.size() == K+1 );
00572      CPPAD_ASSERT_UNKNOWN( row[K] == m );
00573      CPPAD_ASSERT_UNKNOWN( col[K] == n );
00574 
00575      // Point at which we are evaluating the Jacobian
00576      Forward(0, x);
00577 
00578      if( color.size() == 0 )
00579      {    CPPAD_ASSERT_UNKNOWN( p.n_set() ==  m );
00580           CPPAD_ASSERT_UNKNOWN( p.end() ==  n );
00581 
00582           // rows and columns that are in the returned jacobian
00583           VectorSet r_used, c_used;
00584           r_used.resize(n, m);
00585           c_used.resize(m, n);
00586           k = 0;
00587           while( k < K )
00588           {    CPPAD_ASSERT_UNKNOWN( 
00589                     row[sort_row[k]] < m && col[sort_row[k]] < n 
00590                );
00591                CPPAD_ASSERT_UNKNOWN( 
00592                     k == 0 || row[sort_row[k-1]] <= row[sort_row[k]] 
00593                );
00594                CPPAD_ASSERT_KNOWN(
00595                     p.is_element(row[sort_row[k]], col[sort_row[k]]) ,
00596                     "SparseJacobianReverse: "
00597                     "an (row, col) pair is not in sparsity pattern."
00598                );
00599                r_used.add_element(col[sort_row[k]], row[sort_row[k]]);
00600                c_used.add_element(row[sort_row[k]], col[sort_row[k]]);
00601                k++;
00602           }
00603      
00604           // given a column index, which rows are non-zero and not used
00605           VectorSet not_used;
00606           not_used.resize(n, m);
00607           for(i = 0; i < m; i++)
00608           {    p.begin(i);
00609                j = p.next_element();
00610                while( j != p.end() )
00611                {    if( ! r_used.is_element(j , i) )
00612                          not_used.add_element(j, i);
00613                     j = p.next_element();
00614                }
00615           }
00616      
00617           // initial coloring
00618           color.resize(m);
00619           for(i = 0; i < m; i++)
00620                color[i] = i;
00621      
00622           // See GreedyPartialD2Coloring Algorithm Section 3.6.2 of
00623           // Graph Coloring in Optimization Revisited by
00624           // Assefaw Gebremedhin, Fredrik Maane, Alex Pothen
00625           vectorBool forbidden(m);
00626           for(i = 1; i < m; i++)
00627           {
00628                // initial all colors as ok for this row
00629                // (value of forbidden for ell > i does not matter)
00630                for(ell = 0; ell <= i; ell++)
00631                     forbidden[ell] = false;
00632      
00633                // -----------------------------------------------------
00634                // Forbid colors for which this row would destroy results 
00635                // for each column that is non-zero for this row
00636                p.begin(i);
00637                j = p.next_element();
00638                while( j != p.end() )
00639                {    // for each row that this column uses
00640                     r_used.begin(j);
00641                     ell = r_used.next_element();
00642                     while( ell != r_used.end() )
00643                     {    // if this is not the same row, forbid its color 
00644                          if( ell < i )
00645                               forbidden[ color[ell] ] = true;
00646                          ell = r_used.next_element();
00647                     }
00648                     j = p.next_element();
00649                }
00650 
00651      
00652                // -----------------------------------------------------
00653                // Forbid colors that would destroy results for this row.
00654                // for each column that this row uses
00655                c_used.begin(i);
00656                j = c_used.next_element();
00657                while( j != c_used.end() )
00658                {    // For each row that is non-zero for this column
00659                     // (the used rows have already been checked above).
00660                     not_used.begin(j);
00661                     ell = not_used.next_element();
00662                     while( ell != not_used.end() )
00663                     {    // if this is not the same row, forbid its color 
00664                          if( ell < i )
00665                               forbidden[ color[ell] ] = true;
00666                          ell = not_used.next_element();
00667                     }
00668                     j = c_used.next_element();
00669                }
00670 
00671                // pick the color with smallest index
00672                ell = 0;
00673                while( forbidden[ell] )
00674                {    ell++;
00675                     CPPAD_ASSERT_UNKNOWN( ell <= i );
00676                }
00677                color[i] = ell;
00678           }
00679      }
00680      size_t n_color = 1;
00681      for(ell = 0; ell < m; ell++) 
00682           n_color = std::max(n_color, color[ell] + 1);
00683 
00684      // weighting vector for calls to reverse
00685      VectorBase w(m);
00686 
00687      // location for return values from Reverse
00688      VectorBase dw(n);
00689 
00690      // initialize the return value
00691      for(k = 0; k < K; k++)
00692           jac[k] = zero;
00693 
00694      // loop over colors
00695      size_t n_sweep = 0;
00696      for(ell = 0; ell < n_color; ell++)
00697      {    bool any = false;
00698           k = 0;
00699           for(i = 0; i < m; i++) if( color[i] == ell )
00700           {    // find first k such that row[sort_row[k]] has color ell
00701                if( ! any )
00702                {    while( row[sort_row[k]] < i )
00703                          k++;
00704                     any = row[sort_row[k]] == i;
00705                }
00706           }
00707           if( any )
00708           {    n_sweep++;
00709                // combine all the rows with this color
00710                for(i = 0; i < m; i++)
00711                {    w[i] = zero;
00712                     if( color[i] == ell )
00713                          w[i] = one;
00714                }
00715                // call reverse mode for all these rows at once
00716                dw = Reverse(1, w);
00717 
00718                // set the corresponding components of the result
00719                for(i = 0; i < m; i++) if( color[i] == ell )
00720                {    // find first index in r for this row
00721                     while( row[sort_row[k]] < i )
00722                          k++;
00723                     // extract the row results for this row
00724                     while( row[sort_row[k]] == i ) 
00725                     {    jac[ sort_row[k] ] = dw[col[sort_row[k]]];
00726                          k++;
00727                     }
00728                }
00729           }
00730      }
00731      return n_sweep;
00732 }
00733 // ==========================================================================
00734 // Public Member functions
00735 // ==========================================================================
00736 /*!
00737 Compute user specified subset of a sparse Jacobian using forward mode.
00738 
00739 The C++ source code corresponding to this operation is
00740 \verbatim
00741      SparceJacobianForward(x, p, row, col, jac, work)
00742 \endverbatim
00743 
00744 \tparam Base
00745 is the base type for the recording that is stored in this
00746 ADFun<Base object.
00747 
00748 \tparam VectorBase
00749 is a simple vector class with elements of type \a Base.
00750 
00751 \tparam VectorSet
00752 is a simple vector class with elements of type 
00753 \c bool or \c std::set<size_t>.
00754 
00755 \tparam VectorSize
00756 is a simple vector class with elements of type \c size_t.
00757 
00758 \param x
00759 is a vector specifing the point at which to compute the Jacobian.
00760 
00761 \param p
00762 is the sparsity pattern for the Jacobian that we are calculating.
00763 
00764 \param row
00765 is the vector of row indices for the returned Jacobian values.
00766 
00767 \param col
00768 is the vector of columns indices for the returned Jacobian values.
00769 It must have the same size are r.
00770 
00771 \param jac
00772 is the vector of Jacobian values.
00773 It must have the same size are r. 
00774 The return value <code>jac[k]</code> is the partial of the
00775 <code>row[k]</code> component of the function with respect
00776 the the <code>col[k]</code> of its argument.
00777 
00778 \param work
00779 \c work contains information that depends on the function object, 
00780 sparsity pattern, \c row vector, and \c col vector.
00781 These \c work values are computed by \c SparseJacobianForward.
00782 If they are not the same between calls, \c work.clear() must be called
00783 to reinitialize them.
00784 To be more specific:
00785 \c user_row is a copy of \c row with the value \c m at the end
00786 (\c m is the range dimension),
00787 \c user_col is a copy of \c col with the value \c n at the end
00788 (\c n is the domain dimension),
00789 \c sort_col is the index vector that sorts the \c user_col
00790 vector with the extra value \c K at the end 
00791 (\c K is the size of \c jac).
00792 
00793 \return
00794 Is the number of first order forward sweeps used to compute the
00795 requested Jacobian values. The total work, not counting the zero order
00796 forward sweep, or the time to combine computations, is proportional to this
00797 return value.
00798 */
00799 template<class Base>
00800 template <class VectorBase, class VectorSet, class VectorSize>
00801 size_t ADFun<Base>::SparseJacobianForward(
00802      const VectorBase&     x    ,
00803      const VectorSet&      p    ,
00804      const VectorSize&     row  ,
00805      const VectorSize&     col  ,
00806      VectorBase&           jac  ,
00807      sparse_jacobian_work& work )
00808 {
00809      size_t n = Domain();
00810      size_t m = Range();
00811      size_t k, K = jac.size();
00812      if( work.user_row.size() == 0 )
00813      {    // create version of (row, col, k) sorted by column value
00814           work.user_col.resize(K+1);
00815           work.user_row.resize(K+1);
00816           work.sort_col.resize(K+1);
00817 
00818           // put sorted indices in user_row and user_col
00819           for(k = 0; k < K; k++)
00820           {    work.user_row[k] = row[k];
00821                work.user_col[k] = col[k];
00822           }
00823           work.user_row[K] = m;
00824           work.user_col[K] = n;
00825 
00826           // put sorting indices in sort_col
00827           index_sort(work.user_col, work.sort_col);
00828      }
00829 # ifndef NDEBUG
00830      CPPAD_ASSERT_KNOWN(
00831           size_t(x.size()) == n ,
00832           "SparseJacobianForward: size of x not equal domain dimension for f."
00833      ); 
00834      CPPAD_ASSERT_KNOWN(
00835           size_t(row.size()) == K && size_t(col.size()) == K ,
00836           "SparseJacobianForward: either r or c does not have "
00837           "the same size as jac."
00838      ); 
00839      CPPAD_ASSERT_KNOWN(
00840           work.user_row.size() == K+1 &&
00841           work.user_col.size() == K+1 &&
00842           work.sort_col.size() == K+1 ,
00843           "SparseJacobianForward: invalid value in work."
00844      );
00845      CPPAD_ASSERT_KNOWN(
00846           work.color.size() == 0 || work.color.size() == n,
00847           "SparseJacobianForward: invalid value in work."
00848      );
00849      for(k = 0; k < K; k++)
00850      {    CPPAD_ASSERT_KNOWN(
00851                row[k] < m,
00852                "SparseJacobianForward: invalid value in r."
00853           );
00854           CPPAD_ASSERT_KNOWN(
00855                col[k] < n,
00856                "SparseJacobianForward: invalid value in c."
00857           );
00858           CPPAD_ASSERT_KNOWN(
00859                work.sort_col[k] < K,
00860                "SparseJacobianForward: invalid value in work."
00861           );
00862           CPPAD_ASSERT_KNOWN(
00863                work.user_row[k] == row[k],
00864                "SparseJacobianForward: invalid value in work."
00865           );
00866           CPPAD_ASSERT_KNOWN(
00867                work.user_col[k] == col[k],
00868                "SparseJacobianForward: invalid value in work."
00869           );
00870      }
00871      if( work.color.size() != 0 )
00872           for(size_t j = 0; j < n; j++) CPPAD_ASSERT_KNOWN(
00873                work.color[j] < n,
00874                "SparseJacobianForward: invalid value in work."
00875      );
00876 # endif
00877 
00878      typedef typename VectorSet::value_type Set_type;
00879      typedef typename internal_sparsity<Set_type>::pattern_type Pattern_type;
00880      Pattern_type s_transpose;
00881      if( work.color.size() == 0 )
00882      {    bool transpose = true;
00883           sparsity_user2internal(s_transpose, p, m, n, transpose);
00884      }
00885      size_t n_sweep = SparseJacobianFor(x, s_transpose, jac, work);
00886      return n_sweep;
00887 }
00888 /*!
00889 Compute user specified subset of a sparse Jacobian using forward mode.
00890 
00891 The C++ source code corresponding to this operation is
00892 \verbatim
00893      SparceJacobianReverse(x, p, row, col, jac, work)
00894 \endverbatim
00895 
00896 \tparam Base
00897 See \c SparseJacobianForward(x, p, row, col, jac, work).
00898 
00899 \tparam VectorBase
00900 See \c SparseJacobianForward(x, p, row, col, jac, work).
00901 
00902 \tparam VectorSet
00903 See \c SparseJacobianForward(x, p, row, col, jac, work).
00904 
00905 \tparam VectorSize
00906 See \c SparseJacobianForward(x, p, row, col, jac, work).
00907 
00908 \param x
00909 See \c SparseJacobianForward(x, p, row, col, jac, work).
00910 
00911 \param p
00912 See \c SparseJacobianForward(x, p, row, col, jac, work).
00913 
00914 \param row
00915 See \c SparseJacobianForward(x, p, row, col, jac, work).
00916 
00917 \param col
00918 See \c SparseJacobianForward(x, p, row, col, jac, work).
00919 
00920 \param jac
00921 See \c SparseJacobianForward(x, p, row, col, jac, work).
00922 
00923 \param work
00924 \c work contains information that depends on the function object, 
00925 sparsity pattern, \c row vector, and \c col vector.
00926 These \c work values are computed by \c SparseJacobianReverse.
00927 If they are not the same between calls, \c work.clear() must be called
00928 to reinitialize them.
00929 To be more specific:
00930 \c user_row is a copy of \c row with the value \c m at the end
00931 (\c m is the range dimension),
00932 \c user_col is a copy of \c col with the value \c n at the end
00933 (\c n is the domain dimension),
00934 \c sort_row is the index vector that sorts the \c user_row
00935 vector with the extra value \c K at the end 
00936 (\c K is the size of \c jac).
00937 
00938 \return
00939 Is the number of first order reverse sweeps used to compute the
00940 reverse Jacobian values. The total work, not counting the zero order
00941 forward sweep, or the time to combine computations, is proportional to this
00942 return value.
00943 */
00944 template<class Base>
00945 template <class VectorBase, class VectorSet, class VectorSize>
00946 size_t ADFun<Base>::SparseJacobianReverse(
00947      const VectorBase&     x    ,
00948      const VectorSet&      p    ,
00949      const VectorSize&     row  ,
00950      const VectorSize&     col  ,
00951      VectorBase&           jac  ,
00952      sparse_jacobian_work& work )
00953 {
00954      size_t m = Range();
00955      size_t n = Domain();
00956      size_t k, K = jac.size();
00957      if( work.user_row.size() == 0 )
00958      {    // create version of (row, col, k) sorted by row row value
00959           work.user_col.resize(K+1);
00960           work.user_row.resize(K+1);
00961           work.sort_row.resize(K+1);
00962 
00963           // put sorted indices in user_row and user_col
00964           for(k = 0; k < K; k++)
00965           {    work.user_row[k] = row[k];
00966                work.user_col[k] = col[k];
00967           }
00968           work.user_row[K] = m;
00969           work.user_col[K] = n;
00970 
00971           // put sorting indices in sort_row
00972           index_sort(work.user_row, work.sort_row);
00973      }
00974 # ifndef NDEBUG
00975      CPPAD_ASSERT_KNOWN(
00976           size_t(x.size()) == n ,
00977           "SparseJacobianReverse: size of x not equal domain dimension for f."
00978      ); 
00979      CPPAD_ASSERT_KNOWN(
00980           size_t(row.size()) == K && size_t(col.size()) == K ,
00981           "SparseJacobianReverse: either r or c does not have "
00982           "the same size as jac."
00983      ); 
00984      CPPAD_ASSERT_KNOWN(
00985           work.user_row.size() == K+1 &&
00986           work.user_col.size() == K+1 &&
00987           work.sort_row.size() == K+1 ,
00988           "SparseJacobianReverse: invalid value in work."
00989      );
00990      CPPAD_ASSERT_KNOWN(
00991           work.color.size() == 0 || work.color.size() == m,
00992           "SparseJacobianReverse: invalid value in work."
00993      );
00994      for(k = 0; k < K; k++)
00995      {    CPPAD_ASSERT_KNOWN(
00996                row[k] < m,
00997                "SparseJacobianReverse: invalid value in r."
00998           );
00999           CPPAD_ASSERT_KNOWN(
01000                col[k] < n,
01001                "SparseJacobianReverse: invalid value in c."
01002           );
01003           CPPAD_ASSERT_KNOWN(
01004                work.sort_row[k] < K,
01005                "SparseJacobianReverse: invalid value in work."
01006           );
01007           CPPAD_ASSERT_KNOWN(
01008                work.user_row[k] == row[k],
01009                "SparseJacobianReverse: invalid value in work."
01010           );
01011           CPPAD_ASSERT_KNOWN(
01012                work.user_col[k] == col[k],
01013                "SparseJacobianReverse: invalid value in work."
01014           );
01015      }
01016      if( work.color.size() != 0 )
01017           for(size_t i = 0; i < m; i++) CPPAD_ASSERT_KNOWN(
01018                work.color[i] < m,
01019                "SparseJacobianReverse: invalid value in work."
01020      );
01021 # endif
01022  
01023      typedef typename VectorSet::value_type Set_type;
01024      typedef typename internal_sparsity<Set_type>::pattern_type Pattern_type;
01025      Pattern_type s;
01026      if( work.color.size() == 0 )
01027      {    bool transpose = false;
01028           sparsity_user2internal(s, p, m, n, transpose);
01029      }
01030      size_t n_sweep = SparseJacobianRev(x, s, jac, work);
01031      return n_sweep;
01032 }
01033 /*!
01034 Compute a sparse Jacobian.
01035 
01036 The C++ source code corresponding to this operation is
01037 \verbatim
01038      jac = SparseJacobian(x, p)
01039 \endverbatim
01040 
01041 \tparam Base
01042 is the base type for the recording that is stored in this
01043 ADFun<Base object.
01044 
01045 \tparam VectorBase
01046 is a simple vector class with elements of type \a Base.
01047 
01048 \tparam VectorSet
01049 is a simple vector class with elements of type 
01050 \c bool or \c std::set<size_t>.
01051 
01052 \param x
01053 is a vector specifing the point at which to compute the Jacobian.
01054 
01055 \param p
01056 is the sparsity pattern for the Jacobian that we are calculating.
01057 
01058 \return
01059 Will be a vector if size \c m * n containing the Jacobian at the
01060 specified point (in row major order).
01061 */
01062 template <class Base>
01063 template <class VectorBase, class VectorSet>
01064 VectorBase ADFun<Base>::SparseJacobian(
01065      const VectorBase& x, const VectorSet& p
01066 )
01067 {    size_t i, j, k;
01068 
01069      size_t m = Range();
01070      size_t n = Domain();
01071      VectorBase jac(m * n);
01072 
01073      CPPAD_ASSERT_KNOWN(
01074           size_t(x.size()) == n,
01075           "SparseJacobian: size of x not equal domain size for f."
01076      );
01077      CheckSimpleVector<Base, VectorBase>();
01078 
01079      typedef typename VectorSet::value_type Set_type;
01080      typedef typename internal_sparsity<Set_type>::pattern_type Pattern_type;
01081 
01082      // initialize the return value as zero
01083      Base zero(0);
01084      for(i = 0; i < m; i++)
01085           for(j = 0; j < n; j++)
01086                jac[i * n + j] = zero;
01087 
01088      sparse_jacobian_work work;
01089      if( n <= m )
01090      {    CppAD::vector<size_t>& row(work.user_row);
01091           CppAD::vector<size_t>& col(work.user_col);
01092           CppAD::vector<size_t>& sort_col(work.sort_col);
01093 
01094           // forward mode, columns are sorted
01095           assert( row.size() == 0 );
01096           assert( col.size() == 0 );
01097           assert( sort_col.size() == 0 );
01098 
01099           // need an internal copy of sparsity pattern
01100           Pattern_type s_transpose;
01101           bool transpose = true;
01102           sparsity_user2internal(s_transpose, p, m, n, transpose);
01103 
01104           k = 0;
01105           for(j = 0; j < n; j++)
01106           {    s_transpose.begin(j);
01107                i = s_transpose.next_element();
01108                while( i != s_transpose.end() )
01109                {    row.push_back(i);
01110                     col.push_back(j);
01111                     sort_col.push_back(k);
01112                     k++;
01113                     i = s_transpose.next_element();
01114                }
01115           } 
01116           size_t K = k;
01117           VectorBase J(K);
01118           row.push_back(m);
01119           col.push_back(n);
01120           sort_col.push_back(K);
01121      
01122           // now we have folded this into the following case
01123           SparseJacobianFor(x, s_transpose, J, work);
01124 
01125           // now set the non-zero return values
01126           for(k = 0; k < K; k++)
01127                jac[row[sort_col[k]] * n + col[sort_col[k]]] = J[k];
01128      }
01129      else
01130      {    CppAD::vector<size_t>& row(work.user_row);
01131           CppAD::vector<size_t>& col(work.user_col);
01132           CppAD::vector<size_t>& sort_row(work.sort_row);
01133 
01134           // reverse mode, rows are sorted
01135           assert( row.size() == 0 );
01136           assert( col.size() == 0 );
01137           assert( sort_row.size() == 0 );
01138 
01139           // need an internal copy of sparsity pattern
01140           Pattern_type s;
01141           bool transpose = false;
01142           sparsity_user2internal(s, p, m, n, transpose);
01143 
01144           k = 0;
01145           for(i = 0; i < m; i++)
01146           {    s.begin(i);
01147                j = s.next_element();
01148                while( j != s.end() )
01149                {    row.push_back(i);
01150                     col.push_back(j);
01151                     sort_row.push_back(k);
01152                     k++;
01153                     j = s.next_element();
01154                }
01155           } 
01156           size_t K = k;
01157           VectorBase J(K);
01158           row.push_back(m);
01159           col.push_back(n);
01160           sort_row.push_back(K);
01161 
01162           // now we have folded this into the following case
01163           SparseJacobianRev(x, s, J, work);
01164 
01165           // now set the non-zero return values
01166           for(k = 0; k < K; k++)
01167                jac[row[sort_row[k]] * n + col[sort_row[k]]] = J[k];
01168      }
01169 
01170      return jac;
01171 } 
01172 
01173 /*!
01174 Compute a sparse Jacobian.
01175 
01176 The C++ source code corresponding to this operation is
01177 \verbatim
01178      jac = SparseJacobian(x)
01179 \endverbatim
01180 
01181 \tparam Base
01182 is the base type for the recording that is stored in this
01183 ADFun<Base object.
01184 
01185 \tparam VectorBase
01186 is a simple vector class with elements of the \a Base.
01187 
01188 \param x
01189 is a vector specifing the point at which to compute the Jacobian.
01190 
01191 \return
01192 Will be a vector of size \c m * n containing the Jacobian at the
01193 specified point (in row major order).
01194 */
01195 template <class Base>
01196 template <class VectorBase>
01197 VectorBase ADFun<Base>::SparseJacobian( const VectorBase& x )
01198 {    typedef CppAD::vectorBool   VectorBool;
01199 
01200      size_t m = Range();
01201      size_t n = Domain();
01202 
01203      // sparsity pattern for Jacobian
01204      VectorBool p(m * n);
01205 
01206      if( n <= m )
01207      {    size_t j, k;
01208 
01209           // use forward mode 
01210           VectorBool r(n * n);
01211           for(j = 0; j < n; j++)
01212           {    for(k = 0; k < n; k++)
01213                     r[j * n + k] = false;
01214                r[j * n + j] = true;
01215           }
01216           p = ForSparseJac(n, r);
01217      }
01218      else
01219      {    size_t i, k;
01220 
01221           // use reverse mode 
01222           VectorBool s(m * m);
01223           for(i = 0; i < m; i++)
01224           {    for(k = 0; k < m; k++)
01225                     s[i * m + k] = false;
01226                s[i * m + i] = true;
01227           }
01228           p = RevSparseJac(m, s);
01229      }
01230      return SparseJacobian(x, p);
01231 }
01232 
01233 
01234 /*! \} */
01235 CPPAD_END_NAMESPACE
01236 # endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines