source: trunk/cppad/example/eigen_mat_mul.hpp @ 3810

Last change on this file since 3810 was 3810, checked in by bradbell, 4 years ago

merge to branch: trunk
from repository: https://github.com/coin-or/CppAD
start hash code: dd967ef41b8d6731d90ebb8b3e7e8b863565289c
end hash code: f77cd2e443b8d14f354f97efc4c26369227dfe81

commit f77cd2e443b8d14f354f97efc4c26369227dfe81
Author: Brad Bell <bradbell@…>
Date: Sat Mar 26 06:47:33 2016 -0700

  1. Add reverse mode to atomic_eigen_mat_inv.
  2. Advance version to cppad-20160326.


eigen_mat_mul.cpp: impove reverse mode comments.

commit 94c81399191d06ae9cef81c271f520b0e78f079d
Author: Brad Bell <bradbell@…>
Date: Fri Mar 25 12:16:03 2016 -0700

whats_new_16.omh: add comment about changes to eigen_mat_mul.

commit d647288c9e022f16ef94229db612a774808a5de3
Author: Brad Bell <bradbell@…>
Date: Fri Mar 25 11:59:45 2016 -0700

Replace pack and unpack by member function that inverts matrices.

commit 6fba95d90e50590abd438946f47025e0a6e34c2c
Author: Brad Bell <bradbell@…>
Date: Fri Mar 25 11:45:43 2016 -0700

Replace pack and unpack by member function that multiplies matrices.

commit cac97c115abad9fba0d0c3706f4bc0e5a7f617f0
Author: Brad Bell <bradbell@…>
Date: Fri Mar 25 11:11:16 2016 -0700

Change atomic_eigen_mat_div -> atomic_eigen_mat_inv.

commit 1f7e52a4d96986f7be2ad2d6abeafc086c4f7622
Author: Brad Bell <bradbell@…>
Date: Fri Mar 25 08:09:39 2016 -0700

  1. Implement and test forward for atomic_eigen_mat_div.
  2. Advance version to 20160325.


eigen_mat_mul.hpp: edit comments and documentation.

commit e953f669d8e9e122dc2f182c28c051c9ba75c705
Author: Brad Bell <bradbell@…>
Date: Fri Mar 25 04:32:45 2016 -0700

eigen_mat_mul.cpp: No longer under construction.

commit 3a62044eaebb6393676276e5b843fbe5db4f45a2
Author: Brad Bell <bradbell@…>
Date: Fri Mar 25 04:07:18 2016 -0700

eigen_mat_mul.hpp: add assert for f_left_ and f_right_ in reverse.

commit 721d64783547487fe69100c7190628efab7149c3
Author: Brad Bell <bradbell@…>
Date: Fri Mar 25 03:59:05 2016 -0700

eigen_mat_mul.hpp: typo Publice -> Public.

File size: 18.6 KB
Line 
1// $Id$
2# ifndef CPPAD_EXAMPLE_EIGEN_MAT_MUL_HPP
3# define CPPAD_EXAMPLE_EIGEN_MAT_MUL_HPP
4
5/* --------------------------------------------------------------------------
6CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-16 Bradley M. Bell
7
8CppAD is distributed under multiple licenses. This distribution is under
9the terms of the
10                    Eclipse Public License Version 1.0.
11
12A copy of this license is included in the COPYING file of this distribution.
13Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
14-------------------------------------------------------------------------- */
15
16/*
17$begin atomic_eigen_mat_mul.hpp$$
18$spell
19        Eigen
20        Taylor
21$$
22
23$section Atomic Eigen Matrix Multiply Class$$
24
25$head Purpose$$
26For fixed positive integers $latex r$$, $latex m$$, $latex c$$,
27construct and atomic operation that computes the matrix product
28$latex R = A \times B$$ for any $latex A \in \B{R}^{r \times m}$$ and
29$latex B \in \B{R}^{m \times c}$$.
30
31$head Theory$$
32
33$subhead Forward$$
34For $latex k = 0 , \ldots $$, the $th k$$ order Taylor coefficient
35$latex R_k$$ is given by
36$latex \[
37        R_k = \sum_{\ell = 0}^{k} A_\ell B_{k-\ell}
38\] $$
39
40$subhead Reverse$$
41We use $latex \bar{R}_k$$ for the partial of the scalar final result
42with respect to $latex R_k$$.
43The back-propagation algorithm that eliminates $latex R_k$$ is,
44for $latex \ell = 0, \ldots , k-1$$,
45$latex \[
46\bar{A}_\ell     = \bar{A}_\ell     + \bar{R}_k B_{k-\ell}^\R{T}
47\] $$
48$latex \[
49\bar{B}_{k-\ell} =  \bar{B}_{k-\ell} + A_\ell^\R{T} \bar{R}_k
50\] $$
51
52
53$nospell
54
55$head Start Class Definition$$
56$srccode%cpp% */
57# include <cppad/cppad.hpp>
58# include <Eigen/Core>
59
60
61/* %$$
62$head Public$$
63
64$subhead Types$$
65$srccode%cpp% */
66namespace { // BEGIN_EMPTY_NAMESPACE
67
68template <class Base>
69class atomic_eigen_mat_mul : public CppAD::atomic_base<Base> {
70public:
71        // -----------------------------------------------------------
72        // type of elements during calculation of derivatives
73        typedef Base              scalar;
74        // type of elements during taping
75        typedef CppAD::AD<scalar> ad_scalar;
76        // type of matrix during calculation of derivatives
77        typedef Eigen::Matrix<
78                scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>     matrix;
79        // type of matrix during taping
80        typedef Eigen::Matrix<
81                ad_scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > ad_matrix;
82/* %$$
83$subhead Constructor$$
84$srccode%cpp% */
85        // constructor
86        atomic_eigen_mat_mul(
87                // number of rows in left operand
88                const size_t nr_left  ,
89                // number of rows in left and columns in right operand
90                const size_t n_middle   ,
91                // number of columns in right operand
92                const size_t nc_right
93        ) :
94        CppAD::atomic_base<Base>(
95                "atom_eigen_mat_mul"                             ,
96                CppAD::atomic_base<Base>::set_sparsity_enum
97        ) ,
98        nr_left_(  nr_left  )                   ,
99        n_middle_(   n_middle   )                   ,
100        nc_right_( nc_right )                   ,
101        nx_( (nr_left + nc_right) * n_middle_ ) ,
102        ny_( nr_left * nc_right )
103        { }
104/* %$$
105$subhead op$$
106$srccode%cpp% */
107        // use atomic operation to multiply two AD matrices
108        ad_matrix op(
109                const ad_matrix&              left    ,
110                const ad_matrix&              right   )
111        {       assert( nr_left_   == size_t( left.rows() )   );
112                assert( n_middle_  == size_t( left.cols() )   );
113                assert( n_middle_  == size_t( right.rows() )  );
114                assert( nc_right_  == size_t ( right.cols() ) );
115
116                // -----------------------------------------------------------------
117                // packed version of left and right
118                CPPAD_TESTVECTOR(ad_scalar) packed_arg(nx_);
119                size_t n_left   = nr_left_ * n_middle_;
120                size_t n_right  = n_middle_ * nc_right_;
121                size_t n_result = nr_left_ * nc_right_;
122                assert( n_left + n_right == nx_ );
123                assert( n_result == ny_ );
124                //
125                for(size_t i = 0; i < n_left; i++)
126                        packed_arg[i] = left.data()[i];
127                for(size_t i = 0; i < n_right; i++)
128                        packed_arg[ i + n_left ] = right.data()[i];
129
130                // ------------------------------------------------------------------
131                // packed version of result = left * right
132                CPPAD_TESTVECTOR(ad_scalar) packed_result(ny_);
133                (*this)(packed_arg, packed_result);
134
135                // ------------------------------------------------------------------
136                // unpack result matrix
137                ad_matrix result(nr_left_, nc_right_);
138                for(size_t i = 0; i < n_result; i++)
139                        result.data()[i] = packed_result[ i ];
140                //
141                return result;
142        }
143/* %$$
144$head Private$$
145
146$subhead Variables$$
147$srccode%cpp% */
148private:
149        // -------------------------------------------------------------
150        // number of of rows in left operand and result
151        const size_t nr_left_;
152        // number of of columns in left operand and rows in right operand
153        const size_t n_middle_;
154        // number of columns in right operand and result
155        const size_t nc_right_;
156        // size of the domain space
157        const size_t nx_;
158        // size of the range space
159        const size_t ny_;
160        // -------------------------------------------------------------
161        // one forward mode vector of matrices for left, right, and result
162        CppAD::vector<matrix> f_left_, f_right_, f_result_;
163        // one reverse mode vector of matrices for left, right, and result
164        CppAD::vector<matrix> r_left_, r_right_, r_result_;
165        // -------------------------------------------------------------
166/* %$$
167$subhead rows$$
168$srccode%cpp% */
169        // convert from int to size_t
170        static size_t rows(const matrix& x)
171        {       return size_t( x.rows() ); }
172        static size_t rows(const ad_matrix& x)
173        {       return size_t( x.rows() ); }
174/* %$$
175$subhead cols$$
176$srccode%cpp% */
177        // convert from int to size_t
178        static size_t cols(const matrix& x)
179        {       return size_t( x.cols() ); }
180        static size_t cols(const ad_matrix& x)
181        {       return size_t( x.cols() ); }
182/* %$$
183$subhead forward$$
184$srccode%cpp% */
185        // forward mode routine called by CppAD
186        virtual bool forward(
187                // lowest order Taylor coefficient we are evaluating
188                size_t                          p ,
189                // highest order Taylor coefficient we are evaluating
190                size_t                          q ,
191                // which components of x are variables
192                const CppAD::vector<bool>&      vx ,
193                // which components of y are variables
194                CppAD::vector<bool>&            vy ,
195                // tx [ j * (q+1) + k ] is x_j^k
196                const CppAD::vector<scalar>&    tx ,
197                // ty [ i * (q+1) + k ] is y_i^k
198                CppAD::vector<scalar>&          ty
199        )
200        {       size_t n_order = q + 1;
201                assert( vx.size() == 0 || nx_ == vx.size() );
202                assert( vx.size() == 0 || ny_ == vy.size() );
203                assert( nx_ * n_order == tx.size() );
204                assert( ny_ * n_order == ty.size() );
205                //
206                size_t n_left   = nr_left_ * n_middle_;
207                size_t n_right  = n_middle_ * nc_right_;
208                size_t n_result = nr_left_ * nc_right_;
209                assert( n_left + n_right == nx_ );
210                assert( n_result == ny_ );
211                //
212                // -------------------------------------------------------------------
213                // make sure f_left_, f_right_, and f_result_ are large enough
214                assert( f_left_.size() == f_right_.size() );
215                assert( f_left_.size() == f_result_.size() );
216                if( f_left_.size() < n_order )
217                {       f_left_.resize(n_order);
218                        f_right_.resize(n_order);
219                        f_result_.resize(n_order);
220                        //
221                        for(size_t k = 0; k < n_order; k++)
222                        {       f_left_[k].resize(nr_left_, n_middle_);
223                                f_right_[k].resize(n_middle_, nc_right_);
224                                f_result_[k].resize(nr_left_, nc_right_);
225                        }
226                }
227                // -------------------------------------------------------------------
228                // unpack tx into f_left and f_right
229                for(size_t k = 0; k < n_order; k++)
230                {       // unpack left values for this order
231                        for(size_t i = 0; i < n_left; i++)
232                                f_left_[k].data()[i] = tx[ i * n_order + k ];
233                        //
234                        // unpack right values for this order
235                        for(size_t i = 0; i < n_right; i++)
236                                f_right_[k].data()[i] = tx[ (i + n_left) * n_order + k ];
237                }
238                // -------------------------------------------------------------------
239                // result for each order
240                for(size_t k = 0; k < n_order; k++)
241                {       // result[k] = sum_ell left[ell] * right[k-ell]
242                        f_result_[k] = matrix::Zero(nr_left_, nc_right_);
243                        for(size_t ell = 0; ell <= k; ell++)
244                                f_result_[k] += f_left_[ell] * f_right_[k-ell];
245                }
246                // -------------------------------------------------------------------
247                // pack result_ into ty
248                for(size_t k = 0; k < n_order; k++)
249                {       for(size_t i = 0; i < n_result; i++)
250                                ty[ i * n_order + k ] = f_result_[k].data()[i];
251                }
252
253                // check if we are computing vy
254                if( vx.size() == 0 )
255                        return true;
256                // ------------------------------------------------------------------
257                // compute variable information for y; i.e., vy
258                // (note that the constant zero times a variable is a constant)
259                scalar zero(0.0);
260                assert( n_order == 1 );
261                for(size_t i = 0; i < nr_left_; i++)
262                {       for(size_t j = 0; j < nc_right_; j++)
263                        {       bool var = false;
264                                for(size_t ell = 0; ell < n_middle_; ell++)
265                                {       // left information
266                                        size_t index   = i * n_middle_ + ell;
267                                        bool var_left  = vx[index];
268                                        bool nz_left   = var_left | (f_left_[0](i, ell) != zero);
269                                        // right information
270                                        index          = n_left + ell * nc_right_ + j;
271                                        bool var_right = vx[index];
272                                        bool nz_right  = var_right | (f_right_[0](ell, j) != zero);
273                                        // effect of result
274                                        var |= var_left & nz_right;
275                                        var |= nz_left  & var_right;
276                                }
277                                size_t index = i * nc_right_ + j;
278                                vy[index]    = var;
279                        }
280                }
281                return true;
282        }
283/* %$$
284$subhead reverse$$
285$srccode%cpp% */
286        // reverse mode routine called by CppAD
287        virtual bool reverse(
288                // highest order Taylor coefficient that we are computing derivative of
289                size_t                     q ,
290                // forward mode Taylor coefficients for x variables
291                const CppAD::vector<double>&     tx ,
292                // forward mode Taylor coefficients for y variables
293                const CppAD::vector<double>&     ty ,
294                // upon return, derivative of G[ F[ {x_j^k} ] ] w.r.t {x_j^k}
295                CppAD::vector<double>&           px ,
296                // derivative of G[ {y_i^k} ] w.r.t. {y_i^k}
297                const CppAD::vector<double>&     py
298        )
299        {       size_t n_order = q + 1;
300                assert( nx_ * n_order == tx.size() );
301                assert( ny_ * n_order == ty.size() );
302                assert( px.size() == tx.size() );
303                assert( py.size() == ty.size() );
304                //
305                size_t n_left   = nr_left_ * n_middle_;
306                size_t n_right  = n_middle_ * nc_right_;
307                size_t n_result = nr_left_ * nc_right_;
308                assert( n_left + n_right == nx_ );
309                assert( n_result == ny_ );
310                // -------------------------------------------------------------------
311                // make sure f_left_, f_right_ are large enough
312                assert( f_left_.size() == f_right_.size() );
313                assert( f_left_.size() == f_result_.size() );
314                // must have previous run forward with order >= n_order
315                assert( f_left_.size() >= n_order );
316                // -------------------------------------------------------------------
317                // make sure r_left_, r_right_, and r_result_ are large enough
318                assert( r_left_.size() == r_right_.size() );
319                assert( r_left_.size() == r_result_.size() );
320                if( r_left_.size() < n_order )
321                {       r_left_.resize(n_order);
322                        r_right_.resize(n_order);
323                        r_result_.resize(n_order);
324                        //
325                        for(size_t k = 0; k < n_order; k++)
326                        {       r_left_[k].resize(nr_left_, n_middle_);
327                                r_right_[k].resize(n_middle_, nc_right_);
328                                r_result_[k].resize(nr_left_, nc_right_);
329                        }
330                }
331                // -------------------------------------------------------------------
332                // unpack tx into f_left and f_right
333                for(size_t k = 0; k < n_order; k++)
334                {       // unpack left values for this order
335                        for(size_t i = 0; i < n_left; i++)
336                                f_left_[k].data()[i] = tx[ i * n_order + k ];
337                        //
338                        // unpack right values for this order
339                        for(size_t i = 0; i < n_right; i++)
340                                f_right_[k].data()[i] = tx[ (i + n_left) * n_order + k ];
341                }
342                // -------------------------------------------------------------------
343                // unpack py into r_result_
344                for(size_t k = 0; k < n_order; k++)
345                {       for(size_t i = 0; i < n_result; i++)
346                                r_result_[k].data()[i] = py[ i * n_order + k ];
347                }
348                // -------------------------------------------------------------------
349                // initialize r_left_ and r_right_ as zero
350                for(size_t k = 0; k < n_order; k++)
351                {       r_left_[k]   = matrix::Zero(nr_left_, n_middle_);
352                        r_right_[k]  = matrix::Zero(n_middle_, nc_right_);
353                }
354                // -------------------------------------------------------------------
355                // matrix reverse mode calculation
356                for(size_t k1 = n_order; k1 > 0; k1--)
357                {       size_t k = k1 - 1;
358                        for(size_t ell = 0; ell <= k; ell++)
359                        {       // nr x nm       = nr x nc      * nc * nm
360                                r_left_[ell]    += r_result_[k] * f_right_[k-ell].transpose();
361                                // nm x nc       = nm x nr * nr * nc
362                                r_right_[k-ell] += f_left_[ell].transpose() * r_result_[k];
363                        }
364                }
365                // -------------------------------------------------------------------
366                // pack r_left and r_right int px
367                for(size_t k = 0; k < n_order; k++)
368                {       // pack left values for this order
369                        for(size_t i = 0; i < n_left; i++)
370                                px[ i * n_order + k ] = r_left_[k].data()[i];
371                        //
372                        // pack right values for this order
373                        for(size_t i = 0; i < n_right; i++)
374                                px[ (i + n_left) * n_order + k] = r_right_[k].data()[i];
375                }
376                //
377                return true;
378        }
379/* %$$
380$subhead for_sparse_jac$$
381$srccode%cpp% */
382        // forward Jacobian sparsity routine called by CppAD
383        virtual bool for_sparse_jac(
384                // number of columns in the matrix R
385                size_t                                       q ,
386                // sparsity pattern for the matrix R
387                const CppAD::vector< std::set<size_t> >&     r ,
388                // sparsity pattern for the matrix S = f'(x) * R
389                CppAD::vector< std::set<size_t> >&           s )
390        {       assert( nx_ == r.size() );
391                assert( ny_ == s.size() );
392                //
393                size_t n_left = nr_left_ * n_middle_;
394                for(size_t i = 0; i < nr_left_; i++)
395                {       for(size_t j = 0; j < nc_right_; j++)
396                        {       // pack index for entry (i, j) in result
397                                size_t i_result = i * nc_right_ + j;
398                                s[i_result].clear();
399                                for(size_t ell = 0; ell < n_middle_; ell++)
400                                {       // pack index for entry (i, ell) in left
401                                        size_t i_left  = i * n_middle_ + ell;
402                                        // pack index for entry (ell, j) in right
403                                        size_t i_right = n_left + ell * nc_right_ + j;
404                                        //
405                                        s[i_result] = CppAD::set_union(s[i_result], r[i_left] );
406                                        s[i_result] = CppAD::set_union(s[i_result], r[i_right] );
407                                }
408                        }
409                }
410                return true;
411        }
412/* %$$
413$subhead rev_sparse_jac$$
414$srccode%cpp% */
415        // reverse Jacobian sparsity routine called by CppAD
416        virtual bool rev_sparse_jac(
417                // number of columns in the matrix R^T
418                size_t                                      q  ,
419                // sparsity pattern for the matrix R^T
420                const CppAD::vector< std::set<size_t> >&    rt ,
421                // sparsoity pattern for the matrix S^T = f'(x)^T * R^T
422                CppAD::vector< std::set<size_t> >&          st )
423        {       assert( nx_ == st.size() );
424                assert( ny_ == rt.size() );
425
426                // initialize S^T as empty
427                for(size_t i = 0; i < nx_; i++)
428                        st[i].clear();
429
430                // sparsity for S(x)^T = f'(x)^T * R^T
431                size_t n_left = nr_left_ * n_middle_;
432                for(size_t i = 0; i < nr_left_; i++)
433                {       for(size_t j = 0; j < nc_right_; j++)
434                        {       // pack index for entry (i, j) in result
435                                size_t i_result = i * nc_right_ + j;
436                                st[i_result].clear();
437                                for(size_t ell = 0; ell < n_middle_; ell++)
438                                {       // pack index for entry (i, ell) in left
439                                        size_t i_left  = i * n_middle_ + ell;
440                                        // pack index for entry (ell, j) in right
441                                        size_t i_right = n_left + ell * nc_right_ + j;
442                                        //
443                                        st[i_left]  = CppAD::set_union(st[i_left],  rt[i_result]);
444                                        st[i_right] = CppAD::set_union(st[i_right], rt[i_result]);
445                                }
446                        }
447                }
448                return true;
449        }
450/* %$$
451$subhead for_sparse_hes$$
452$srccode%cpp% */
453        virtual bool for_sparse_hes(
454                // which components of x are variables for this call
455                const CppAD::vector<bool>&                   vx,
456                // sparsity pattern for the diagonal of R
457                const CppAD::vector<bool>&                   r ,
458                // sparsity pattern for the vector S
459                const CppAD::vector<bool>&                   s ,
460                // sparsity patternfor the Hessian H(x)
461                CppAD::vector< std::set<size_t> >&           h )
462        {       assert( vx.size() == nx_ );
463                assert( r.size()  == nx_ );
464                assert( s.size()  == ny_ );
465                assert( h.size()  == nx_ );
466                //
467                // initilize h as empty
468                for(size_t i = 0; i < nx_; i++)
469                        h[i].clear();
470                //
471                size_t n_left = nr_left_ * n_middle_;
472                for(size_t i = 0; i < nr_left_; i++)
473                {       for(size_t j = 0; j < nc_right_; j++)
474                        {       // pack index for entry (i, j) in result
475                                size_t i_result = i * nc_right_ + j;
476                                if( s[i_result] )
477                                {       for(size_t ell = 0; ell < n_middle_; ell++)
478                                        {       // pack index for entry (i, ell) in left
479                                                size_t i_left  = i * n_middle_ + ell;
480                                                // pack index for entry (ell, j) in right
481                                                size_t i_right = n_left + ell * nc_right_ + j;
482                                                if( r[i_left] & r[i_right] )
483                                                {       h[i_left].insert(i_right);
484                                                        h[i_right].insert(i_left);
485                                                }
486                                        }
487                                }
488                        }
489                }
490                return true;
491        }
492/* %$$
493$subhead rev_sparse_hes$$
494$srccode%cpp% */
495        // reverse Hessian sparsity routine called by CppAD
496        virtual bool rev_sparse_hes(
497                // which components of x are variables for this call
498                const CppAD::vector<bool>&                   vx,
499                // sparsity pattern for S(x) = g'[f(x)]
500                const CppAD::vector<bool>&                   s ,
501                // sparsity pattern for d/dx g[f(x)] = S(x) * f'(x)
502                CppAD::vector<bool>&                         t ,
503                // number of columns in R, U(x), and V(x)
504                size_t                                       q ,
505                // sparsity pattern for R
506                const CppAD::vector< std::set<size_t> >&     r ,
507                // sparsity pattern for U(x) = g^{(2)} [ f(x) ] * f'(x) * R
508                const CppAD::vector< std::set<size_t> >&     u ,
509                // sparsity pattern for
510                // V(x) = f'(x)^T * U(x) + sum_{i=0}^{m-1} S_i(x) f_i^{(2)} (x) * R
511                CppAD::vector< std::set<size_t> >&           v )
512        {       assert( vx.size() == nx_ );
513                assert( s.size()  == ny_ );
514                assert( t.size()  == nx_ );
515                assert( r.size()  == nx_ );
516                assert( v.size()  == nx_ );
517                //
518                // initilaize return sparsity patterns as false
519                for(size_t j = 0; j < nx_; j++)
520                {       t[j] = false;
521                        v[j].clear();
522                }
523                //
524                size_t n_left = nr_left_ * n_middle_;
525                for(size_t i = 0; i < nr_left_; i++)
526                {       for(size_t j = 0; j < nc_right_; j++)
527                        {       // pack index for entry (i, j) in result
528                                size_t i_result = i * nc_right_ + j;
529                                for(size_t ell = 0; ell < n_middle_; ell++)
530                                {       // pack index for entry (i, ell) in left
531                                        size_t i_left  = i * n_middle_ + ell;
532                                        // pack index for entry (ell, j) in right
533                                        size_t i_right = n_left + ell * nc_right_ + j;
534                                        //
535                                        // back propagate T(x) = S(x) * f'(x).
536                                        t[i_left]  |= bool( s[i_result] );
537                                        t[i_right] |= bool( s[i_result] );
538                                        //
539                                        // V(x) = f'(x)^T * U(x) +  sum_i S_i(x) * f_i''(x) * R
540                                        // U(x)   = g''[ f(x) ] * f'(x) * R
541                                        // S_i(x) = g_i'[ f(x) ]
542                                        //
543                                        // back propagate f'(x)^T * U(x)
544                                        v[i_left]  = CppAD::set_union(v[i_left],  u[i_result] );
545                                        v[i_right] = CppAD::set_union(v[i_right], u[i_result] );
546                                        //
547                                        // back propagate S_i(x) * f_i''(x) * R
548                                        // (here is where we use vx to check for cross terms)
549                                        if( s[i_result] & vx[i_left] & vx[i_right] )
550                                        {       v[i_left]  = CppAD::set_union(v[i_left],  r[i_right] );
551                                                v[i_right] = CppAD::set_union(v[i_right], r[i_left]  );
552                                        }
553                                }
554                        }
555                }
556                return true;
557        }
558/* %$$
559$head End Class Definition$$
560$srccode%cpp% */
561}; // End of atomic_eigen_mat_mul class
562
563}  // END_EMPTY_NAMESPACE
564/* %$$
565$$ $comment end nospell$$
566$end
567*/
568
569
570# endif
Note: See TracBrowser for help on using the repository browser.