CppAD: A C++ Algorithmic Differentiation Package  20171217
eigen_mat_mul.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 /*
16 $begin atomic_eigen_mat_mul.hpp$$17 spell 18 Eigen 19 Taylor 20 nr 21 nc 22$$ 23 24$section Atomic Eigen Matrix Multiply Class$$25 26 head See Also$$
27 $cref atomic_mat_mul.hpp$$28 29 head Purpose$$ 30 Construct an atomic operation that computes the matrix product, 31$latex R = A \times B$$32 for any positive integers latex r$$, $latex m$$, latex c$$, 33 and any$latex A \in \B{R}^{r \times m}$$, 34 latex B \in \B{R}^{m \times c}$$.
35
36 $head Matrix Dimensions$$37 This example puts the matrix dimensions in the atomic function arguments, 38 instead of the cref/constructor/atomic_ctor/$$, so that they can 39 be different for different calls to the atomic function. 40 These dimensions are: 41$table
42 $icode nr_left$$43 cnext number of rows in the left matrix; i.e, latex r$$$rend
44 $icode n_middle$$45 cnext rows in the left matrix and columns in right; i.e, latex m$$$rend
46 $icode nc_right$$47 cnext number of columns in the right matrix; i.e., latex c$$ 48$tend
49
50 $head Theory$$51 52 subhead Forward$$ 53 For$latex k = 0 , \ldots $$, the th k$$ order Taylor coefficient
54 $latex R_k$$is given by 55 latex $56 R_k = \sum_{\ell = 0}^{k} A_\ell B_{k-\ell} 57$$$ 58 59$subhead Product of Two Matrices$$60 Suppose latex \bar{E}$$ is the derivative of the
61 scalar value function $latex s(E)$$with respect to latex E$$; i.e., 62$latex $63 \bar{E}_{i,j} = \frac{ \partial s } { \partial E_{i,j} } 64$ $$65 Also suppose that latex t$$ is a scalar valued argument and
66 $latex $67 E(t) = C(t) D(t) 68$ $$69 It follows that 70 latex $71 E'(t) = C'(t) D(t) + C(t) D'(t) 72$$$ 73 74$latex $75 (s \circ E)'(t) 76 = 77 \R{tr} [ \bar{E}^\R{T} E'(t) ] 78$ $$79 latex $80 = 81 \R{tr} [ \bar{E}^\R{T} C'(t) D(t) ] + 82 \R{tr} [ \bar{E}^\R{T} C(t) D'(t) ] 83$$$
84 $latex $85 = 86 \R{tr} [ D(t) \bar{E}^\R{T} C'(t) ] + 87 \R{tr} [ \bar{E}^\R{T} C(t) D'(t) ] 88$ $$89 latex $90 \bar{C} = \bar{E} D^\R{T} \W{,} 91 \bar{D} = C^\R{T} \bar{E} 92$$$ 93 94$subhead Reverse$$95 Reverse mode eliminates latex R_k$$ as follows:
96 for $latex \ell = 0, \ldots , k-1$$, 97 latex $98 \bar{A}_\ell = \bar{A}_\ell + \bar{R}_k B_{k-\ell}^\R{T} 99$$$ 100$latex $101 \bar{B}_{k-\ell} = \bar{B}_{k-\ell} + A_\ell^\R{T} \bar{R}_k 102$ $$103 104 105 nospell 106 107 head Start Class Definition$$
108 $srccode%cpp% */ 109 # include <cppad/cppad.hpp> 110 # include <Eigen/Core> 111 112 113 /* %$$114 head Public$$ 115 116$subhead Types$$117 srccode%cpp% */ 118 namespace { // BEGIN_EMPTY_NAMESPACE 119 120 template <class Base> 122 public: 123 // ----------------------------------------------------------- 124 // type of elements during calculation of derivatives 125 typedef Base scalar; 126 // type of elements during taping 128 // type of matrix during calculation of derivatives 129 typedef Eigen::Matrix< 130 scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> matrix; 131 // type of matrix during taping 132 typedef Eigen::Matrix< 133 ad_scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > ad_matrix; 134 /* %$$
135 $subhead Constructor$$136 srccode%cpp% */ 137 // constructor 138 atomic_eigen_mat_mul(void) : CppAD::atomic_base<Base>( 139 "atom_eigen_mat_mul" , 140 CppAD::atomic_base<Base>::set_sparsity_enum 141 ) 142 { } 143 /* %$$ 144$subhead op$$145 srccode%cpp% */ 146 // use atomic operation to multiply two AD matrices 148 const ad_matrix& left , 149 const ad_matrix& right ) 150 { size_t nr_left = size_t( left.rows() ); 151 size_t n_middle = size_t( left.cols() ); 152 size_t nc_right = size_t( right.cols() ); 153 assert( n_middle == size_t( right.rows() ) ); 154 size_t nx = 3 + (nr_left + nc_right) * n_middle; 155 size_t ny = nr_left * nc_right; 156 size_t n_left = nr_left * n_middle; 157 size_t n_right = n_middle * nc_right; 158 size_t n_result = nr_left * nc_right; 159 // 160 assert( 3 + n_left + n_right == nx ); 161 assert( n_result == ny ); 162 // ----------------------------------------------------------------- 163 // packed version of left and right 164 CPPAD_TESTVECTOR(ad_scalar) packed_arg(nx); 165 // 166 packed_arg[0] = ad_scalar( nr_left ); 167 packed_arg[1] = ad_scalar( n_middle ); 168 packed_arg[2] = ad_scalar( nc_right ); 169 for(size_t i = 0; i < n_left; i++) 170 packed_arg[3 + i] = left.data()[i]; 171 for(size_t i = 0; i < n_right; i++) 172 packed_arg[ 3 + n_left + i ] = right.data()[i]; 173 // ------------------------------------------------------------------ 174 // Packed version of result = left * right. 175 // This as an atomic_base funciton call that CppAD uses 176 // to store the atomic operation on the tape. 177 CPPAD_TESTVECTOR(ad_scalar) packed_result(ny); 178 (*this)(packed_arg, packed_result); 179 // ------------------------------------------------------------------ 180 // unpack result matrix 181 ad_matrix result(nr_left, nc_right); 182 for(size_t i = 0; i < n_result; i++) 183 result.data()[i] = packed_result[ i ]; 184 // 185 return result; 186 } 187 /* %$$
188 $head Private$$189 190 subhead Variables$$ 191$srccode%cpp% */
192 private:
193  // -------------------------------------------------------------
194  // one forward mode vector of matrices for left, right, and result
196  // one reverse mode vector of matrices for left, right, and result
198  // -------------------------------------------------------------
199 /* %$$200 subhead forward$$
201 $srccode%cpp% */ 202 // forward mode routine called by CppAD 203 virtual bool forward( 204 // lowest order Taylor coefficient we are evaluating 205 size_t p , 206 // highest order Taylor coefficient we are evaluating 207 size_t q , 208 // which components of x are variables 209 const CppAD::vector<bool>& vx , 210 // which components of y are variables 211 CppAD::vector<bool>& vy , 212 // tx [ 3 + j * (q+1) + k ] is x_j^k 213 const CppAD::vector<scalar>& tx , 214 // ty [ i * (q+1) + k ] is y_i^k 216 ) 217 { size_t n_order = q + 1; 218 size_t nr_left = size_t( CppAD::Integer( tx[ 0 * n_order + 0 ] ) ); 219 size_t n_middle = size_t( CppAD::Integer( tx[ 1 * n_order + 0 ] ) ); 220 size_t nc_right = size_t( CppAD::Integer( tx[ 2 * n_order + 0 ] ) ); 221 # ifndef NDEBUG 222 size_t nx = 3 + (nr_left + nc_right) * n_middle; 223 size_t ny = nr_left * nc_right; 224 # endif 225 // 226 assert( vx.size() == 0 || nx == vx.size() ); 227 assert( vx.size() == 0 || ny == vy.size() ); 228 assert( nx * n_order == tx.size() ); 229 assert( ny * n_order == ty.size() ); 230 // 231 size_t n_left = nr_left * n_middle; 232 size_t n_right = n_middle * nc_right; 233 size_t n_result = nr_left * nc_right; 234 assert( 3 + n_left + n_right == nx ); 235 assert( n_result == ny ); 236 // 237 // ------------------------------------------------------------------- 238 // make sure f_left_, f_right_, and f_result_ are large enough 239 assert( f_left_.size() == f_right_.size() ); 240 assert( f_left_.size() == f_result_.size() ); 241 if( f_left_.size() < n_order ) 242 { f_left_.resize(n_order); 243 f_right_.resize(n_order); 244 f_result_.resize(n_order); 245 // 246 for(size_t k = 0; k < n_order; k++) 247 { f_left_[k].resize(nr_left, n_middle); 248 f_right_[k].resize(n_middle, nc_right); 249 f_result_[k].resize(nr_left, nc_right); 250 } 251 } 252 // ------------------------------------------------------------------- 253 // unpack tx into f_left and f_right 254 for(size_t k = 0; k < n_order; k++) 255 { // unpack left values for this order 256 for(size_t i = 0; i < n_left; i++) 257 f_left_[k].data()[i] = tx[ (3 + i) * n_order + k ]; 258 // 259 // unpack right values for this order 260 for(size_t i = 0; i < n_right; i++) 261 f_right_[k].data()[i] = tx[ ( 3 + n_left + i) * n_order + k ]; 262 } 263 // ------------------------------------------------------------------- 264 // result for each order 265 // (we could avoid recalculting f_result_[k] for k=0,...,p-1) 266 for(size_t k = 0; k < n_order; k++) 267 { // result[k] = sum_ell left[ell] * right[k-ell] 268 f_result_[k] = matrix::Zero(nr_left, nc_right); 269 for(size_t ell = 0; ell <= k; ell++) 270 f_result_[k] += f_left_[ell] * f_right_[k-ell]; 271 } 272 // ------------------------------------------------------------------- 273 // pack result_ into ty 274 for(size_t k = 0; k < n_order; k++) 275 { for(size_t i = 0; i < n_result; i++) 276 ty[ i * n_order + k ] = f_result_[k].data()[i]; 277 } 278 // ------------------------------------------------------------------ 279 // check if we are computing vy 280 if( vx.size() == 0 ) 281 return true; 282 // ------------------------------------------------------------------ 283 // compute variable information for y; i.e., vy 284 // (note that the constant zero times a variable is a constant) 285 scalar zero(0.0); 286 assert( n_order == 1 ); 287 for(size_t i = 0; i < nr_left; i++) 288 { for(size_t j = 0; j < nc_right; j++) 289 { bool var = false; 290 for(size_t ell = 0; ell < n_middle; ell++) 291 { // left information 292 size_t index = 3 + i * n_middle + ell; 293 bool var_left = vx[index]; 294 bool nz_left = var_left | (f_left_[0](i, ell) != zero); 295 // right information 296 index = 3 + n_left + ell * nc_right + j; 297 bool var_right = vx[index]; 298 bool nz_right = var_right | (f_right_[0](ell, j) != zero); 299 // effect of result 300 var |= var_left & nz_right; 301 var |= nz_left & var_right; 302 } 303 size_t index = i * nc_right + j; 304 vy[index] = var; 305 } 306 } 307 return true; 308 } 309 /* %$$310 subhead reverse$$ 311$srccode%cpp% */
312  // reverse mode routine called by CppAD
313  virtual bool reverse(
314  // highest order Taylor coefficient that we are computing derivative of
315  size_t q ,
316  // forward mode Taylor coefficients for x variables
318  // forward mode Taylor coefficients for y variables
320  // upon return, derivative of G[ F[ {x_j^k} ] ] w.r.t {x_j^k}
322  // derivative of G[ {y_i^k} ] w.r.t. {y_i^k}
324  )
325  { size_t n_order = q + 1;
326  size_t nr_left = size_t( CppAD::Integer( tx[ 0 * n_order + 0 ] ) );
327  size_t n_middle = size_t( CppAD::Integer( tx[ 1 * n_order + 0 ] ) );
328  size_t nc_right = size_t( CppAD::Integer( tx[ 2 * n_order + 0 ] ) );
329 # ifndef NDEBUG
330  size_t nx = 3 + (nr_left + nc_right) * n_middle;
331  size_t ny = nr_left * nc_right;
332 # endif
333  //
334  assert( nx * n_order == tx.size() );
335  assert( ny * n_order == ty.size() );
336  assert( px.size() == tx.size() );
337  assert( py.size() == ty.size() );
338  //
339  size_t n_left = nr_left * n_middle;
340  size_t n_right = n_middle * nc_right;
341  size_t n_result = nr_left * nc_right;
342  assert( 3 + n_left + n_right == nx );
343  assert( n_result == ny );
344  // -------------------------------------------------------------------
345  // make sure f_left_, f_right_ are large enough
346  assert( f_left_.size() == f_right_.size() );
347  assert( f_left_.size() == f_result_.size() );
348  // must have previous run forward with order >= n_order
349  assert( f_left_.size() >= n_order );
350  // -------------------------------------------------------------------
351  // make sure r_left_, r_right_, and r_result_ are large enough
352  assert( r_left_.size() == r_right_.size() );
353  assert( r_left_.size() == r_result_.size() );
354  if( r_left_.size() < n_order )
355  { r_left_.resize(n_order);
356  r_right_.resize(n_order);
357  r_result_.resize(n_order);
358  //
359  for(size_t k = 0; k < n_order; k++)
360  { r_left_[k].resize(nr_left, n_middle);
361  r_right_[k].resize(n_middle, nc_right);
362  r_result_[k].resize(nr_left, nc_right);
363  }
364  }
365  // -------------------------------------------------------------------
366  // unpack tx into f_left and f_right
367  for(size_t k = 0; k < n_order; k++)
368  { // unpack left values for this order
369  for(size_t i = 0; i < n_left; i++)
370  f_left_[k].data()[i] = tx[ (3 + i) * n_order + k ];
371  //
372  // unpack right values for this order
373  for(size_t i = 0; i < n_right; i++)
374  f_right_[k].data()[i] = tx[ (3 + n_left + i) * n_order + k ];
375  }
376  // -------------------------------------------------------------------
377  // unpack py into r_result_
378  for(size_t k = 0; k < n_order; k++)
379  { for(size_t i = 0; i < n_result; i++)
380  r_result_[k].data()[i] = py[ i * n_order + k ];
381  }
382  // -------------------------------------------------------------------
383  // initialize r_left_ and r_right_ as zero
384  for(size_t k = 0; k < n_order; k++)
385  { r_left_[k] = matrix::Zero(nr_left, n_middle);
386  r_right_[k] = matrix::Zero(n_middle, nc_right);
387  }
388  // -------------------------------------------------------------------
389  // matrix reverse mode calculation
390  for(size_t k1 = n_order; k1 > 0; k1--)
391  { size_t k = k1 - 1;
392  for(size_t ell = 0; ell <= k; ell++)
393  { // nr x nm = nr x nc * nc * nm
394  r_left_[ell] += r_result_[k] * f_right_[k-ell].transpose();
395  // nm x nc = nm x nr * nr * nc
396  r_right_[k-ell] += f_left_[ell].transpose() * r_result_[k];
397  }
398  }
399  // -------------------------------------------------------------------
400  // pack r_left and r_right int px
401  for(size_t k = 0; k < n_order; k++)
402  { // dimensions are integer constants
403  px[ 0 * n_order + k ] = 0.0;
404  px[ 1 * n_order + k ] = 0.0;
405  px[ 2 * n_order + k ] = 0.0;
406  //
407  // pack left values for this order
408  for(size_t i = 0; i < n_left; i++)
409  px[ (3 + i) * n_order + k ] = r_left_[k].data()[i];
410  //
411  // pack right values for this order
412  for(size_t i = 0; i < n_right; i++)
413  px[ (3 + i + n_left) * n_order + k] = r_right_[k].data()[i];
414  }
415  //
416  return true;
417  }
418 /* %$$419 subhead for_sparse_jac$$
420 $srccode%cpp% */ 421 // forward Jacobian sparsity routine called by CppAD 422 virtual bool for_sparse_jac( 423 // number of columns in the matrix R 424 size_t q , 425 // sparsity pattern for the matrix R 426 const CppAD::vector< std::set<size_t> >& r , 427 // sparsity pattern for the matrix S = f'(x) * R 428 CppAD::vector< std::set<size_t> >& s , 429 const CppAD::vector<Base>& x ) 430 { 431 size_t nr_left = size_t( CppAD::Integer( x[0] ) ); 432 size_t n_middle = size_t( CppAD::Integer( x[1] ) ); 433 size_t nc_right = size_t( CppAD::Integer( x[2] ) ); 434 # ifndef NDEBUG 435 size_t nx = 3 + (nr_left + nc_right) * n_middle; 436 size_t ny = nr_left * nc_right; 437 # endif 438 // 439 assert( nx == r.size() ); 440 assert( ny == s.size() ); 441 // 442 size_t n_left = nr_left * n_middle; 443 for(size_t i = 0; i < nr_left; i++) 444 { for(size_t j = 0; j < nc_right; j++) 445 { // pack index for entry (i, j) in result 446 size_t i_result = i * nc_right + j; 447 s[i_result].clear(); 448 for(size_t ell = 0; ell < n_middle; ell++) 449 { // pack index for entry (i, ell) in left 450 size_t i_left = 3 + i * n_middle + ell; 451 // pack index for entry (ell, j) in right 452 size_t i_right = 3 + n_left + ell * nc_right + j; 453 // check if result of for this product is alwasy zero 454 // note that x is nan for commponents that are variables 455 bool zero = x[i_left] == Base(0.0) || x[i_right] == Base(0); 456 if( ! zero ) 457 { s[i_result] = 458 CppAD::set_union(s[i_result], r[i_left] ); 459 s[i_result] = 460 CppAD::set_union(s[i_result], r[i_right] ); 461 } 462 } 463 } 464 } 465 return true; 466 } 467 /* %$$468 subhead rev_sparse_jac$$ 469$srccode%cpp% */
470  // reverse Jacobian sparsity routine called by CppAD
471  virtual bool rev_sparse_jac(
472  // number of columns in the matrix R^T
473  size_t q ,
474  // sparsity pattern for the matrix R^T
475  const CppAD::vector< std::set<size_t> >& rt ,
476  // sparsoity pattern for the matrix S^T = f'(x)^T * R^T
477  CppAD::vector< std::set<size_t> >& st ,
479  {
480  size_t nr_left = size_t( CppAD::Integer( x[0] ) );
481  size_t n_middle = size_t( CppAD::Integer( x[1] ) );
482  size_t nc_right = size_t( CppAD::Integer( x[2] ) );
483  size_t nx = 3 + (nr_left + nc_right) * n_middle;
484 # ifndef NDEBUG
485  size_t ny = nr_left * nc_right;
486 # endif
487  //
488  assert( nx == st.size() );
489  assert( ny == rt.size() );
490  //
491  // initialize S^T as empty
492  for(size_t i = 0; i < nx; i++)
493  st[i].clear();
494
495  // sparsity for S(x)^T = f'(x)^T * R^T
496  size_t n_left = nr_left * n_middle;
497  for(size_t i = 0; i < nr_left; i++)
498  { for(size_t j = 0; j < nc_right; j++)
499  { // pack index for entry (i, j) in result
500  size_t i_result = i * nc_right + j;
501  st[i_result].clear();
502  for(size_t ell = 0; ell < n_middle; ell++)
503  { // pack index for entry (i, ell) in left
504  size_t i_left = 3 + i * n_middle + ell;
505  // pack index for entry (ell, j) in right
506  size_t i_right = 3 + n_left + ell * nc_right + j;
507  //
510  }
511  }
512  }
513  return true;
514  }
515 /* %$$516 subhead for_sparse_hes$$
517 $srccode%cpp% */ 518 virtual bool for_sparse_hes( 519 // which components of x are variables for this call 520 const CppAD::vector<bool>& vx, 521 // sparsity pattern for the diagonal of R 522 const CppAD::vector<bool>& r , 523 // sparsity pattern for the vector S 524 const CppAD::vector<bool>& s , 525 // sparsity patternfor the Hessian H(x) 526 CppAD::vector< std::set<size_t> >& h , 527 const CppAD::vector<Base>& x ) 528 { 529 size_t nr_left = size_t( CppAD::Integer( x[0] ) ); 530 size_t n_middle = size_t( CppAD::Integer( x[1] ) ); 531 size_t nc_right = size_t( CppAD::Integer( x[2] ) ); 532 size_t nx = 3 + (nr_left + nc_right) * n_middle; 533 # ifndef NDEBUG 534 size_t ny = nr_left * nc_right; 535 # endif 536 // 537 assert( vx.size() == nx ); 538 assert( r.size() == nx ); 539 assert( s.size() == ny ); 540 assert( h.size() == nx ); 541 // 542 // initilize h as empty 543 for(size_t i = 0; i < nx; i++) 544 h[i].clear(); 545 // 546 size_t n_left = nr_left * n_middle; 547 for(size_t i = 0; i < nr_left; i++) 548 { for(size_t j = 0; j < nc_right; j++) 549 { // pack index for entry (i, j) in result 550 size_t i_result = i * nc_right + j; 551 if( s[i_result] ) 552 { for(size_t ell = 0; ell < n_middle; ell++) 553 { // pack index for entry (i, ell) in left 554 size_t i_left = 3 + i * n_middle + ell; 555 // pack index for entry (ell, j) in right 556 size_t i_right = 3 + n_left + ell * nc_right + j; 557 if( r[i_left] & r[i_right] ) 558 { h[i_left].insert(i_right); 559 h[i_right].insert(i_left); 560 } 561 } 562 } 563 } 564 } 565 return true; 566 } 567 /* %$$568 subhead rev_sparse_hes$$ 569$srccode%cpp% */
570  // reverse Hessian sparsity routine called by CppAD
571  virtual bool rev_sparse_hes(
572  // which components of x are variables for this call
574  // sparsity pattern for S(x) = g'[f(x)]
576  // sparsity pattern for d/dx g[f(x)] = S(x) * f'(x)
578  // number of columns in R, U(x), and V(x)
579  size_t q ,
580  // sparsity pattern for R
581  const CppAD::vector< std::set<size_t> >& r ,
582  // sparsity pattern for U(x) = g^{(2)} [ f(x) ] * f'(x) * R
583  const CppAD::vector< std::set<size_t> >& u ,
584  // sparsity pattern for
585  // V(x) = f'(x)^T * U(x) + sum_{i=0}^{m-1} S_i(x) f_i^{(2)} (x) * R
586  CppAD::vector< std::set<size_t> >& v ,
587  // parameters as integers
589  {
590  size_t nr_left = size_t( CppAD::Integer( x[0] ) );
591  size_t n_middle = size_t( CppAD::Integer( x[1] ) );
592  size_t nc_right = size_t( CppAD::Integer( x[2] ) );
593  size_t nx = 3 + (nr_left + nc_right) * n_middle;
594 # ifndef NDEBUG
595  size_t ny = nr_left * nc_right;
596 # endif
597  //
598  assert( vx.size() == nx );
599  assert( s.size() == ny );
600  assert( t.size() == nx );
601  assert( r.size() == nx );
602  assert( v.size() == nx );
603  //
604  // initilaize return sparsity patterns as false
605  for(size_t j = 0; j < nx; j++)
606  { t[j] = false;
607  v[j].clear();
608  }
609  //
610  size_t n_left = nr_left * n_middle;
611  for(size_t i = 0; i < nr_left; i++)
612  { for(size_t j = 0; j < nc_right; j++)
613  { // pack index for entry (i, j) in result
614  size_t i_result = i * nc_right + j;
615  for(size_t ell = 0; ell < n_middle; ell++)
616  { // pack index for entry (i, ell) in left
617  size_t i_left = 3 + i * n_middle + ell;
618  // pack index for entry (ell, j) in right
619  size_t i_right = 3 + n_left + ell * nc_right + j;
620  //
621  // back propagate T(x) = S(x) * f'(x).
622  t[i_left] |= bool( s[i_result] );
623  t[i_right] |= bool( s[i_result] );
624  //
625  // V(x) = f'(x)^T * U(x) + sum_i S_i(x) * f_i''(x) * R
626  // U(x) = g''[ f(x) ] * f'(x) * R
627  // S_i(x) = g_i'[ f(x) ]
628  //
629  // back propagate f'(x)^T * U(x)
630  v[i_left] = CppAD::set_union(v[i_left], u[i_result] );
631  v[i_right] = CppAD::set_union(v[i_right], u[i_result] );
632  //
633  // back propagate S_i(x) * f_i''(x) * R
634  // (here is where we use vx to check for cross terms)
635  if( s[i_result] & vx[i_left] & vx[i_right] )
636  { v[i_left] = CppAD::set_union(v[i_left], r[i_right] );
637  v[i_right] = CppAD::set_union(v[i_right], r[i_left] );
638  }
639  }
640  }
641  }
642  return true;
643  }
644 /* %$$645 head End Class Definition$$
646 $srccode%cpp% */ 647 }; // End of atomic_eigen_mat_mul class 648 649 } // END_EMPTY_NAMESPACE 650 /* %$$651$$$comment end nospell
652 \$end
653 */
654
655
656 # endif
includes the entire CppAD package in the necessary order.
std::set< Element > set_union(const std::set< Element > &left, const std::set< Element > &right)
Definition: set_union.hpp:73
void clear(void)
free memory and set number of elements to zero
Definition: vector.hpp:418
Type * data(void)
raw pointer to the data
Definition: vector.hpp:391