CppAD: A C++ Algorithmic Differentiation Package  20171217
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
eigen_mat_mul.hpp
Go to the documentation of this file.
1 # ifndef CPPAD_EXAMPLE_EIGEN_MAT_MUL_HPP
2 # define CPPAD_EXAMPLE_EIGEN_MAT_MUL_HPP
3 
4 /* --------------------------------------------------------------------------
5 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-17 Bradley M. Bell
6 
7 CppAD is distributed under multiple licenses. This distribution is under
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.
12 Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
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>
121 class atomic_eigen_mat_mul : public CppAD::atomic_base<Base> {
122 public:
123  // -----------------------------------------------------------
124  // type of elements during calculation of derivatives
125  typedef Base scalar;
126  // type of elements during taping
127  typedef CppAD::AD<scalar> ad_scalar;
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
147  ad_matrix op(
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
195  CppAD::vector<matrix> f_left_, f_right_, f_result_;
196  // one reverse mode vector of matrices for left, right, and result
197  CppAD::vector<matrix> r_left_, r_right_, r_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
215  CppAD::vector<scalar>& ty
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
317  const CppAD::vector<double>& tx ,
318  // forward mode Taylor coefficients for y variables
319  const CppAD::vector<double>& ty ,
320  // upon return, derivative of G[ F[ {x_j^k} ] ] w.r.t {x_j^k}
321  CppAD::vector<double>& px ,
322  // derivative of G[ {y_i^k} ] w.r.t. {y_i^k}
323  const CppAD::vector<double>& py
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 ,
478  const CppAD::vector<Base>& x )
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  //
508  st[i_left] = CppAD::set_union(st[i_left], rt[i_result]);
509  st[i_right] = CppAD::set_union(st[i_right], rt[i_result]);
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
573  const CppAD::vector<bool>& vx,
574  // sparsity pattern for S(x) = g'[f(x)]
575  const CppAD::vector<bool>& s ,
576  // sparsity pattern for d/dx g[f(x)] = S(x) * f'(x)
577  CppAD::vector<bool>& t ,
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
588  const CppAD::vector<Base>& x )
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
cppad.hpp
includes the entire CppAD package in the necessary order.
CppAD::AD
Definition: ad.hpp:34
CppAD::set_union
std::set< Element > set_union(const std::set< Element > &left, const std::set< Element > &right)
Definition: set_union.hpp:73
CppAD::vector::clear
void clear(void)
free memory and set number of elements to zero
Definition: vector.hpp:418
CppAD::vector::data
Type * data(void)
raw pointer to the data
Definition: vector.hpp:391
anonymous_namespace{eigen_mat_mul.hpp}::atomic_eigen_mat_mul::atomic_eigen_mat_mul
atomic_eigen_mat_mul(void)
Definition: eigen_mat_mul.hpp:138
anonymous_namespace{eigen_mat_mul.hpp}::atomic_eigen_mat_mul::op
ad_matrix op(const ad_matrix &left, const ad_matrix &right)
Definition: eigen_mat_mul.hpp:147
CppAD::vector::size
size_t size(void) const
number of elements currently in this vector.
Definition: vector.hpp:387
anonymous_namespace{eigen_mat_mul.hpp}::atomic_eigen_mat_mul::r_right_
CppAD::vector< matrix > r_right_
Definition: eigen_mat_mul.hpp:197
anonymous_namespace{eigen_mat_mul.hpp}::atomic_eigen_mat_mul
Definition: eigen_mat_mul.hpp:121
anonymous_namespace{eigen_mat_mul.hpp}::atomic_eigen_mat_mul::forward
virtual bool forward(size_t p, size_t q, const CppAD::vector< bool > &vx, CppAD::vector< bool > &vy, const CppAD::vector< scalar > &tx, CppAD::vector< scalar > &ty)
Definition: eigen_mat_mul.hpp:203
anonymous_namespace{eigen_mat_mul.hpp}::atomic_eigen_mat_mul::matrix
Eigen::Matrix< scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > matrix
Definition: eigen_mat_mul.hpp:130
anonymous_namespace{eigen_mat_mul.hpp}::atomic_eigen_mat_mul::ad_matrix
Eigen::Matrix< ad_scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > ad_matrix
Definition: eigen_mat_mul.hpp:133
anonymous_namespace{eigen_mat_mul.hpp}::atomic_eigen_mat_mul::for_sparse_hes
virtual bool for_sparse_hes(const CppAD::vector< bool > &vx, const CppAD::vector< bool > &r, const CppAD::vector< bool > &s, CppAD::vector< std::set< size_t > > &h, const CppAD::vector< Base > &x)
Definition: eigen_mat_mul.hpp:518
anonymous_namespace{eigen_mat_mul.hpp}::atomic_eigen_mat_mul::rev_sparse_jac
virtual bool rev_sparse_jac(size_t q, const CppAD::vector< std::set< size_t > > &rt, CppAD::vector< std::set< size_t > > &st, const CppAD::vector< Base > &x)
Definition: eigen_mat_mul.hpp:471
CppAD::atomic_base
Definition: atomic_base.hpp:28
CppAD::Integer
int Integer(const std::complex< double > &x)
Definition: base_complex.hpp:175
anonymous_namespace{eigen_mat_mul.hpp}::atomic_eigen_mat_mul::for_sparse_jac
virtual bool for_sparse_jac(size_t q, const CppAD::vector< std::set< size_t > > &r, CppAD::vector< std::set< size_t > > &s, const CppAD::vector< Base > &x)
Definition: eigen_mat_mul.hpp:422
anonymous_namespace{eigen_mat_mul.hpp}::atomic_eigen_mat_mul::reverse
virtual bool reverse(size_t q, const CppAD::vector< double > &tx, const CppAD::vector< double > &ty, CppAD::vector< double > &px, const CppAD::vector< double > &py)
Definition: eigen_mat_mul.hpp:313
anonymous_namespace{eigen_mat_mul.hpp}::atomic_eigen_mat_mul::rev_sparse_hes
virtual bool rev_sparse_hes(const CppAD::vector< bool > &vx, const CppAD::vector< bool > &s, CppAD::vector< bool > &t, size_t q, const CppAD::vector< std::set< size_t > > &r, const CppAD::vector< std::set< size_t > > &u, CppAD::vector< std::set< size_t > > &v, const CppAD::vector< Base > &x)
Definition: eigen_mat_mul.hpp:571
anonymous_namespace{eigen_mat_mul.hpp}::atomic_eigen_mat_mul::ad_scalar
CppAD::AD< scalar > ad_scalar
Definition: eigen_mat_mul.hpp:127
anonymous_namespace{eigen_mat_mul.hpp}::atomic_eigen_mat_mul::scalar
Base scalar
Definition: eigen_mat_mul.hpp:125
anonymous_namespace{eigen_mat_mul.hpp}::atomic_eigen_mat_mul::f_right_
CppAD::vector< matrix > f_right_
Definition: eigen_mat_mul.hpp:195