CppAD: A C++ Algorithmic Differentiation Package  20171217
sparse_jacobian.hpp
Go to the documentation of this file.
3
4 /* --------------------------------------------------------------------------
6
8 the terms of the
9  Eclipse Public License Version 1.0.
10
11 A copy of this license is included in the COPYING file of this distribution.
13 -------------------------------------------------------------------------- */
14
15 // maximum number of sparse directions to compute at the same time
16
17 // # define CPPAD_SPARSE_JACOBIAN_MAX_MULTIPLE_DIRECTION 1
19
20 /*
21 $begin sparse_jacobian$$22 spell 23 cppad 24 colpack 25 cmake 26 recomputed 27 valarray 28 std 29 CppAD 30 Bool 31 jac 32 Jacobian 33 Jacobians 34 const 35 Taylor 36$$ 37 38$section Sparse Jacobian$$39 40 head Syntax$$
41 $icode%jac% = %f%.SparseJacobian(%x%) 42 %jac% = %f%.SparseJacobian(%x%, %p%) 43 %n_sweep% = %f%.SparseJacobianForward(%x%, %p%, %row%, %col%, %jac%, %work%) 44 %n_sweep% = %f%.SparseJacobianReverse(%x%, %p%, %row%, %col%, %jac%, %work%) 45 %$$46 47 head Purpose$$ 48 We use$latex n$$for the cref/domain/seq_property/Domain/$$ size,
49 and $latex m$$for the cref/range/seq_property/Range/$$ size of$icode f$$. 50 We use latex F : \B{R}^n \rightarrow \B{R}^m$$ do denote the
51 $cref/AD function/glossary/AD Function/$$52 corresponding to icode f$$. 53 The syntax above sets$icode jac$$to the Jacobian 54 latex $55 jac = F^{(1)} (x) 56$$$
57 This routine takes advantage of the sparsity of the Jacobian
58 in order to reduce the amount of computation necessary.
59 If $icode row$$and icode col$$ are present, it also takes 60 advantage of the reduced set of elements of the Jacobian that 61 need to be computed. 62 One can use speed tests (e.g.$cref speed_test$$) 63 to verify that results are computed faster 64 than when using the routine cref Jacobian$$.
65
66 $head f$$67 The object icode f$$ has prototype 68$codei%
70 %$$71 Note that the cref ADFun$$ object $icode f$$is not code const$$ 72 (see$cref/Uses Forward/sparse_jacobian/Uses Forward/$$below). 73 74 head x$$
75 The argument $icode x$$has prototype 76 codei% 77 const %VectorBase%& %x% 78 %$$ 79 (see$cref/VectorBase/sparse_jacobian/VectorBase/$$below) 80 and its size 81 must be equal to icode n$$, the dimension of the
82 $cref/domain/seq_property/Domain/$$space for icode f$$. 83 It specifies 84 that point at which to evaluate the Jacobian. 85 86$head p$$87 The argument icode p$$ is optional and has prototype
88 $codei% 89 const %VectorSet%& %p% 90 %$$91 (see cref/VectorSet/sparse_jacobian/VectorSet/$$ below). 92 If it has elements of type$code bool$$, 93 its size is latex m * n$$.
94 If it has elements of type $code std::set<size_t>$$, 95 its size is latex m$$ and all its set elements are between 96 zero and$latex n - 1$$. 97 It specifies a 98 cref/sparsity pattern/glossary/Sparsity Pattern/$$
99 for the Jacobian $latex F^{(1)} (x)$$. 100 pre 101 102$$ 103 If this sparsity pattern does not change between calls to 104$codei SparseJacobian$$, it should be faster to calculate icode p$$ once
105 (using $cref ForSparseJac$$or cref RevSparseJac$$) 106 and then pass$icode p$$to codei SparseJacobian$$.
107 Furthermore, if you specify $icode work$$in the calling sequence, 108 it is not necessary to keep the sparsity pattern; see the heading 109 cref/p/sparse_jacobian/work/p/$$ under the$icode work$$description. 110 pre 111 112$$
114 if you specify $icode p$$, CppAD will use the same 115 type of sparsity representation 116 (vectors of code bool$$ or vectors of$code std::set<size_t>$$) 117 for its internal calculations. 118 Otherwise, the representation 119 for the internal calculations is unspecified. 120 121 head row, col$$
122 The arguments $icode row$$and icode col$$ are optional and have prototype 123$codei%
124  const %VectorSize%& %row%
125  const %VectorSize%& %col%
126 %$$127 (see cref/VectorSize/sparse_jacobian/VectorSize/$$ below).
128 They specify which rows and columns of $latex F^{(1)} (x)$$are 129 computes and in what order. 130 Not all the non-zero entries in latex F^{(1)} (x)$$ need be computed, 131 but all the entries specified by$icode row$$and icode col$$
132 must be possibly non-zero in the sparsity pattern.
133 We use $latex K$$to denote the value icode%jac%.size()%$$ 134 which must also equal the size of$icode row$$and icode col$$.
135 Furthermore,
136 for $latex k = 0 , \ldots , K-1$$, it must hold that 137 latex row[k] < m$$ and$latex col[k] < n$$. 138 139 head jac$$
140 The result $icode jac$$has prototype 141 codei% 142 %VectorBase%& %jac% 143 %$$ 144 In the case where the arguments$icode row$$and icode col$$ are not present,
145 the size of $icode jac$$is latex m * n$$ and 146 for$latex i = 0 , \ldots , m-1$$, 147 latex j = 0 , \ldots , n-1$$,
148 $latex $149 jac [ i * n + j ] = \D{ F_i }{ x_j } (x) 150$ $$151 pre 152 153$$ 154 In the case where the arguments$icode row$$and icode col$$ are present,
155 we use $latex K$$to denote the size of icode jac$$. 156 The input value of its elements does not matter. 157 Upon return, for$latex k = 0 , \ldots , K - 1$$, 158 latex $159 jac [ k ] = \D{ F_i }{ x_j } (x) 160 \; , \; 161 \; {\rm where} \; 162 i = row[k] 163 \; {\rm and } \; 164 j = col[k] 165$$$
166
167 $head work$$168 If this argument is present, it has prototype 169 codei% 170 sparse_jacobian_work& %work% 171 %$$ 172 This object can only be used with the routines 173$code SparseJacobianForward$$and code SparseJacobianReverse$$.
174 During its the first use, information is stored in $icode work$$. 175 This is used to reduce the work done by future calls to the same mode 176 (forward or reverse), 177 the same icode f$$,$icode p$$, icode row$$, and $icode col$$. 178 If a future call is for a different mode, 179 or any of these values have changed, 180 you must first call icode%work%.clear()%$$ 181 to inform CppAD that this information needs to be recomputed. 182 183$subhead color_method$$184 The coloring algorithm determines which columns (forward mode) 185 or rows (reverse mode) can be computed during the same sweep. 186 This field has prototype 187 codei% 188 std::string %work%.color_method 189 %$$
190 and its default value (after a constructor or $code clear()$$) 191 is code "cppad"$$. 192 If$cref colpack_prefix$$is specified on the 193 cref/cmake command/cmake/CMake Command/$$ line,
194 you can set this method to $code "colpack"$$. 195 This value only matters on the first call to code sparse_jacobian$$ 196 that follows the$icode work$$constructor or a call to 197 icode%work%.clear()%$$.
198
199 $subhead p$$200 If icode work$$ is present, and it is not the first call after 201 its construction or a clear, 202 the sparsity pattern$icode p$$is not used. 203 This enables one to free the sparsity pattern 204 and still compute corresponding sparse Jacobians. 205 206 head n_sweep$$
207 The return value $icode n_sweep$$has prototype 208 codei% 209 size_t %n_sweep% 210 %$$ 211 If$code SparseJacobianForward$$(code SparseJacobianReverse$$) is used,
212 $icode n_sweep$$is the number of first order forward (reverse) sweeps 213 used to compute the requested Jacobian values. 214 (This is also the number of colors determined by the coloring method 215 mentioned above). 216 This is proportional to the total work that code SparseJacobian$$ does, 217 not counting the zero order forward sweep, 218 or the work to combine multiple columns (rows) into a single sweep. 219 220$head VectorBase$$221 The type icode VectorBase$$ must be a $cref SimpleVector$$class with 222 cref/elements of type/SimpleVector/Elements of Specified Type/$$ 223$icode Base$$. 224 The routine cref CheckSimpleVector$$ will generate an error message
225 if this is not the case.
226
227 $head VectorSet$$228 The type icode VectorSet$$ must be a$cref SimpleVector$$class with 229 cref/elements of type/SimpleVector/Elements of Specified Type/$$
230 $code bool$$or code std::set<size_t>$$; 231 see$cref/sparsity pattern/glossary/Sparsity Pattern/$$for a discussion 232 of the difference. 233 The routine cref CheckSimpleVector$$ will generate an error message
234 if this is not the case.
235
236 $subhead Restrictions$$237 If icode VectorSet$$ has elements of$code std::set<size_t>$$, 238 then icode%p%[%i%]%$$ must return a reference (not a copy) to the
239 corresponding set.
240 According to section 26.3.2.3 of the 1998 C++ standard,
241 $code std::valarray< std::set<size_t> >$$does not satisfy 242 this condition. 243 244 head VectorSize$$ 245 The type$icode VectorSize$$must be a cref SimpleVector$$ class with
246 $cref/elements of type/SimpleVector/Elements of Specified Type/$$247 code size_t$$. 248 The routine$cref CheckSimpleVector$$will generate an error message 249 if this is not the case. 250 251 head Uses Forward$$
252 After each call to $cref Forward$$, 253 the object icode f$$ contains the corresponding 254$cref/Taylor coefficients/glossary/Taylor Coefficient/$$. 255 After a call to any of the sparse Jacobian routines, 256 the zero order Taylor coefficients correspond to 257 icode%f%.Forward(0, %x%)%$$
258 and the other coefficients are unspecified.
259
260 After $code SparseJacobian$$, 261 the previous calls to cref Forward$$ are undefined. 262 263$head Example$$264 children% 265 example/sparse/sparse_jacobian.cpp 266 %$$
267 The routine
268 $cref sparse_jacobian.cpp$$269 is examples and tests of code sparse_jacobian$$. 270 It return$code true$$, if it succeeds and code false$$ otherwise.
271
272 \$end
273 ==============================================================================
274 */
277
279 /*!
280 \file sparse_jacobian.hpp
281 Sparse Jacobian driver routine and helper functions.
282 */
283 // ===========================================================================
284 /*!
285 class used by SparseJacobian to hold information so it does not need to be
286 recomputed.
287 */
289  public:
290  /// Coloring method: "cppad", or "colpack"
291  /// (this field is set by user)
292  std::string color_method;
293  /// indices that sort the user row and col arrays by color
295  /// results of the coloring algorithm
297
298  /// constructor
300  { }
301  /// reset coloring method to its default and
302  /// inform CppAD that color and order need to be recomputed
303  void clear(void)
305  order.clear();
306  color.clear();
307  }
308 };
309 // ===========================================================================
310 /*!
311 Private helper function forward mode cases
312
313 \tparam Base
314 is the base type for the recording that is stored in this
316
317 \tparam VectorBase
318 is a simple vector class with elements of type \a Base.
319
320 \tparam VectorSet
321 is either sparse_pack or sparse_list.
322
323 \tparam VectorSize
324 is a simple vector class with elements of type \c size_t.
325
326 \param x [in]
327 is a vector specifing the point at which to compute the Jacobian.
328
329 \param p_transpose [in]
330 If <code>work.color.size() != 0</code>,
331 then \c p_transpose is not used.
332 Otherwise, it is a
333 sparsity pattern for the transpose of the Jacobian of this ADFun<Base> object.
334 Note that we do not change the values in \c p_transpose,
335 but is not \c const because we use its iterator facility.
336
337 \param row [in]
338 is the vector of row indices for the returned Jacobian values.
339
340 \param col [in]
341 is the vector of columns indices for the returned Jacobian values.
342 It must have the same size as \c row.
343
344 \param jac [out]
345 is the vector of Jacobian values. We use \c K to denote the size of \c jac.
346 The return value <code>jac[k]</code> is the partial of the
347 <code>row[k]</code> range component of the function with respect
348 the the <code>col[k]</code> domain component of its argument.
349
350 \param work
351 <code>work.color_method</code> is an input. The rest of
352 this structure contains information that is computed by \c SparseJacobainFor.
353 If the sparsity pattern, \c row vector, or \c col vectors
354 are not the same between calls to \c SparseJacobianFor,
355 \c work.clear() must be called to reinitialize \c work.
356
357 \return
358 Is the number of first order forward sweeps used to compute the
359 requested Jacobian values. The total work, not counting the zero order
360 forward sweep, or the time to combine computations, is proportional to this
361 return value.
362 */
363 template<class Base>
364 template <class VectorBase, class VectorSet, class VectorSize>
366  const VectorBase& x ,
367  VectorSet& p_transpose ,
368  const VectorSize& row ,
369  const VectorSize& col ,
370  VectorBase& jac ,
371  sparse_jacobian_work& work )
372 {
373  size_t j, k, ell;
374
377
378  size_t m = Range();
379  size_t n = Domain();
380
381  // some values
382  const Base zero(0);
383  const Base one(1);
384
385  // check VectorBase is Simple Vector class with Base type elements
386  CheckSimpleVector<Base, VectorBase>();
387
388  CPPAD_ASSERT_UNKNOWN( size_t(x.size()) == n );
389  CPPAD_ASSERT_UNKNOWN( color.size() == 0 || color.size() == n );
390
391  // number of components of Jacobian that are required
392  size_t K = size_t(jac.size());
393  CPPAD_ASSERT_UNKNOWN( size_t( row.size() ) == K );
394  CPPAD_ASSERT_UNKNOWN( size_t( col.size() ) == K );
395
396  // Point at which we are evaluating the Jacobian
397  Forward(0, x);
398
399  // check for case where nothing (except Forward above) to do
400  if( K == 0 )
401  return 0;
402
403  if( color.size() == 0 )
404  {
405  CPPAD_ASSERT_UNKNOWN( p_transpose.n_set() == n );
406  CPPAD_ASSERT_UNKNOWN( p_transpose.end() == m );
407
408  // execute coloring algorithm
409  color.resize(n);
410  if( work.color_method == "cppad" )
412  else if( work.color_method == "colpack" )
413  {
415  local::color_general_colpack(p_transpose, col, row, color);
416 # else
418  false,
419  "SparseJacobianForward: work.color_method = colpack "
420  "and colpack_prefix missing from cmake command line."
421  );
422 # endif
423  }
425  false,
426  "SparseJacobianForward: work.color_method is not valid."
427  );
428
429  // put sorting indices in color order
430  VectorSize key(K);
431  order.resize(K);
432  for(k = 0; k < K; k++)
433  key[k] = color[ col[k] ];
434  index_sort(key, order);
435  }
436  size_t n_color = 1;
437  for(j = 0; j < n; j++) if( color[j] < n )
438  n_color = std::max(n_color, color[j] + 1);
439
440  // initialize the return value
441  for(k = 0; k < K; k++)
442  jac[k] = zero;
443
444 # if CPPAD_SPARSE_JACOBIAN_MAX_MULTIPLE_DIRECTION == 1
445  // direction vector and return values for calls to forward
446  VectorBase dx(n), dy(m);
447
448  // loop over colors
449  k = 0;
450  for(ell = 0; ell < n_color; ell++)
451  { CPPAD_ASSERT_UNKNOWN( color[ col[ order[k] ] ] == ell );
452
453  // combine all columns with this color
454  for(j = 0; j < n; j++)
455  { dx[j] = zero;
456  if( color[j] == ell )
457  dx[j] = one;
458  }
459  // call forward mode for all these columns at once
460  dy = Forward(1, dx);
461
462  // set the corresponding components of the result
463  while( k < K && color[ col[order[k]] ] == ell )
464  { jac[ order[k] ] = dy[row[order[k]]];
465  k++;
466  }
467  }
468 # else
469  // abbreviation for this value
471  CPPAD_ASSERT_UNKNOWN( max_r > 1 );
472
473  // count the number of colors done so far
474  size_t count_color = 0;
475  // count the sparse matrix entries done so far
476  k = 0;
477  while( count_color < n_color )
478  { // number of colors we will do this time
479  size_t r = std::min(max_r , n_color - count_color);
480  VectorBase dx(n * r), dy(m * r);
481
482  // loop over colors we will do this tme
483  for(ell = 0; ell < r; ell++)
484  { // combine all columns with this color
485  for(j = 0; j < n; j++)
486  { dx[j * r + ell] = zero;
487  if( color[j] == ell + count_color )
488  dx[j * r + ell] = one;
489  }
490  }
491  size_t q = 1;
492  dy = Forward(q, r, dx);
493
494  // store results
495  for(ell = 0; ell < r; ell++)
496  { // set the components of the result for this color
497  while( k < K && color[ col[order[k]] ] == ell + count_color )
498  { jac[ order[k] ] = dy[ row[order[k]] * r + ell ];
499  k++;
500  }
501  }
502  count_color += r;
503  }
504 # endif
505  return n_color;
506 }
507 /*!
508 Private helper function for reverse mode cases.
509
510 \tparam Base
511 is the base type for the recording that is stored in this
513
514 \tparam VectorBase
515 is a simple vector class with elements of type \a Base.
516
517 \tparam VectorSet
518 is either sparse_pack or sparse_list.
519
520 \tparam VectorSize
521 is a simple vector class with elements of type \c size_t.
522
523 \param x [in]
524 is a vector specifing the point at which to compute the Jacobian.
525
526 \param p [in]
527 If <code>work.color.size() != 0</code>, then \c p is not used.
528 Otherwise, it is a
529 sparsity pattern for the Jacobian of this ADFun<Base> object.
530 Note that we do not change the values in \c p,
531 but is not \c const because we use its iterator facility.
532
533 \param row [in]
534 is the vector of row indices for the returned Jacobian values.
535
536 \param col [in]
537 is the vector of columns indices for the returned Jacobian values.
538 It must have the same size as \c row.
539
540 \param jac [out]
541 is the vector of Jacobian values.
542 It must have the same size as \c row.
543 The return value <code>jac[k]</code> is the partial of the
544 <code>row[k]</code> range component of the function with respect
545 the the <code>col[k]</code> domain component of its argument.
546
547 \param work
548 <code>work.color_method</code> is an input. The rest of
549 This structure contains information that is computed by \c SparseJacobainRev.
550 If the sparsity pattern, \c row vector, or \c col vectors
551 are not the same between calls to \c SparseJacobianRev,
552 \c work.clear() must be called to reinitialize \c work.
553
554 \return
555 Is the number of first order reverse sweeps used to compute the
556 reverse Jacobian values. The total work, not counting the zero order
557 forward sweep, or the time to combine computations, is proportional to this
558 return value.
559 */
560 template<class Base>
561 template <class VectorBase, class VectorSet, class VectorSize>
563  const VectorBase& x ,
564  VectorSet& p ,
565  const VectorSize& row ,
566  const VectorSize& col ,
567  VectorBase& jac ,
568  sparse_jacobian_work& work )
569 {
570  size_t i, k, ell;
571
574
575  size_t m = Range();
576  size_t n = Domain();
577
578  // some values
579  const Base zero(0);
580  const Base one(1);
581
582  // check VectorBase is Simple Vector class with Base type elements
583  CheckSimpleVector<Base, VectorBase>();
584
585  CPPAD_ASSERT_UNKNOWN( size_t(x.size()) == n );
586  CPPAD_ASSERT_UNKNOWN (color.size() == m || color.size() == 0 );
587
588  // number of components of Jacobian that are required
589  size_t K = size_t(jac.size());
590  CPPAD_ASSERT_UNKNOWN( size_t( size_t( row.size() ) ) == K );
591  CPPAD_ASSERT_UNKNOWN( size_t( size_t( col.size() ) ) == K );
592
593  // Point at which we are evaluating the Jacobian
594  Forward(0, x);
595
596  // check for case where nothing (except Forward above) to do
597  if( K == 0 )
598  return 0;
599
600  if( color.size() == 0 )
601  {
602  CPPAD_ASSERT_UNKNOWN( p.n_set() == m );
603  CPPAD_ASSERT_UNKNOWN( p.end() == n );
604
605  // execute the coloring algorithm
606  color.resize(m);
607  if( work.color_method == "cppad" )
609  else if( work.color_method == "colpack" )
610  {
612  local::color_general_colpack(p, row, col, color);
613 # else
615  false,
616  "SparseJacobianReverse: work.color_method = colpack "
617  "and colpack_prefix missing from cmake command line."
618  );
619 # endif
620  }
622  false,
623  "SparseJacobianReverse: work.color_method is not valid."
624  );
625
626  // put sorting indices in color order
627  VectorSize key(K);
628  order.resize(K);
629  for(k = 0; k < K; k++)
630  key[k] = color[ row[k] ];
631  index_sort(key, order);
632  }
633  size_t n_color = 1;
634  for(i = 0; i < m; i++) if( color[i] < m )
635  n_color = std::max(n_color, color[i] + 1);
636
637  // weighting vector for calls to reverse
638  VectorBase w(m);
639
640  // location for return values from Reverse
641  VectorBase dw(n);
642
643  // initialize the return value
644  for(k = 0; k < K; k++)
645  jac[k] = zero;
646
647  // loop over colors
648  k = 0;
649  for(ell = 0; ell < n_color; ell++)
650  { CPPAD_ASSERT_UNKNOWN( color[ row[ order[k] ] ] == ell );
651
652  // combine all the rows with this color
653  for(i = 0; i < m; i++)
654  { w[i] = zero;
655  if( color[i] == ell )
656  w[i] = one;
657  }
658  // call reverse mode for all these rows at once
659  dw = Reverse(1, w);
660
661  // set the corresponding components of the result
662  while( k < K && color[ row[order[k]] ] == ell )
663  { jac[ order[k] ] = dw[col[order[k]]];
664  k++;
665  }
666  }
667  return n_color;
668 }
669 // ==========================================================================
670 // Public Member functions
671 // ==========================================================================
672 /*!
673 Compute user specified subset of a sparse Jacobian using forward mode.
674
675 The C++ source code corresponding to this operation is
676 \verbatim
677  SparceJacobianForward(x, p, row, col, jac, work)
678 \endverbatim
679
680 \tparam Base
681 is the base type for the recording that is stored in this
683
684 \tparam VectorBase
685 is a simple vector class with elements of type \a Base.
686
687 \tparam VectorSet
688 is a simple vector class with elements of type
689 \c bool or \c std::set<size_t>.
690
691 \tparam VectorSize
692 is a simple vector class with elements of type \c size_t.
693
694 \param x [in]
695 is a vector specifing the point at which to compute the Jacobian.
696
697 \param p [in]
698 is the sparsity pattern for the Jacobian that we are calculating.
699
700 \param row [in]
701 is the vector of row indices for the returned Jacobian values.
702
703 \param col [in]
704 is the vector of columns indices for the returned Jacobian values.
705 It must have the same size as \c row.
706
707 \param jac [out]
708 is the vector of Jacobian values.
709 It must have the same size as \c row.
710 The return value <code>jac[k]</code> is the partial of the
711 <code>row[k]</code> range component of the function with respect
712 the the <code>col[k]</code> domain component of its argument.
713
714 \param work [in,out]
715 this structure contains information that depends on the function object,
716 sparsity pattern, \c row vector, and \c col vector.
717 If they are not the same between calls to \c SparseJacobianForward,
718 \c work.clear() must be called to reinitialize them.
719
720 \return
721 Is the number of first order forward sweeps used to compute the
722 requested Jacobian values. The total work, not counting the zero order
723 forward sweep, or the time to combine computations, is proportional to this
724 return value.
725 */
726 template<class Base>
727 template <class VectorBase, class VectorSet, class VectorSize>
729  const VectorBase& x ,
730  const VectorSet& p ,
731  const VectorSize& row ,
732  const VectorSize& col ,
733  VectorBase& jac ,
734  sparse_jacobian_work& work )
735 {
736  size_t n = Domain();
737  size_t m = Range();
738  size_t K = jac.size();
739 # ifndef NDEBUG
740  size_t k;
742  size_t(x.size()) == n ,
743  "SparseJacobianForward: size of x not equal domain dimension for f."
744  );
746  size_t(row.size()) == K && size_t(col.size()) == K ,
747  "SparseJacobianForward: either r or c does not have "
748  "the same size as jac."
749  );
751  work.color.size() == 0 || work.color.size() == n,
752  "SparseJacobianForward: invalid value in work."
753  );
754  for(k = 0; k < K; k++)
756  row[k] < m,
757  "SparseJacobianForward: invalid value in r."
758  );
760  col[k] < n,
761  "SparseJacobianForward: invalid value in c."
762  );
763  }
764  if( work.color.size() != 0 )
765  for(size_t j = 0; j < n; j++) CPPAD_ASSERT_KNOWN(
766  work.color[j] <= n,
767  "SparseJacobianForward: invalid value in work."
768  );
769 # endif
770  // check for case where there is nothing to compute
771  size_t n_sweep = 0;
772  if( K == 0 )
773  return n_sweep;
774
775  typedef typename VectorSet::value_type Set_type;
776  typedef typename local::internal_sparsity<Set_type>::pattern_type Pattern_type;
777  Pattern_type s_transpose;
778  if( work.color.size() == 0 )
779  { bool transpose = true;
780  const char* error_msg = "SparseJacobianForward: transposed sparsity"
781  " pattern does not have proper row or column dimension";
782  sparsity_user2internal(s_transpose, p, n, m, transpose, error_msg);
783  }
784  n_sweep = SparseJacobianFor(x, s_transpose, row, col, jac, work);
785  return n_sweep;
786 }
787 /*!
788 Compute user specified subset of a sparse Jacobian using forward mode.
789
790 The C++ source code corresponding to this operation is
791 \verbatim
792  SparceJacobianReverse(x, p, row, col, jac, work)
793 \endverbatim
794
795 \tparam Base
796 is the base type for the recording that is stored in this
798
799 \tparam VectorBase
800 is a simple vector class with elements of type \a Base.
801
802 \tparam VectorSet
803 is a simple vector class with elements of type
804 \c bool or \c std::set<size_t>.
805
806 \tparam VectorSize
807 is a simple vector class with elements of type \c size_t.
808
809 \param x [in]
810 is a vector specifing the point at which to compute the Jacobian.
811
812 \param p [in]
813 is the sparsity pattern for the Jacobian that we are calculating.
814
815 \param row [in]
816 is the vector of row indices for the returned Jacobian values.
817
818 \param col [in]
819 is the vector of columns indices for the returned Jacobian values.
820 It must have the same size as \c row.
821
822 \param jac [out]
823 is the vector of Jacobian values.
824 It must have the same size as \c row.
825 The return value <code>jac[k]</code> is the partial of the
826 <code>row[k]</code> range component of the function with respect
827 the the <code>col[k]</code> domain component of its argument.
828
829 \param work [in,out]
830 this structure contains information that depends on the function object,
831 sparsity pattern, \c row vector, and \c col vector.
832 If they are not the same between calls to \c SparseJacobianReverse,
833 \c work.clear() must be called to reinitialize them.
834
835 \return
836 Is the number of first order reverse sweeps used to compute the
837 reverse Jacobian values. The total work, not counting the zero order
838 forward sweep, or the time to combine computations, is proportional to this
839 return value.
840 */
841 template<class Base>
842 template <class VectorBase, class VectorSet, class VectorSize>
844  const VectorBase& x ,
845  const VectorSet& p ,
846  const VectorSize& row ,
847  const VectorSize& col ,
848  VectorBase& jac ,
849  sparse_jacobian_work& work )
850 {
851  size_t m = Range();
852  size_t n = Domain();
853  size_t K = jac.size();
854 # ifndef NDEBUG
855  size_t k;
857  size_t(x.size()) == n ,
858  "SparseJacobianReverse: size of x not equal domain dimension for f."
859  );
861  size_t(row.size()) == K && size_t(col.size()) == K ,
862  "SparseJacobianReverse: either r or c does not have "
863  "the same size as jac."
864  );
866  work.color.size() == 0 || work.color.size() == m,
867  "SparseJacobianReverse: invalid value in work."
868  );
869  for(k = 0; k < K; k++)
871  row[k] < m,
872  "SparseJacobianReverse: invalid value in r."
873  );
875  col[k] < n,
876  "SparseJacobianReverse: invalid value in c."
877  );
878  }
879  if( work.color.size() != 0 )
880  for(size_t i = 0; i < m; i++) CPPAD_ASSERT_KNOWN(
881  work.color[i] <= m,
882  "SparseJacobianReverse: invalid value in work."
883  );
884 # endif
885  // check for case where there is nothing to compute
886  size_t n_sweep = 0;
887  if( K == 0 )
888  return n_sweep;
889
890  typedef typename VectorSet::value_type Set_type;
891  typedef typename local::internal_sparsity<Set_type>::pattern_type Pattern_type;
892  Pattern_type s;
893  if( work.color.size() == 0 )
894  { bool transpose = false;
895  const char* error_msg = "SparseJacobianReverse: sparsity"
896  " pattern does not have proper row or column dimension";
897  sparsity_user2internal(s, p, m, n, transpose, error_msg);
898  }
899  n_sweep = SparseJacobianRev(x, s, row, col, jac, work);
900  return n_sweep;
901 }
902 /*!
903 Compute a sparse Jacobian.
904
905 The C++ source code corresponding to this operation is
906 \verbatim
907  jac = SparseJacobian(x, p)
908 \endverbatim
909
910 \tparam Base
911 is the base type for the recording that is stored in this
913
914 \tparam VectorBase
915 is a simple vector class with elements of type \a Base.
916
917 \tparam VectorSet
918 is a simple vector class with elements of type
919 \c bool or \c std::set<size_t>.
920
921 \param x [in]
922 is a vector specifing the point at which to compute the Jacobian.
923
924 \param p [in]
925 is the sparsity pattern for the Jacobian that we are calculating.
926
927 \return
928 Will be a vector if size \c m * n containing the Jacobian at the
929 specified point (in row major order).
930 */
931 template <class Base>
932 template <class VectorBase, class VectorSet>
934  const VectorBase& x, const VectorSet& p
935 )
936 { size_t i, j, k;
937
938  size_t m = Range();
939  size_t n = Domain();
940  VectorBase jac(m * n);
941
943  size_t(x.size()) == n,
944  "SparseJacobian: size of x not equal domain size for f."
945  );
946  CheckSimpleVector<Base, VectorBase>();
947
948  typedef typename VectorSet::value_type Set_type;
949  typedef typename local::internal_sparsity<Set_type>::pattern_type Pattern_type;
950
951  // initialize the return value as zero
952  Base zero(0);
953  for(i = 0; i < m; i++)
954  for(j = 0; j < n; j++)
955  jac[i * n + j] = zero;
956
960  if( n <= m )
961  {
962  // need an internal copy of sparsity pattern
963  Pattern_type s_transpose;
964  bool transpose = true;
965  const char* error_msg = "SparseJacobian: transposed sparsity"
966  " pattern does not have proper row or column dimension";
967  sparsity_user2internal(s_transpose, p, n, m, transpose, error_msg);
968
969  k = 0;
970  for(j = 0; j < n; j++)
971  { typename Pattern_type::const_iterator itr(s_transpose, j);
972  i = *itr;
973  while( i != s_transpose.end() )
974  { row.push_back(i);
975  col.push_back(j);
976  k++;
977  i = *(++itr);
978  }
979  }
980  size_t K = k;
981  VectorBase J(K);
982
983  // now we have folded this into the following case
984  SparseJacobianFor(x, s_transpose, row, col, J, work);
985
986  // now set the non-zero return values
987  for(k = 0; k < K; k++)
988  jac[ row[k] * n + col[k] ] = J[k];
989  }
990  else
991  {
992  // need an internal copy of sparsity pattern
993  Pattern_type s;
994  bool transpose = false;
995  const char* error_msg = "SparseJacobian: sparsity"
996  " pattern does not have proper row or column dimension";
997  sparsity_user2internal(s, p, m, n, transpose, error_msg);
998
999  k = 0;
1000  for(i = 0; i < m; i++)
1001  { typename Pattern_type::const_iterator itr(s, i);
1002  j = *itr;
1003  while( j != s.end() )
1004  { row.push_back(i);
1005  col.push_back(j);
1006  k++;
1007  j = *(++itr);
1008  }
1009  }
1010  size_t K = k;
1011  VectorBase J(K);
1012
1013  // now we have folded this into the following case
1014  SparseJacobianRev(x, s, row, col, J, work);
1015
1016  // now set the non-zero return values
1017  for(k = 0; k < K; k++)
1018  jac[ row[k] * n + col[k] ] = J[k];
1019  }
1020
1021  return jac;
1022 }
1023
1024 /*!
1025 Compute a sparse Jacobian.
1026
1027 The C++ source code corresponding to this operation is
1028 \verbatim
1029  jac = SparseJacobian(x)
1030 \endverbatim
1031
1032 \tparam Base
1033 is the base type for the recording that is stored in this
1035
1036 \tparam VectorBase
1037 is a simple vector class with elements of the \a Base.
1038
1039 \param x [in]
1040 is a vector specifing the point at which to compute the Jacobian.
1041
1042 \return
1043 Will be a vector of size \c m * n containing the Jacobian at the
1044 specified point (in row major order).
1045 */
1046 template <class Base>
1047 template <class VectorBase>
1048 VectorBase ADFun<Base>::SparseJacobian( const VectorBase& x )
1050
1051  size_t m = Range();
1052  size_t n = Domain();
1053
1054  // sparsity pattern for Jacobian
1055  VectorBool p(m * n);
1056
1057  if( n <= m )
1058  { size_t j, k;
1059
1060  // use forward mode
1061  VectorBool r(n * n);
1062  for(j = 0; j < n; j++)
1063  { for(k = 0; k < n; k++)
1064  r[j * n + k] = false;
1065  r[j * n + j] = true;
1066  }
1067  p = ForSparseJac(n, r);
1068  }
1069  else
1070  { size_t i, k;
1071
1072  // use reverse mode
1073  VectorBool s(m * m);
1074  for(i = 0; i < m; i++)
1075  { for(k = 0; k < m; k++)
1076  s[i * m + k] = false;
1077  s[i * m + i] = true;
1078  }
1079  p = RevSparseJac(m, s);
1080  }
1081  return SparseJacobian(x, p);
1082 }
1083
1086 # endif
Check that exp is true, if not print msg and terminate execution.
sparse_jacobian_work(void)
constructor
size_t SparseJacobianReverse(const VectorBase &x, const VectorSet &p, const VectorSize &r, const VectorSize &c, VectorBase &jac, sparse_jacobian_work &work)
Compute user specified subset of a sparse Jacobian using forward mode.
void color_general_cppad(const VectorSet &pattern, const VectorSize &row, const VectorSize &col, CppAD::vector< size_t > &color)
Determine which rows of a general sparse matrix can be computed together; i.e., do not have non-zero ...
std::string color_method
Coloring method: &quot;cppad&quot;, or &quot;colpack&quot; (this field is set by user)
results of the coloring algorithm
void clear(void)
free memory and set number of elements to zero
Definition: vector.hpp:418
VectorBase SparseJacobian(const VectorBase &x)
calculate sparse Jacobians
size_t SparseJacobianForward(const VectorBase &x, const VectorSet &p, const VectorSize &r, const VectorSize &c, VectorBase &jac, sparse_jacobian_work &work)
Compute user specified subset of a sparse Jacobian using forward mode.
Coloring algorithm for a general sparse matrix.
void resize(size_t n)
change the number of elements in this vector.
Definition: vector.hpp:399
size_t size(void) const
number of elements currently in this vector.
Definition: vector.hpp:387
void clear(void)
reset coloring method to its default and inform CppAD that color and order need to be recomputed ...
Two constant standard sets (currently used for concept checking).
void index_sort(const VectorKey &keys, VectorSize &ind)
Compute the indices that sort a vector of keys.
Definition: index_sort.hpp:140
size_t SparseJacobianFor(const VectorBase &x, VectorSet &p_transpose, const VectorSize &row, const VectorSize &col, VectorBase &jac, sparse_jacobian_work &work)
Private helper function forward mode cases.
void push_back(const Type &s)
add an element to the back of this vector
Definition: vector.hpp:491
class used by SparseJacobian to hold information so it does not need to be recomputed.