CppAD: A C++ Algorithmic Differentiation Package  20171217
sparse_jac.hpp
Go to the documentation of this file.
1 # ifndef CPPAD_CORE_SPARSE_JAC_HPP
2 # define CPPAD_CORE_SPARSE_JAC_HPP
3 /* --------------------------------------------------------------------------
4 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-17 Bradley M. Bell
5
7 the terms of the
8  Eclipse Public License Version 1.0.
9
10 A copy of this license is included in the COPYING file of this distribution.
12 -------------------------------------------------------------------------- */
13
14 /*
15 $begin sparse_jac$$16 spell 17 Jacobian 18 Jacobians 19 const 20 jac 21 Taylor 22 rc 23 rcv 24 nr 25 nc 26 std 27 Cppad 28 Colpack 29 cmake 30$$ 31 32$section Computing Sparse Jacobians$$33 34 head Syntax$$
35 $icode%n_sweep% = %f%.sparse_jac_for( 36 %group_max%, %x%, %subset%, %pattern%, %coloring%, %work% 37 ) 38 %$$39 icode%n_sweep% = %f%.sparse_jac_rev( 40 %x%, %subset%, %pattern%, %coloring%, %work% 41 )%$$ 42 43 44$head Purpose$$45 We use latex F : \B{R}^n \rightarrow \B{R}^m$$ to denote the
46 function corresponding to $icode f$$. 47 Here icode n$$ is the$cref/domain/seq_property/Domain/$$size, 48 and icode m$$ is the $cref/range/seq_property/Range/$$size, or icode f$$. 49 The syntax above takes advantage of sparsity when computing the Jacobian 50$latex $51 J(x) = F^{(1)} (x) 52$ $$53 In the sparse case, this should be faster and take less memory than 54 cref Jacobian$$.
55 We use the notation $latex J_{i,j} (x)$$to denote the partial of 56 latex F_i (x)$$ with respect to$latex x_j$$. 57 58 head SizeVector$$
59 The type $icode SizeVector$$is a cref SimpleVector$$ class with 60$cref/elements of type/SimpleVector/Elements of Specified Type/$$61 code size_t$$.
62
63 $head BaseVector$$64 The type icode BaseVector$$ is a$cref SimpleVector$$class with 65 cref/elements of type/SimpleVector/Elements of Specified Type/$$
66 $code size_t$$. 67 68 head sparse_jac_for$$ 69 This function uses first order forward mode sweeps$cref forward_one$$70 to compute multiple columns of the Jacobian at the same time. 71 72 head sparse_jac_rev$$
73 This uses function first order reverse mode sweeps $cref reverse_one$$74 to compute multiple rows of the Jacobian at the same time. 75 76 head f$$ 77 This object has prototype 78$codei%
80 %$$81 Note that the Taylor coefficients stored in icode f$$ are affected
82 by this operation; see
83 $cref/uses forward/sparse_jac/Uses Forward/$$below. 84 85 head group_max$$ 86 This argument has prototype 87$codei%
88  size_t %group_max%
89 %$$90 and must be greater than zero. 91 It specifies the maximum number of colors to group during 92 a single forward sweep. 93 If a single color is in a group, 94 a single direction for of first order forward mode 95 cref forward_one$$ is used for each color.
96 If multiple colors are in a group,
97 the multiple direction for of first order forward mode
98 $cref forward_dir$$is used with one direction for each color. 99 This uses separate memory for each direction (more memory), 100 but my be significantly faster. 101 102 head x$$ 103 This argument has prototype 104$codei%
105  const %BaseVector%& %x%
106 %$$107 and its size is icode n$$.
108 It specifies the point at which to evaluate the Jacobian
109 $latex J(x)$$. 110 111 head subset$$ 112 This argument has prototype 113$codei%
114  sparse_rcv<%SizeVector%, %BaseVector%>& %subset%
115 %$$116 Its row size is icode%subset%.nr() == %m%$$,
117 and its column size is $icode%subset%.nc() == %n%$$. 118 It specifies which elements of the Jacobian are computed. 119 The input value of its value vector 120 icode%subset%.val()%$$ does not matter. 121 Upon return it contains the value of the corresponding elements 122 of the Jacobian. 123 All of the row, column pairs in$icode subset$$must also appear in 124 icode pattern$$; i.e., they must be possibly non-zero.
125
126 $head pattern$$127 This argument has prototype 128 codei% 129 const sparse_rc<%SizeVector%>& %pattern% 130 %$$ 131 Its row size is$icode%pattern%.nr() == %m%$$, 132 and its column size is icode%pattern%.nc() == %n%$$.
133 It is a sparsity pattern for the Jacobian $latex J(x)$$. 134 This argument is not used (and need not satisfy any conditions), 135 when cref/work/sparse_jac/work/$$ is non-empty. 136 137$head coloring$$138 The coloring algorithm determines which rows (reverse) or columns (forward) 139 can be computed during the same sweep. 140 This field has prototype 141 codei% 142 const std::string& %coloring% 143 %$$
144 This value only matters when work is empty; i.e.,
145 after the $icode work$$constructor or icode%work%.clear()%$$. 146 147$subhead cppad$$148 This uses a general purpose coloring algorithm written for Cppad. 149 150 subhead colpack$$
151 If $cref colpack_prefix$$is specified on the 152 cref/cmake command/cmake/CMake Command/$$ line, 153 you can set$icode coloring$$to code colpack$$.
154 This uses a general purpose coloring algorithm that is part of Colpack.
155
156 $head work$$157 This argument has prototype 158 codei% 159 sparse_jac_work& %work% 160 %$$ 161 We refer to its initial value, 162 and its value after$icode%work%.clear()%$$, as empty. 163 If it is empty, information is stored in icode work$$.
164 This can be used to reduce computation when
165 a future call is for the same object $icode f$$, 166 the same member function code sparse_jac_for$$ or$code sparse_jac_rev$$, 167 and the same subset of the Jacobian. 168 If any of these values change, use icode%work%.clear()%$$ to
169 empty this structure.
170
171 $head n_sweep$$172 The return value icode n_sweep$$ has prototype 173$codei%
174  size_t %n_sweep%
175 %$$176 If code sparse_jac_for$$ ($code sparse_jac_rev$$) is used, 177 icode n_sweep$$ is the number of first order forward (reverse) sweeps 178 used to compute the requested Jacobian values. 179 It is also the number of colors determined by the coloring method 180 mentioned above. 181 This is proportional to the total computational work, 182 not counting the zero order forward sweep, 183 or combining multiple columns (rows) into a single sweep. 184 185$head Uses Forward$$186 After each call to cref Forward$$,
187 the object $icode f$$contains the corresponding 188 cref/Taylor coefficients/glossary/Taylor Coefficient/$$. 189 After a call to$code sparse_jac_forward$$or code sparse_jac_rev$$,
190 the zero order coefficients correspond to
191 $codei% 192 %f%.Forward(0, %x%) 193 %$$194 All the other forward mode coefficients are unspecified. 195 196 head Example$$ 197$children%
198  example/sparse/sparse_jac_for.cpp%
199  example/sparse/sparse_jac_rev.cpp
200 %$$201 The files cref sparse_jac_for.cpp$$ and $cref sparse_jac_rev.cpp$$202 are examples and tests of code sparse_jac_for$$ and$code sparse_jac_rev$$. 203 They return code true$$, if they succeed, and $code false$$otherwise. 204 205$end
206 */
210 # include <cppad/utility/vector.hpp>
211
212 /*!
213 \file sparse_jac.hpp
214 Sparse Jacobian calculation routines.
215 */
217 /*!
218 Class used to hold information used by Sparse Jacobian routines in this file,
219 so they do not need to be recomputed every time.
220 */
222  public:
223  /// indices that sort the user row and col arrays by color
225  /// results of the coloring algorithm
227  //
228  /// constructor
230  { }
231  /// reset work to empty.
232  /// This informs CppAD that color and order need to be recomputed
233  void clear(void)
234  { order.clear();
235  color.clear();
236  }
237 };
238 // ----------------------------------------------------------------------------
239 /*!
240 Calculate sparse Jacobains using forward mode
241
242 \tparam Base
243 the base type for the recording that is stored in the ADFun object.
244
245 \tparam SizeVector
246 a simple vector class with elements of type size_t.
247
248 \tparam BaseVector
249 a simple vector class with elements of type Base.
250
251 \param group_max
252 specifies the maximum number of colors to group during a single forward sweep.
253 This must be greater than zero and group_max = 1 minimizes memory usage.
254
255 \param x
256 a vector of length n, the number of independent variables in f
257 (this ADFun object).
258
259 \param subset
260 specifices the subset of the sparsity pattern where the Jacobian is evaluated.
261 subset.nr() == m,
262 subset.nc() == n.
263
264 \param pattern
265 is a sparsity pattern for the Jacobian of f;
266 pattern.nr() == m,
267 pattern.nc() == n,
268 where m is number of dependent variables in f.
269
270 \param coloring
271 determines which coloring algorithm is used.
272 This must be cppad or colpack.
273
274 \param work
275 this structure must be empty, or contain the information stored
276 by a previous call to sparse_jac_for.
277 The previous call must be for the same ADFun object f
278 and the same subset.
279
280 \return
281 This is the number of first order forward sweeps used to compute
282 the Jacobian.
283 */
284 template <class Base>
285 template <class SizeVector, class BaseVector>
287  size_t group_max ,
288  const BaseVector& x ,
290  const sparse_rc<SizeVector>& pattern ,
291  const std::string& coloring ,
292  sparse_jac_work& work )
293 { size_t m = Range();
294  size_t n = Domain();
295  //
297  subset.nr() == m,
298  "sparse_jac_for: subset.nr() not equal range dimension for f"
299  );
301  subset.nc() == n,
302  "sparse_jac_for: subset.nc() not equal domain dimension for f"
303  );
304  //
305  // row and column vectors in subset
306  const SizeVector& row( subset.row() );
307  const SizeVector& col( subset.col() );
308  //
309  vector<size_t>& color(work.color);
310  vector<size_t>& order(work.order);
312  color.size() == 0 || color.size() == n,
313  "sparse_jac_for: work is non-empty and conditions have changed"
314  );
315  //
316  // point at which we are evaluationg the Jacobian
317  Forward(0, x);
318  //
319  // number of elements in the subset
320  size_t K = subset.nnz();
321  //
322  // check for case were there is nothing to do
323  // (except for call to Forward(0, x)
324  if( K == 0 )
325  return 0;
326  //
327  // check for case where input work is empty
328  if( color.size() == 0 )
329  { // compute work color and order vectors
331  pattern.nr() == m,
332  "sparse_jac_for: pattern.nr() not equal range dimension for f"
333  );
335  pattern.nc() == n,
336  "sparse_jac_for: pattern.nc() not equal domain dimension for f"
337  );
338  //
339  // convert pattern to an internal version of its transpose
340  vector<size_t> internal_index(n);
341  for(size_t j = 0; j < n; j++)
342  internal_index[j] = j;
343  bool transpose = true;
344  bool zero_empty = false;
345  bool input_empty = true;
346  local::sparse_list pattern_transpose;
347  pattern_transpose.resize(n, m);
348  local::set_internal_sparsity(zero_empty, input_empty,
349  transpose, internal_index, pattern_transpose, pattern
350  );
351  //
352  // execute coloring algorithm
353  // (we are using transpose because coloring groups rows, not columns).
354  color.resize(n);
355  if( coloring == "cppad" )
356  local::color_general_cppad(pattern_transpose, col, row, color);
357  else if( coloring == "colpack" )
358  {
359 # if CPPAD_HAS_COLPACK
360  local::color_general_colpack(pattern_transpose, col, row, color);
361 # else
363  false,
364  "sparse_jac_for: coloring = colpack "
365  "and colpack_prefix missing from cmake command line."
366  );
367 # endif
368  }
370  false,
371  "sparse_jac_for: coloring is not valid."
372  );
373  //
374  // put sorting indices in color order
375  SizeVector key(K);
376  order.resize(K);
377  for(size_t k = 0; k < K; k++)
378  key[k] = color[ col[k] ];
379  index_sort(key, order);
380  }
381  // Base versions of zero and one
382  Base one(1.0);
383  Base zero(0.0);
384  //
385  size_t n_color = 1;
386  for(size_t j = 0; j < n; j++) if( color[j] < n )
387  n_color = std::max(n_color, color[j] + 1);
388  //
389  // initialize the return Jacobian values as zero
390  for(size_t k = 0; k < K; k++)
391  subset.set(k, zero);
392  //
393  // index in subset
394  size_t k = 0;
395  // number of colors computed so far
396  size_t color_count = 0;
397  //
398  while( color_count < n_color )
399  { // number of colors that will be in this group
400  size_t group_size = std::min(group_max, n_color - color_count);
401  //
402  // forward mode values for independent and dependent variables
403  BaseVector dx(n * group_size), dy(m * group_size);
404  //
405  // set dx
406  for(size_t ell = 0; ell < group_size; ell++)
407  { // combine all columns with this color
408  for(size_t j = 0; j < n; j++)
409  { dx[j * group_size + ell] = zero;
410  if( color[j] == ell + color_count )
411  dx[j * group_size + ell] = one;
412  }
413  }
414  if( group_size == 1 )
415  dy = Forward(1, dx);
416  else
417  dy = Forward(1, group_size, dx);
418  //
419  // store results in subset
420  for(size_t ell = 0; ell < group_size; ell++)
421  { // color with index ell + color_count is in this group
422  while(k < K && color[ col[ order[k] ] ] == ell + color_count )
423  { // subset element with index order[k] is included in this color
424  size_t r = row[ order[k] ];
425  subset.set( order[k], dy[ r * group_size + ell ] );
426  ++k;
427  }
428  }
429  // advance color count
430  color_count += group_size;
431  }
432  CPPAD_ASSERT_UNKNOWN( color_count == n_color );
433  //
434  return n_color;
435 }
436 // ----------------------------------------------------------------------------
437 /*!
438 Calculate sparse Jacobains using reverse mode
439
440 \tparam Base
441 the base type for the recording that is stored in the ADFun object.
442
443 \tparam SizeVector
444 a simple vector class with elements of type size_t.
445
446 \tparam BaseVector
447 a simple vector class with elements of type Base.
448
449 \param x
450 a vector of length n, the number of independent variables in f
451 (this ADFun object).
452
453 \param subset
454 specifices the subset of the sparsity pattern where the Jacobian is evaluated.
455 subset.nr() == m,
456 subset.nc() == n.
457
458 \param pattern
459 is a sparsity pattern for the Jacobian of f;
460 pattern.nr() == m,
461 pattern.nc() == n,
462 where m is number of dependent variables in f.
463
464 \param coloring
465 determines which coloring algorithm is used.
466 This must be cppad or colpack.
467
468 \param work
469 this structure must be empty, or contain the information stored
470 by a previous call to sparse_jac_rev.
471 The previous call must be for the same ADFun object f
472 and the same subset.
473
474 \return
475 This is the number of first order reverse sweeps used to compute
476 the Jacobian.
477 */
478 template <class Base>
479 template <class SizeVector, class BaseVector>
481  const BaseVector& x ,
483  const sparse_rc<SizeVector>& pattern ,
484  const std::string& coloring ,
485  sparse_jac_work& work )
486 { size_t m = Range();
487  size_t n = Domain();
488  //
490  subset.nr() == m,
491  "sparse_jac_rev: subset.nr() not equal range dimension for f"
492  );
494  subset.nc() == n,
495  "sparse_jac_rev: subset.nc() not equal domain dimension for f"
496  );
497  //
498  // row and column vectors in subset
499  const SizeVector& row( subset.row() );
500  const SizeVector& col( subset.col() );
501  //
502  vector<size_t>& color(work.color);
503  vector<size_t>& order(work.order);
505  color.size() == 0 || color.size() == m,
506  "sparse_jac_rev: work is non-empty and conditions have changed"
507  );
508  //
509  // point at which we are evaluationg the Jacobian
510  Forward(0, x);
511  //
512  // number of elements in the subset
513  size_t K = subset.nnz();
514  //
515  // check for case were there is nothing to do
516  // (except for call to Forward(0, x)
517  if( K == 0 )
518  return 0;
519  //
520  // check for case where input work is empty
521  if( color.size() == 0 )
522  { // compute work color and order vectors
524  pattern.nr() == m,
525  "sparse_jac_rev: pattern.nr() not equal range dimension for f"
526  );
528  pattern.nc() == n,
529  "sparse_jac_rev: pattern.nc() not equal domain dimension for f"
530  );
531  //
532  // convert pattern to an internal version
533  vector<size_t> internal_index(m);
534  for(size_t i = 0; i < m; i++)
535  internal_index[i] = i;
536  bool transpose = false;
537  bool zero_empty = false;
538  bool input_empty = true;
539  local::sparse_list internal_pattern;
540  internal_pattern.resize(m, n);
541  local::set_internal_sparsity(zero_empty, input_empty,
542  transpose, internal_index, internal_pattern, pattern
543  );
544  //
545  // execute coloring algorithm
546  color.resize(m);
547  if( coloring == "cppad" )
548  local::color_general_cppad(internal_pattern, row, col, color);
549  else if( coloring == "colpack" )
550  {
551 # if CPPAD_HAS_COLPACK
552  local::color_general_colpack(internal_pattern, row, col, color);
553 # else
555  false,
556  "sparse_jac_rev: coloring = colpack "
557  "and colpack_prefix missing from cmake command line."
558  );
559 # endif
560  }
562  false,
563  "sparse_jac_rev: coloring is not valid."
564  );
565  //
566  // put sorting indices in color order
567  SizeVector key(K);
568  order.resize(K);
569  for(size_t k = 0; k < K; k++)
570  key[k] = color[ row[k] ];
571  index_sort(key, order);
572  }
573  // Base versions of zero and one
574  Base one(1.0);
575  Base zero(0.0);
576  //
577  size_t n_color = 1;
578  for(size_t i = 0; i < m; i++) if( color[i] < m )
579  n_color = std::max(n_color, color[i] + 1);
580  //
581  // initialize the return Jacobian values as zero
582  for(size_t k = 0; k < K; k++)
583  subset.set(k, zero);
584  //
585  // weighting vector and return values for calls to Reverse
586  BaseVector w(m), dw(n);
587  //
588  // loop over colors
589  size_t k = 0;
590  for(size_t ell = 0; ell < n_color; ell++)
591  if( k == K )
592  { // kludge because colpack returns colors that are not used
593  // (it does not know about the subset corresponding to row, col)
594  CPPAD_ASSERT_UNKNOWN( coloring == "colpack" );
595  }
596  else if( color[ row[ order[k] ] ] != ell )
597  { // kludge because colpack returns colors that are not used
598  // (it does not know about the subset corresponding to row, col)
599  CPPAD_ASSERT_UNKNOWN( coloring == "colpack" );
600  }
601  else
602  { CPPAD_ASSERT_UNKNOWN( color[ row[ order[k] ] ] == ell );
603  //
604  // combine all rows with this color
605  for(size_t i = 0; i < m; i++)
606  { w[i] = zero;
607  if( color[i] == ell )
608  w[i] = one;
609  }
610  // call reverse mode for all these rows at once
611  dw = Reverse(1, w);
612  //
613  // set the corresponding components of the result
614  while( k < K && color[ row[order[k]] ] == ell )
615  { subset.set(order[k], dw[col[order[k]]] );
616  k++;
617  }
618  }
619  return n_color;
620 }
621
622 } // END_CPPAD_NAMESPACE
623 # endif
Check that exp is true, if not print msg and terminate execution.
Vector of sets of positive integers, each set stored as a singly linked list.
Definition: sparse_list.hpp:35
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 ...
void set(size_t k, const value_type &v)
Definition: sparse_rcv.hpp:240
size_t nnz(void) const
number of possibly non-zero elements in matrix
Definition: sparse_rcv.hpp:254
const SizeVector & col(void) const
column indices
Definition: sparse_rcv.hpp:260
CppAD::vector< size_t > color
results of the coloring algorithm
Definition: sparse_jac.hpp:226
void clear(void)
free memory and set number of elements to zero
Definition: vector.hpp:418
Define the CppAD error checking macros (all of which begin with CPPAD_ASSERT_)
size_t nc(void) const
number of columns in matrix
Definition: sparse_rcv.hpp:251
Sparse matrices with elements of type Scalar.
Definition: sparse_rcv.hpp:212
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 nr(void) const
number of rows in matrix
Definition: sparse_rcv.hpp:248
size_t nr(void) const
number of rows in matrix
Definition: sparse_rc.hpp:284
CppAD::vector< size_t > order
indices that sort the user row and col arrays by color
Definition: sparse_jac.hpp:224
sparsity pattern for a matrix with indices of type size_t
Definition: sparse_rc.hpp:212
void index_sort(const VectorKey &keys, VectorSize &ind)
Compute the indices that sort a vector of keys.
Definition: index_sort.hpp:140
size_t sparse_jac_for(size_t group_max, const BaseVector &x, sparse_rcv< SizeVector, BaseVector > &subset, const sparse_rc< SizeVector > &pattern, const std::string &coloring, sparse_jac_work &work)
Calculate sparse Jacobains using forward mode.
Definition: sparse_jac.hpp:286
Check that exp is true, if not terminate execution.
void resize(size_t n_set, size_t end)
Start a new vector of sets.
Class used to hold information used by Sparse Jacobian routines in this file, so they do not need to ...
Definition: sparse_jac.hpp:221
File used to define CppAD::vector and CppAD::vectorBool.
size_t nc(void) const
number of columns in matrix
Definition: sparse_rc.hpp:287
void set_internal_sparsity(bool zero_empty, bool input_empty, bool transpose, const vector< size_t > &internal_index, InternalSparsity &internal_pattern, const sparse_rc< SizeVector > &pattern_in)
Update the internal sparsity pattern for a sub-set of rows.
Routines that enable code to be independent of which internal spasity pattern is used.
const SizeVector & row(void) const
row indices
Definition: sparse_rcv.hpp:257
size_t sparse_jac_rev(const BaseVector &x, sparse_rcv< SizeVector, BaseVector > &subset, const sparse_rc< SizeVector > &pattern, const std::string &coloring, sparse_jac_work &work)
Calculate sparse Jacobains using reverse mode.
Definition: sparse_jac.hpp:480
sparse_jac_work(void)
constructor
Definition: sparse_jac.hpp:229
void clear(void)
reset work to empty. This informs CppAD that color and order need to be recomputed ...
Definition: sparse_jac.hpp:233